isl_coalesce.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  * Copyright 2010      INRIA Saclay
   4  *
   5  * Use of this software is governed by the GNU LGPLv2.1 license
   6  *
   7  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   8  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   9  * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
  10  * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
  11  */
  12
  13 #include "isl_map_private.h"
  14 #include "isl_seq.h"
  15 #include "isl_tab.h"
  16
  17 #define STATUS_ERROR            -1
  18 #define STATUS_REDUNDANT         1
  19 #define STATUS_VALID             2
  20 #define STATUS_SEPARATE          3
  21 #define STATUS_CUT               4
  22 #define STATUS_ADJ_EQ            5
  23 #define STATUS_ADJ_INEQ          6
  24
  25 static int status_in(isl_int *ineq, struct isl_tab *tab)
  26 {
  27         enum isl_ineq_type type = isl_tab_ineq_type(tab, ineq);
  28         switch (type) {
  29         case isl_ineq_error:            return STATUS_ERROR;
  30         case isl_ineq_redundant:        return STATUS_VALID;
  31         case isl_ineq_separate:         return STATUS_SEPARATE;
  32         case isl_ineq_cut:              return STATUS_CUT;
  33         case isl_ineq_adj_eq:           return STATUS_ADJ_EQ;
  34         case isl_ineq_adj_ineq:         return STATUS_ADJ_INEQ;
  35         }
  36 }
  37
  38 /* Compute the position of the equalities of basic map "i"
  39  * with respect to basic map "j".
  40  * The resulting array has twice as many entries as the number
  41  * of equalities corresponding to the two inequalties to which
  42  * each equality corresponds.
  43  */
  44 static int *eq_status_in(struct isl_map *map, int i, int j,
  45         struct isl_tab **tabs)
  46 {
  47         int k, l;
  48         int *eq = isl_calloc_array(map->ctx, int, 2 * map->p[i]->n_eq);
  49         unsigned dim;
  50
  51         dim = isl_basic_map_total_dim(map->p[i]);
  52         for (k = 0; k < map->p[i]->n_eq; ++k) {
  53                 for (l = 0; l < 2; ++l) {
  54                         isl_seq_neg(map->p[i]->eq[k], map->p[i]->eq[k], 1+dim);
  55                         eq[2 * k + l] = status_in(map->p[i]->eq[k], tabs[j]);
  56                         if (eq[2 * k + l] == STATUS_ERROR)
  57                                 goto error;
  58                 }
  59                 if (eq[2 * k] == STATUS_SEPARATE ||
  60                     eq[2 * k + 1] == STATUS_SEPARATE)
  61                         break;
  62         }
  63
  64         return eq;
  65 error:
  66         free(eq);
  67         return NULL;
  68 }
  69
  70 /* Compute the position of the inequalities of basic map "i"
  71  * with respect to basic map "j".
  72  */
  73 static int *ineq_status_in(struct isl_map *map, int i, int j,
  74         struct isl_tab **tabs)
  75 {
  76         int k;
  77         unsigned n_eq = map->p[i]->n_eq;
  78         int *ineq = isl_calloc_array(map->ctx, int, map->p[i]->n_ineq);
  79
  80         for (k = 0; k < map->p[i]->n_ineq; ++k) {
  81                 if (isl_tab_is_redundant(tabs[i], n_eq + k)) {
  82                         ineq[k] = STATUS_REDUNDANT;
  83                         continue;
  84                 }
  85                 ineq[k] = status_in(map->p[i]->ineq[k], tabs[j]);
  86                 if (ineq[k] == STATUS_ERROR)
  87                         goto error;
  88                 if (ineq[k] == STATUS_SEPARATE)
  89                         break;
  90         }
  91
  92         return ineq;
  93 error:
  94         free(ineq);
  95         return NULL;
  96 }
  97
  98 static int any(int *con, unsigned len, int status)
  99 {
 100         int i;
 101
 102         for (i = 0; i < len ; ++i)
 103                 if (con[i] == status)
 104                         return 1;
 105         return 0;
 106 }
 107
 108 static int count(int *con, unsigned len, int status)
 109 {
 110         int i;
 111         int c = 0;
 112
 113         for (i = 0; i < len ; ++i)
 114                 if (con[i] == status)
 115                         c++;
 116         return c;
 117 }
 118
 119 static int all(int *con, unsigned len, int status)
 120 {
 121         int i;
 122
 123         for (i = 0; i < len ; ++i) {
 124                 if (con[i] == STATUS_REDUNDANT)
 125                         continue;
 126                 if (con[i] != status)
 127                         return 0;
 128         }
 129         return 1;
 130 }
 131
 132 static void drop(struct isl_map *map, int i, struct isl_tab **tabs)
 133 {
 134         isl_basic_map_free(map->p[i]);
 135         isl_tab_free(tabs[i]);
 136
 137         if (i != map->n - 1) {
 138                 map->p[i] = map->p[map->n - 1];
 139                 tabs[i] = tabs[map->n - 1];
 140         }
 141         tabs[map->n - 1] = NULL;
 142         map->n--;
 143 }
 144
 145 /* Replace the pair of basic maps i and j by the basic map bounded
 146  * by the valid constraints in both basic maps and the constraint
 147  * in extra (if not NULL).
 148  */
 149 static int fuse(struct isl_map *map, int i, int j,
 150         struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j,
 151         __isl_keep isl_mat *extra)
 152 {
 153         int k, l;
 154         struct isl_basic_map *fused = NULL;
 155         struct isl_tab *fused_tab = NULL;
 156         unsigned total = isl_basic_map_total_dim(map->p[i]);
 157         unsigned extra_rows = extra ? extra->n_row : 0;
 158
 159         fused = isl_basic_map_alloc_dim(isl_dim_copy(map->p[i]->dim),
 160                         map->p[i]->n_div,
 161                         map->p[i]->n_eq + map->p[j]->n_eq,
 162                         map->p[i]->n_ineq + map->p[j]->n_ineq + extra_rows);
 163         if (!fused)
 164                 goto error;
 165
 166         for (k = 0; k < map->p[i]->n_eq; ++k) {
 167                 if (eq_i && (eq_i[2 * k] != STATUS_VALID ||
 168                              eq_i[2 * k + 1] != STATUS_VALID))
 169                         continue;
 170                 l = isl_basic_map_alloc_equality(fused);
 171                 if (l < 0)
 172                         goto error;
 173                 isl_seq_cpy(fused->eq[l], map->p[i]->eq[k], 1 + total);
 174         }
 175
 176         for (k = 0; k < map->p[j]->n_eq; ++k) {
 177                 if (eq_j && (eq_j[2 * k] != STATUS_VALID ||
 178                              eq_j[2 * k + 1] != STATUS_VALID))
 179                         continue;
 180                 l = isl_basic_map_alloc_equality(fused);
 181                 if (l < 0)
 182                         goto error;
 183                 isl_seq_cpy(fused->eq[l], map->p[j]->eq[k], 1 + total);
 184         }
 185
 186         for (k = 0; k < map->p[i]->n_ineq; ++k) {
 187                 if (ineq_i[k] != STATUS_VALID)
 188                         continue;
 189                 l = isl_basic_map_alloc_inequality(fused);
 190                 if (l < 0)
 191                         goto error;
 192                 isl_seq_cpy(fused->ineq[l], map->p[i]->ineq[k], 1 + total);
 193         }
 194
 195         for (k = 0; k < map->p[j]->n_ineq; ++k) {
 196                 if (ineq_j[k] != STATUS_VALID)
 197                         continue;
 198                 l = isl_basic_map_alloc_inequality(fused);
 199                 if (l < 0)
 200                         goto error;
 201                 isl_seq_cpy(fused->ineq[l], map->p[j]->ineq[k], 1 + total);
 202         }
 203
 204         for (k = 0; k < map->p[i]->n_div; ++k) {
 205                 int l = isl_basic_map_alloc_div(fused);
 206                 if (l < 0)
 207                         goto error;
 208                 isl_seq_cpy(fused->div[l], map->p[i]->div[k], 1 + 1 + total);
 209         }
 210
 211         for (k = 0; k < extra_rows; ++k) {
 212                 l = isl_basic_map_alloc_inequality(fused);
 213                 if (l < 0)
 214                         goto error;
 215                 isl_seq_cpy(fused->ineq[l], extra->row[k], 1 + total);
 216         }
 217
 218         fused = isl_basic_map_gauss(fused, NULL);
 219         ISL_F_SET(fused, ISL_BASIC_MAP_FINAL);
 220         if (ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_RATIONAL) &&
 221             ISL_F_ISSET(map->p[j], ISL_BASIC_MAP_RATIONAL))
 222                 ISL_F_SET(fused, ISL_BASIC_MAP_RATIONAL);
 223
 224         fused_tab = isl_tab_from_basic_map(fused);
 225         if (isl_tab_detect_redundant(fused_tab) < 0)
 226                 goto error;
 227
 228         isl_basic_map_free(map->p[i]);
 229         map->p[i] = fused;
 230         isl_tab_free(tabs[i]);
 231         tabs[i] = fused_tab;
 232         drop(map, j, tabs);
 233
 234         return 1;
 235 error:
 236         isl_tab_free(fused_tab);
 237         isl_basic_map_free(fused);
 238         return -1;
 239 }
 240
 241 /* Given a pair of basic maps i and j such that all constraints are either
 242  * "valid" or "cut", check if the facets corresponding to the "cut"
 243  * constraints of i lie entirely within basic map j.
 244  * If so, replace the pair by the basic map consisting of the valid
 245  * constraints in both basic maps.
 246  *
 247  * To see that we are not introducing any extra points, call the
 248  * two basic maps A and B and the resulting map U and let x
 249  * be an element of U \setminus ( A \cup B ).
 250  * Then there is a pair of cut constraints c_1 and c_2 in A and B such that x
 251  * violates them.  Let X be the intersection of U with the opposites
 252  * of these constraints.  Then x \in X.
 253  * The facet corresponding to c_1 contains the corresponding facet of A.
 254  * This facet is entirely contained in B, so c_2 is valid on the facet.
 255  * However, since it is also (part of) a facet of X, -c_2 is also valid
 256  * on the facet.  This means c_2 is saturated on the facet, so c_1 and
 257  * c_2 must be opposites of each other, but then x could not violate
 258  * both of them.
 259  */
 260 static int check_facets(struct isl_map *map, int i, int j,
 261         struct isl_tab **tabs, int *ineq_i, int *ineq_j)
 262 {
 263         int k, l;
 264         struct isl_tab_undo *snap;
 265         unsigned n_eq = map->p[i]->n_eq;
 266
 267         snap = isl_tab_snap(tabs[i]);
 268
 269         for (k = 0; k < map->p[i]->n_ineq; ++k) {
 270                 if (ineq_i[k] != STATUS_CUT)
 271                         continue;
 272                 tabs[i] = isl_tab_select_facet(tabs[i], n_eq + k);
 273                 for (l = 0; l < map->p[j]->n_ineq; ++l) {
 274                         int stat;
 275                         if (ineq_j[l] != STATUS_CUT)
 276                                 continue;
 277                         stat = status_in(map->p[j]->ineq[l], tabs[i]);
 278                         if (stat != STATUS_VALID)
 279                                 break;
 280                 }
 281                 if (isl_tab_rollback(tabs[i], snap) < 0)
 282                         return -1;
 283                 if (l < map->p[j]->n_ineq)
 284                         break;
 285         }
 286
 287         if (k < map->p[i]->n_ineq)
 288                 /* BAD CUT PAIR */
 289                 return 0;
 290         return fuse(map, i, j, tabs, NULL, ineq_i, NULL, ineq_j, NULL);
 291 }
 292
 293 /* Both basic maps have at least one inequality with and adjacent
 294  * (but opposite) inequality in the other basic map.
 295  * Check that there are no cut constraints and that there is only
 296  * a single pair of adjacent inequalities.
 297  * If so, we can replace the pair by a single basic map described
 298  * by all but the pair of adjacent inequalities.
 299  * Any additional points introduced lie strictly between the two
 300  * adjacent hyperplanes and can therefore be integral.
 301  *
 302  *        ____                    _____
 303  *       /    ||\                /     \
 304  *      /     || \              /       \
 305  *      \     ||  \     =>      \        \
 306  *       \    ||  /              \       /
 307  *        \___||_/                \_____/
 308  *
 309  * The test for a single pair of adjancent inequalities is important
 310  * for avoiding the combination of two basic maps like the following
 311  *
 312  *       /|
 313  *      / |
 314  *     /__|
 315  *         _____
 316  *         |   |
 317  *         |   |
 318  *         |___|
 319  */
 320 static int check_adj_ineq(struct isl_map *map, int i, int j,
 321         struct isl_tab **tabs, int *ineq_i, int *ineq_j)
 322 {
 323         int changed = 0;
 324
 325         if (any(ineq_i, map->p[i]->n_ineq, STATUS_CUT) ||
 326             any(ineq_j, map->p[j]->n_ineq, STATUS_CUT))
 327                 /* ADJ INEQ CUT */
 328                 ;
 329         else if (count(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ) == 1 &&
 330                  count(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ) == 1)
 331                 changed = fuse(map, i, j, tabs, NULL, ineq_i, NULL, ineq_j, NULL);
 332         /* else ADJ INEQ TOO MANY */
 333
 334         return changed;
 335 }
 336
 337 /* Check if basic map "i" contains the basic map represented
 338  * by the tableau "tab".
 339  */
 340 static int contains(struct isl_map *map, int i, int *ineq_i,
 341         struct isl_tab *tab)
 342 {
 343         int k, l;
 344         unsigned dim;
 345
 346         dim = isl_basic_map_total_dim(map->p[i]);
 347         for (k = 0; k < map->p[i]->n_eq; ++k) {
 348                 for (l = 0; l < 2; ++l) {
 349                         int stat;
 350                         isl_seq_neg(map->p[i]->eq[k], map->p[i]->eq[k], 1+dim);
 351                         stat = status_in(map->p[i]->eq[k], tab);
 352                         if (stat != STATUS_VALID)
 353                                 return 0;
 354                 }
 355         }
 356
 357         for (k = 0; k < map->p[i]->n_ineq; ++k) {
 358                 int stat;
 359                 if (ineq_i[k] == STATUS_REDUNDANT)
 360                         continue;
 361                 stat = status_in(map->p[i]->ineq[k], tab);
 362                 if (stat != STATUS_VALID)
 363                         return 0;
 364         }
 365         return 1;
 366 }
 367
 368 /* Basic map "i" has an inequality "k" that is adjacent to some equality
 369  * of basic map "j".  All the other inequalities are valid for "j".
 370  * Check if basic map "j" forms an extension of basic map "i".
 371  *
 372  * In particular, we relax constraint "k", compute the corresponding
 373  * facet and check whether it is included in the other basic map.
 374  * If so, we know that relaxing the constraint extends the basic
 375  * map with exactly the other basic map (we already know that this
 376  * other basic map is included in the extension, because there
 377  * were no "cut" inequalities in "i") and we can replace the
 378  * two basic maps by thie extension.
 379  *        ____                    _____
 380  *       /    ||                 /     |
 381  *      /     ||                /      |
 382  *      \     ||        =>      \      |
 383  *       \    ||                 \     |
 384  *        \___||                  \____|
 385  */
 386 static int is_extension(struct isl_map *map, int i, int j, int k,
 387         struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)
 388 {
 389         int changed = 0;
 390         int super;
 391         struct isl_tab_undo *snap, *snap2;
 392         unsigned n_eq = map->p[i]->n_eq;
 393
 394         snap = isl_tab_snap(tabs[i]);
 395         tabs[i] = isl_tab_relax(tabs[i], n_eq + k);
 396         snap2 = isl_tab_snap(tabs[i]);
 397         tabs[i] = isl_tab_select_facet(tabs[i], n_eq + k);
 398         super = contains(map, j, ineq_j, tabs[i]);
 399         if (super) {
 400                 if (isl_tab_rollback(tabs[i], snap2) < 0)
 401                         return -1;
 402                 map->p[i] = isl_basic_map_cow(map->p[i]);
 403                 if (!map->p[i])
 404                         return -1;
 405                 isl_int_add_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);
 406                 ISL_F_SET(map->p[i], ISL_BASIC_MAP_FINAL);
 407                 drop(map, j, tabs);
 408                 changed = 1;
 409         } else
 410                 if (isl_tab_rollback(tabs[i], snap) < 0)
 411                         return -1;
 412
 413         return changed;
 414 }
 415
 416 /* For each non-redundant constraint in "bmap" (as determined by "tab"),
 417  * wrap the constraint around "bound" such that it includes the whole
 418  * set "set" and append the resulting constraint to "wraps".
 419  * "wraps" is assumed to have been pre-allocated to the appropriate size.
 420  * wraps->n_row is the number of actual wrapped constraints that have
 421  * been added.
 422  * If any of the wrapping problems results in a constraint that is
 423  * identical to "bound", then this means that "set" is unbounded in such
 424  * way that no wrapping is possible.  If this happens then wraps->n_row
 425  * is reset to zero.
 426  */
 427 static int add_wraps(__isl_keep isl_mat *wraps, __isl_keep isl_basic_map *bmap,
 428         struct isl_tab *tab, isl_int *bound, __isl_keep isl_set *set)
 429 {
 430         int l;
 431         int w;
 432         unsigned total = isl_basic_map_total_dim(bmap);
 433
 434         w = wraps->n_row;
 435
 436         for (l = 0; l < bmap->n_ineq; ++l) {
 437                 if (isl_seq_is_neg(bound, bmap->ineq[l], 1 + total))
 438                         continue;
 439                 if (isl_seq_eq(bound, bmap->ineq[l], 1 + total))
 440                         continue;
 441                 if (isl_tab_is_redundant(tab, bmap->n_eq + l))
 442                         continue;
 443
 444                 isl_seq_cpy(wraps->row[w], bound, 1 + total);
 445                 if (!isl_set_wrap_facet(set, wraps->row[w], bmap->ineq[l]))
 446                         return -1;
 447                 if (isl_seq_eq(wraps->row[w], bound, 1 + total))
 448                         goto unbounded;
 449                 ++w;
 450         }
 451         for (l = 0; l < bmap->n_eq; ++l) {
 452                 if (isl_seq_is_neg(bound, bmap->eq[l], 1 + total))
 453                         continue;
 454                 if (isl_seq_eq(bound, bmap->eq[l], 1 + total))
 455                         continue;
 456
 457                 isl_seq_cpy(wraps->row[w], bound, 1 + total);
 458                 isl_seq_neg(wraps->row[w + 1], bmap->eq[l], 1 + total);
 459                 if (!isl_set_wrap_facet(set, wraps->row[w], wraps->row[w + 1]))
 460                         return -1;
 461                 if (isl_seq_eq(wraps->row[w], bound, 1 + total))
 462                         goto unbounded;
 463                 ++w;
 464
 465                 isl_seq_cpy(wraps->row[w], bound, 1 + total);
 466                 if (!isl_set_wrap_facet(set, wraps->row[w], bmap->eq[l]))
 467                         return -1;
 468                 if (isl_seq_eq(wraps->row[w], bound, 1 + total))
 469                         goto unbounded;
 470                 ++w;
 471         }
 472
 473         wraps->n_row = w;
 474         return 0;
 475 unbounded:
 476         wraps->n_row = 0;
 477         return 0;
 478 }
 479
 480 /* Given a basic set i with a constraint k that is adjacent to either the
 481  * whole of basic set j or a facet of basic set j, check if we can wrap
 482  * both the facet corresponding to k and the facet of j (or the whole of j)
 483  * around their ridges to include the other set.
 484  * If so, replace the pair of basic sets by their union.
 485  *
 486  * All constraints of i (except k) are assumed to be valid for j.
 487  *
 488  * In the case where j has a facet adjacent to i, tab[j] is assumed
 489  * to have been restricted to this facet, so that the non-redundant
 490  * constraints in tab[j] are the ridges of the facet.
 491  * Note that for the purpose of wrapping, it does not matter whether
 492  * we wrap the ridges of i around the whole of j or just around
 493  * the facet since all the other constraints are assumed to be valid for j.
 494  * In practice, we wrap to include the whole of j.
 495  *        ____                    _____
 496  *       /    |                  /     \
 497  *      /     ||                /      |
 498  *      \     ||        =>      \      |
 499  *       \    ||                 \     |
 500  *        \___||                  \____|
 501  *
 502  */
 503 static int can_wrap_in_facet(struct isl_map *map, int i, int j, int k,
 504         struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)
 505 {
 506         int changed = 0;
 507         struct isl_mat *wraps = NULL;
 508         struct isl_set *set_i = NULL;
 509         struct isl_set *set_j = NULL;
 510         struct isl_vec *bound = NULL;
 511         unsigned total = isl_basic_map_total_dim(map->p[i]);
 512         struct isl_tab_undo *snap;
 513
 514         snap = isl_tab_snap(tabs[i]);
 515
 516         set_i = isl_set_from_basic_set(
 517                     isl_basic_map_underlying_set(isl_basic_map_copy(map->p[i])));
 518         set_j = isl_set_from_basic_set(
 519                     isl_basic_map_underlying_set(isl_basic_map_copy(map->p[j])));
 520         wraps = isl_mat_alloc(map->ctx, 2 * (map->p[i]->n_eq + map->p[j]->n_eq) +
 521                                         map->p[i]->n_ineq + map->p[j]->n_ineq,
 522                                         1 + total);
 523         bound = isl_vec_alloc(map->ctx, 1 + total);
 524         if (!set_i || !set_j || !wraps || !bound)
 525                 goto error;
 526
 527         isl_seq_cpy(bound->el, map->p[i]->ineq[k], 1 + total);
 528         isl_int_add_ui(bound->el[0], bound->el[0], 1);
 529
 530         isl_seq_cpy(wraps->row[0], bound->el, 1 + total);
 531         wraps->n_row = 1;
 532
 533         if (add_wraps(wraps, map->p[j], tabs[j], bound->el, set_i) < 0)
 534                 goto error;
 535         if (!wraps->n_row)
 536                 goto unbounded;
 537
 538         tabs[i] = isl_tab_select_facet(tabs[i], map->p[i]->n_eq + k);
 539         if (isl_tab_detect_redundant(tabs[i]) < 0)
 540                 goto error;
 541
 542         isl_seq_neg(bound->el, map->p[i]->ineq[k], 1 + total);
 543
 544         if (add_wraps(wraps, map->p[i], tabs[i], bound->el, set_j) < 0)
 545                 goto error;
 546         if (!wraps->n_row)
 547                 goto unbounded;
 548
 549         changed = fuse(map, i, j, tabs, eq_i, ineq_i, eq_j, ineq_j, wraps);
 550
 551         if (!changed) {
 552 unbounded:
 553                 if (isl_tab_rollback(tabs[i], snap) < 0)
 554                         goto error;
 555         }
 556
 557         isl_mat_free(wraps);
 558
 559         isl_set_free(set_i);
 560         isl_set_free(set_j);
 561
 562         isl_vec_free(bound);
 563
 564         return changed;
 565 error:
 566         isl_vec_free(bound);
 567         isl_mat_free(wraps);
 568         isl_set_free(set_i);
 569         isl_set_free(set_j);
 570         return -1;
 571 }
 572
 573 /* Given two basic sets i and j such that i has exactly one cut constraint,
 574  * check if we can wrap the corresponding facet around its ridges to include
 575  * the other basic set (and nothing else).
 576  * If so, replace the pair by their union.
 577  *
 578  * We first check if j has a facet adjacent to the cut constraint of i.
 579  * If so, we try to wrap in the facet.
 580  *        ____                    _____
 581  *       / ___|_                 /     \
 582  *      / |    |                /      |
 583  *      \ |    |        =>      \      |
 584  *       \|____|                 \     |
 585  *        \___|                   \____/
 586  */
 587 static int can_wrap_in_set(struct isl_map *map, int i, int j,
 588         struct isl_tab **tabs, int *ineq_i, int *ineq_j)
 589 {
 590         int changed = 0;
 591         int k, l;
 592         unsigned total = isl_basic_map_total_dim(map->p[i]);
 593         struct isl_tab_undo *snap;
 594
 595         for (k = 0; k < map->p[i]->n_ineq; ++k)
 596                 if (ineq_i[k] == STATUS_CUT)
 597                         break;
 598
 599         isl_assert(map->ctx, k < map->p[i]->n_ineq, return -1);
 600
 601         isl_int_add_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);
 602         for (l = 0; l < map->p[j]->n_ineq; ++l) {
 603                 if (isl_tab_is_redundant(tabs[j], map->p[j]->n_eq + l))
 604                         continue;
 605                 if (isl_seq_eq(map->p[i]->ineq[k],
 606                                 map->p[j]->ineq[l], 1 + total))
 607                         break;
 608         }
 609         isl_int_sub_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);
 610
 611         if (l >= map->p[j]->n_ineq)
 612                 return 0;
 613
 614         snap = isl_tab_snap(tabs[j]);
 615         tabs[j] = isl_tab_select_facet(tabs[j], map->p[j]->n_eq + l);
 616         if (isl_tab_detect_redundant(tabs[j]) < 0)
 617                 return -1;
 618
 619         changed = can_wrap_in_facet(map, i, j, k, tabs, NULL, ineq_i, NULL, ineq_j);
 620
 621         if (!changed && isl_tab_rollback(tabs[j], snap) < 0)
 622                 return -1;
 623
 624         return changed;
 625 }
 626
 627 /* Check if either i or j has a single cut constraint that can
 628  * be used to wrap in (a facet of) the other basic set.
 629  * if so, replace the pair by their union.
 630  */
 631 static int check_wrap(struct isl_map *map, int i, int j,
 632         struct isl_tab **tabs, int *ineq_i, int *ineq_j)
 633 {
 634         int changed = 0;
 635
 636         if (count(ineq_i, map->p[i]->n_ineq, STATUS_CUT) == 1)
 637                 changed = can_wrap_in_set(map, i, j, tabs, ineq_i, ineq_j);
 638         if (changed)
 639                 return changed;
 640
 641         if (count(ineq_j, map->p[j]->n_ineq, STATUS_CUT) == 1)
 642                 changed = can_wrap_in_set(map, j, i, tabs, ineq_j, ineq_i);
 643         return changed;
 644 }
 645
 646 /* At least one of the basic maps has an equality that is adjacent
 647  * to inequality.  Make sure that only one of the basic maps has
 648  * such an equality and that the other basic map has exactly one
 649  * inequality adjacent to an equality.
 650  * We call the basic map that has the inequality "i" and the basic
 651  * map that has the equality "j".
 652  * If "i" has any "cut" inequality, then relaxing the inequality
 653  * by one would not result in a basic map that contains the other
 654  * basic map.
 655  */
 656 static int check_adj_eq(struct isl_map *map, int i, int j,
 657         struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)
 658 {
 659         int changed = 0;
 660         int k;
 661
 662         if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_INEQ) &&
 663             any(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_INEQ))
 664                 /* ADJ EQ TOO MANY */
 665                 return 0;
 666
 667         if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_INEQ))
 668                 return check_adj_eq(map, j, i, tabs,
 669                                         eq_j, ineq_j, eq_i, ineq_i);
 670
 671         /* j has an equality adjacent to an inequality in i */
 672
 673         if (any(ineq_i, map->p[i]->n_ineq, STATUS_CUT))
 674                 /* ADJ EQ CUT */
 675                 return 0;
 676         if (count(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_INEQ) != 1 ||
 677             count(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_EQ) != 1 ||
 678             any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_EQ) ||
 679             any(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ) ||
 680             any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ))
 681                 /* ADJ EQ TOO MANY */
 682                 return 0;
 683
 684         for (k = 0; k < map->p[i]->n_ineq ; ++k)
 685                 if (ineq_i[k] == STATUS_ADJ_EQ)
 686                         break;
 687
 688         changed = is_extension(map, i, j, k, tabs, eq_i, ineq_i, eq_j, ineq_j);
 689         if (changed)
 690                 return changed;
 691
 692         changed = can_wrap_in_facet(map, i, j, k, tabs, eq_i, ineq_i, eq_j, ineq_j);
 693
 694         return changed;
 695 }
 696
 697 /* Check if the union of the given pair of basic maps
 698  * can be represented by a single basic map.
 699  * If so, replace the pair by the single basic map and return 1.
 700  * Otherwise, return 0;
 701  *
 702  * We first check the effect of each constraint of one basic map
 703  * on the other basic map.
 704  * The constraint may be
 705  *      redundant       the constraint is redundant in its own
 706  *                      basic map and should be ignore and removed
 707  *                      in the end
 708  *      valid           all (integer) points of the other basic map
 709  *                      satisfy the constraint
 710  *      separate        no (integer) point of the other basic map
 711  *                      satisfies the constraint
 712  *      cut             some but not all points of the other basic map
 713  *                      satisfy the constraint
 714  *      adj_eq          the given constraint is adjacent (on the outside)
 715  *                      to an equality of the other basic map
 716  *      adj_ineq        the given constraint is adjacent (on the outside)
 717  *                      to an inequality of the other basic map
 718  *
 719  * We consider six cases in which we can replace the pair by a single
 720  * basic map.  We ignore all "redundant" constraints.
 721  *
 722  *      1. all constraints of one basic map are valid
 723  *              => the other basic map is a subset and can be removed
 724  *
 725  *      2. all constraints of both basic maps are either "valid" or "cut"
 726  *         and the facets corresponding to the "cut" constraints
 727  *         of one of the basic maps lies entirely inside the other basic map
 728  *              => the pair can be replaced by a basic map consisting
 729  *                 of the valid constraints in both basic maps
 730  *
 731  *      3. there is a single pair of adjacent inequalities
 732  *         (all other constraints are "valid")
 733  *              => the pair can be replaced by a basic map consisting
 734  *                 of the valid constraints in both basic maps
 735  *
 736  *      4. there is a single adjacent pair of an inequality and an equality,
 737  *         the other constraints of the basic map containing the inequality are
 738  *         "valid".  Moreover, if the inequality the basic map is relaxed
 739  *         and then turned into an equality, then resulting facet lies
 740  *         entirely inside the other basic map
 741  *              => the pair can be replaced by the basic map containing
 742  *                 the inequality, with the inequality relaxed.
 743  *
 744  *      5. there is a single adjacent pair of an inequality and an equality,
 745  *         the other constraints of the basic map containing the inequality are
 746  *         "valid".  Moreover, the facets corresponding to both
 747  *         the inequality and the equality can be wrapped around their
 748  *         ridges to include the other basic map
 749  *              => the pair can be replaced by a basic map consisting
 750  *                 of the valid constraints in both basic maps together
 751  *                 with all wrapping constraints
 752  *
 753  *      6. one of the basic maps has a single cut constraint and
 754  *         the other basic map has a constraint adjacent to this constraint.
 755  *         Moreover, the facets corresponding to both constraints
 756  *         can be wrapped around their ridges to include the other basic map
 757  *              => the pair can be replaced by a basic map consisting
 758  *                 of the valid constraints in both basic maps together
 759  *                 with all wrapping constraints
 760  *
 761  * Throughout the computation, we maintain a collection of tableaus
 762  * corresponding to the basic maps.  When the basic maps are dropped
 763  * or combined, the tableaus are modified accordingly.
 764  */
 765 static int coalesce_pair(struct isl_map *map, int i, int j,
 766         struct isl_tab **tabs)
 767 {
 768         int changed = 0;
 769         int *eq_i = NULL;
 770         int *eq_j = NULL;
 771         int *ineq_i = NULL;
 772         int *ineq_j = NULL;
 773
 774         eq_i = eq_status_in(map, i, j, tabs);
 775         if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ERROR))
 776                 goto error;
 777         if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_SEPARATE))
 778                 goto done;
 779
 780         eq_j = eq_status_in(map, j, i, tabs);
 781         if (any(eq_j, 2 * map->p[j]->n_eq, STATUS_ERROR))
 782                 goto error;
 783         if (any(eq_j, 2 * map->p[j]->n_eq, STATUS_SEPARATE))
 784                 goto done;
 785
 786         ineq_i = ineq_status_in(map, i, j, tabs);
 787         if (any(ineq_i, map->p[i]->n_ineq, STATUS_ERROR))
 788                 goto error;
 789         if (any(ineq_i, map->p[i]->n_ineq, STATUS_SEPARATE))
 790                 goto done;
 791
 792         ineq_j = ineq_status_in(map, j, i, tabs);
 793         if (any(ineq_j, map->p[j]->n_ineq, STATUS_ERROR))
 794                 goto error;
 795         if (any(ineq_j, map->p[j]->n_ineq, STATUS_SEPARATE))
 796                 goto done;
 797
 798         if (all(eq_i, 2 * map->p[i]->n_eq, STATUS_VALID) &&
 799             all(ineq_i, map->p[i]->n_ineq, STATUS_VALID)) {
 800                 drop(map, j, tabs);
 801                 changed = 1;
 802         } else if (all(eq_j, 2 * map->p[j]->n_eq, STATUS_VALID) &&
 803                    all(ineq_j, map->p[j]->n_ineq, STATUS_VALID)) {
 804                 drop(map, i, tabs);
 805                 changed = 1;
 806         } else if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_CUT) ||
 807                    any(eq_j, 2 * map->p[j]->n_eq, STATUS_CUT)) {
 808                 /* BAD CUT */
 809         } else if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_EQ) ||
 810                    any(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_EQ)) {
 811                 /* ADJ EQ PAIR */
 812         } else if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_INEQ) ||
 813                    any(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_INEQ)) {
 814                 changed = check_adj_eq(map, i, j, tabs,
 815                                         eq_i, ineq_i, eq_j, ineq_j);
 816         } else if (any(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_EQ) ||
 817                    any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_EQ)) {
 818                 /* Can't happen */
 819                 /* BAD ADJ INEQ */
 820         } else if (any(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ) ||
 821                    any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ)) {
 822                 changed = check_adj_ineq(map, i, j, tabs, ineq_i, ineq_j);
 823         } else {
 824                 changed = check_facets(map, i, j, tabs, ineq_i, ineq_j);
 825                 if (!changed)
 826                         changed = check_wrap(map, i, j, tabs, ineq_i, ineq_j);
 827         }
 828
 829 done:
 830         free(eq_i);
 831         free(eq_j);
 832         free(ineq_i);
 833         free(ineq_j);
 834         return changed;
 835 error:
 836         free(eq_i);
 837         free(eq_j);
 838         free(ineq_i);
 839         free(ineq_j);
 840         return -1;
 841 }
 842
 843 static struct isl_map *coalesce(struct isl_map *map, struct isl_tab **tabs)
 844 {
 845         int i, j;
 846
 847         for (i = 0; i < map->n - 1; ++i)
 848                 for (j = i + 1; j < map->n; ++j) {
 849                         int changed;
 850                         changed = coalesce_pair(map, i, j, tabs);
 851                         if (changed < 0)
 852                                 goto error;
 853                         if (changed)
 854                                 return coalesce(map, tabs);
 855                 }
 856         return map;
 857 error:
 858         isl_map_free(map);
 859         return NULL;
 860 }
 861
 862 /* For each pair of basic maps in the map, check if the union of the two
 863  * can be represented by a single basic map.
 864  * If so, replace the pair by the single basic map and start over.
 865  */
 866 struct isl_map *isl_map_coalesce(struct isl_map *map)
 867 {
 868         int i;
 869         unsigned n;
 870         struct isl_tab **tabs = NULL;
 871
 872         if (!map)
 873                 return NULL;
 874
 875         if (map->n <= 1)
 876                 return map;
 877
 878         map = isl_map_align_divs(map);
 879
 880         tabs = isl_calloc_array(map->ctx, struct isl_tab *, map->n);
 881         if (!tabs)
 882                 goto error;
 883
 884         n = map->n;
 885         for (i = 0; i < map->n; ++i) {
 886                 tabs[i] = isl_tab_from_basic_map(map->p[i]);
 887                 if (!tabs[i])
 888                         goto error;
 889                 if (!ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_NO_IMPLICIT))
 890                         tabs[i] = isl_tab_detect_implicit_equalities(tabs[i]);
 891                 if (!ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_NO_REDUNDANT))
 892                         if (isl_tab_detect_redundant(tabs[i]) < 0)
 893                                 goto error;
 894         }
 895         for (i = map->n - 1; i >= 0; --i)
 896                 if (tabs[i]->empty)
 897                         drop(map, i, tabs);
 898
 899         map = coalesce(map, tabs);
 900
 901         if (map)
 902                 for (i = 0; i < map->n; ++i) {
 903                         map->p[i] = isl_basic_map_update_from_tab(map->p[i],
 904                                                                     tabs[i]);
 905                         map->p[i] = isl_basic_map_finalize(map->p[i]);
 906                         if (!map->p[i])
 907                                 goto error;
 908                         ISL_F_SET(map->p[i], ISL_BASIC_MAP_NO_IMPLICIT);
 909                         ISL_F_SET(map->p[i], ISL_BASIC_MAP_NO_REDUNDANT);
 910                 }
 911
 912         for (i = 0; i < n; ++i)
 913                 isl_tab_free(tabs[i]);
 914
 915         free(tabs);
 916
 917         return map;
 918 error:
 919         if (tabs)
 920                 for (i = 0; i < n; ++i)
 921                         isl_tab_free(tabs[i]);
 922         free(tabs);
 923         return NULL;
 924 }
 925
 926 /* For each pair of basic sets in the set, check if the union of the two
 927  * can be represented by a single basic set.
 928  * If so, replace the pair by the single basic set and start over.
 929  */
 930 struct isl_set *isl_set_coalesce(struct isl_set *set)
 931 {
 932         return (struct isl_set *)isl_map_coalesce((struct isl_map *)set);
 933 }