add isl_basic_set_positive_orthant
[isl.git] / isl_coalesce.c
blob1ccd208b2e58da36845b9f1a08b8b288e7ea052e
1 #include "isl_map_private.h"
2 #include "isl_tab.h"
4 #define STATUS_ERROR -1
5 #define STATUS_REDUNDANT 1
6 #define STATUS_VALID 2
7 #define STATUS_SEPARATE 3
8 #define STATUS_CUT 4
9 #define STATUS_ADJ_EQ 5
10 #define STATUS_ADJ_INEQ 6
12 static int status_in(isl_int *ineq, struct isl_tab *tab)
14 enum isl_ineq_type type = isl_tab_ineq_type(tab, ineq);
15 switch (type) {
16 case isl_ineq_error: return STATUS_ERROR;
17 case isl_ineq_redundant: return STATUS_VALID;
18 case isl_ineq_separate: return STATUS_SEPARATE;
19 case isl_ineq_cut: return STATUS_CUT;
20 case isl_ineq_adj_eq: return STATUS_ADJ_EQ;
21 case isl_ineq_adj_ineq: return STATUS_ADJ_INEQ;
25 /* Compute the position of the equalities of basic map "i"
26 * with respect to basic map "j".
27 * The resulting array has twice as many entries as the number
28 * of equalities corresponding to the two inequalties to which
29 * each equality corresponds.
31 static int *eq_status_in(struct isl_map *map, int i, int j,
32 struct isl_tab **tabs)
34 int k, l;
35 int *eq = isl_calloc_array(map->ctx, int, 2 * map->p[i]->n_eq);
36 unsigned dim;
38 dim = isl_basic_map_total_dim(map->p[i]);
39 for (k = 0; k < map->p[i]->n_eq; ++k) {
40 for (l = 0; l < 2; ++l) {
41 isl_seq_neg(map->p[i]->eq[k], map->p[i]->eq[k], 1+dim);
42 eq[2 * k + l] = status_in(map->p[i]->eq[k], tabs[j]);
43 if (eq[2 * k + l] == STATUS_ERROR)
44 goto error;
46 if (eq[2 * k] == STATUS_SEPARATE ||
47 eq[2 * k + 1] == STATUS_SEPARATE)
48 break;
51 return eq;
52 error:
53 free(eq);
54 return NULL;
57 /* Compute the position of the inequalities of basic map "i"
58 * with respect to basic map "j".
60 static int *ineq_status_in(struct isl_map *map, int i, int j,
61 struct isl_tab **tabs)
63 int k;
64 unsigned n_eq = map->p[i]->n_eq;
65 int *ineq = isl_calloc_array(map->ctx, int, map->p[i]->n_ineq);
67 for (k = 0; k < map->p[i]->n_ineq; ++k) {
68 if (isl_tab_is_redundant(tabs[i], n_eq + k)) {
69 ineq[k] = STATUS_REDUNDANT;
70 continue;
72 ineq[k] = status_in(map->p[i]->ineq[k], tabs[j]);
73 if (ineq[k] == STATUS_ERROR)
74 goto error;
75 if (ineq[k] == STATUS_SEPARATE)
76 break;
79 return ineq;
80 error:
81 free(ineq);
82 return NULL;
85 static int any(int *con, unsigned len, int status)
87 int i;
89 for (i = 0; i < len ; ++i)
90 if (con[i] == status)
91 return 1;
92 return 0;
95 static int count(int *con, unsigned len, int status)
97 int i;
98 int c = 0;
100 for (i = 0; i < len ; ++i)
101 if (con[i] == status)
102 c++;
103 return c;
106 static int all(int *con, unsigned len, int status)
108 int i;
110 for (i = 0; i < len ; ++i) {
111 if (con[i] == STATUS_REDUNDANT)
112 continue;
113 if (con[i] != status)
114 return 0;
116 return 1;
119 static void drop(struct isl_map *map, int i, struct isl_tab **tabs)
121 isl_basic_map_free(map->p[i]);
122 isl_tab_free(tabs[i]);
124 if (i != map->n - 1) {
125 map->p[i] = map->p[map->n - 1];
126 tabs[i] = tabs[map->n - 1];
128 tabs[map->n - 1] = NULL;
129 map->n--;
132 /* Replace the pair of basic maps i and j but the basic map bounded
133 * by the valid constraints in both basic maps.
135 static int fuse(struct isl_map *map, int i, int j, struct isl_tab **tabs,
136 int *ineq_i, int *ineq_j)
138 int k, l;
139 struct isl_basic_map *fused = NULL;
140 struct isl_tab *fused_tab = NULL;
141 unsigned total = isl_basic_map_total_dim(map->p[i]);
143 fused = isl_basic_map_alloc_dim(isl_dim_copy(map->p[i]->dim),
144 map->p[i]->n_div,
145 map->p[i]->n_eq + map->p[j]->n_eq,
146 map->p[i]->n_ineq + map->p[j]->n_ineq);
147 if (!fused)
148 goto error;
150 for (k = 0; k < map->p[i]->n_eq; ++k) {
151 int l = isl_basic_map_alloc_equality(fused);
152 isl_seq_cpy(fused->eq[l], map->p[i]->eq[k], 1 + total);
155 for (k = 0; k < map->p[j]->n_eq; ++k) {
156 int l = isl_basic_map_alloc_equality(fused);
157 isl_seq_cpy(fused->eq[l], map->p[j]->eq[k], 1 + total);
160 for (k = 0; k < map->p[i]->n_ineq; ++k) {
161 if (ineq_i[k] != STATUS_VALID)
162 continue;
163 l = isl_basic_map_alloc_inequality(fused);
164 isl_seq_cpy(fused->ineq[l], map->p[i]->ineq[k], 1 + total);
167 for (k = 0; k < map->p[j]->n_ineq; ++k) {
168 if (ineq_j[k] != STATUS_VALID)
169 continue;
170 l = isl_basic_map_alloc_inequality(fused);
171 isl_seq_cpy(fused->ineq[l], map->p[j]->ineq[k], 1 + total);
174 for (k = 0; k < map->p[i]->n_div; ++k) {
175 int l = isl_basic_map_alloc_div(fused);
176 isl_seq_cpy(fused->div[l], map->p[i]->div[k], 1 + 1 + total);
179 fused = isl_basic_map_gauss(fused, NULL);
180 ISL_F_SET(fused, ISL_BASIC_MAP_FINAL);
181 if (ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_RATIONAL) &&
182 ISL_F_ISSET(map->p[j], ISL_BASIC_MAP_RATIONAL))
183 ISL_F_SET(fused, ISL_BASIC_MAP_RATIONAL);
185 fused_tab = isl_tab_from_basic_map(fused);
186 fused_tab = isl_tab_detect_redundant(fused_tab);
187 if (!fused_tab)
188 goto error;
190 isl_basic_map_free(map->p[i]);
191 map->p[i] = fused;
192 isl_tab_free(tabs[i]);
193 tabs[i] = fused_tab;
194 drop(map, j, tabs);
196 return 1;
197 error:
198 isl_basic_map_free(fused);
199 return -1;
202 /* Given a pair of basic maps i and j such that all constraints are either
203 * "valid" or "cut", check if the facets corresponding to the "cut"
204 * constraints of i lie entirely within basic map j.
205 * If so, replace the pair by the basic map consisting of the valid
206 * constraints in both basic maps.
208 * To see that we are not introducing any extra points, call the
209 * two basic maps A and B and the resulting map U and let x
210 * be an element of U \setminus ( A \cup B ).
211 * Then there is a pair of cut constraints c_1 and c_2 in A and B such that x
212 * violates them. Let X be the intersection of U with the opposites
213 * of these constraints. Then x \in X.
214 * The facet corresponding to c_1 contains the corresponding facet of A.
215 * This facet is entirely contained in B, so c_2 is valid on the facet.
216 * However, since it is also (part of) a facet of X, -c_2 is also valid
217 * on the facet. This means c_2 is saturated on the facet, so c_1 and
218 * c_2 must be opposites of each other, but then x could not violate
219 * both of them.
221 static int check_facets(struct isl_map *map, int i, int j,
222 struct isl_tab **tabs, int *ineq_i, int *ineq_j)
224 int k, l;
225 struct isl_tab_undo *snap;
226 unsigned n_eq = map->p[i]->n_eq;
228 snap = isl_tab_snap(tabs[i]);
230 for (k = 0; k < map->p[i]->n_ineq; ++k) {
231 if (ineq_i[k] != STATUS_CUT)
232 continue;
233 tabs[i] = isl_tab_select_facet(tabs[i], n_eq + k);
234 for (l = 0; l < map->p[j]->n_ineq; ++l) {
235 int stat;
236 if (ineq_j[l] != STATUS_CUT)
237 continue;
238 stat = status_in(map->p[j]->ineq[l], tabs[i]);
239 if (stat != STATUS_VALID)
240 break;
242 isl_tab_rollback(tabs[i], snap);
243 if (l < map->p[j]->n_ineq)
244 break;
247 if (k < map->p[i]->n_ineq)
248 /* BAD CUT PAIR */
249 return 0;
250 return fuse(map, i, j, tabs, ineq_i, ineq_j);
253 /* Both basic maps have at least one inequality with and adjacent
254 * (but opposite) inequality in the other basic map.
255 * Check that there are no cut constraints and that there is only
256 * a single pair of adjacent inequalities.
257 * If so, we can replace the pair by a single basic map described
258 * by all but the pair of adjacent inequalities.
259 * Any additional points introduced lie strictly between the two
260 * adjacent hyperplanes and can therefore be integral.
262 * ____ _____
263 * / ||\ / \
264 * / || \ / \
265 * \ || \ => \ \
266 * \ || / \ /
267 * \___||_/ \_____/
269 * The test for a single pair of adjancent inequalities is important
270 * for avoiding the combination of two basic maps like the following
272 * /|
273 * / |
274 * /__|
275 * _____
276 * | |
277 * | |
278 * |___|
280 static int check_adj_ineq(struct isl_map *map, int i, int j,
281 struct isl_tab **tabs, int *ineq_i, int *ineq_j)
283 int changed = 0;
285 if (any(ineq_i, map->p[i]->n_ineq, STATUS_CUT) ||
286 any(ineq_j, map->p[j]->n_ineq, STATUS_CUT))
287 /* ADJ INEQ CUT */
289 else if (count(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ) == 1 &&
290 count(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ) == 1)
291 changed = fuse(map, i, j, tabs, ineq_i, ineq_j);
292 /* else ADJ INEQ TOO MANY */
294 return changed;
297 /* Check if basic map "i" contains the basic map represented
298 * by the tableau "tab".
300 static int contains(struct isl_map *map, int i, int *ineq_i,
301 struct isl_tab *tab)
303 int k, l;
304 unsigned dim;
306 dim = isl_basic_map_total_dim(map->p[i]);
307 for (k = 0; k < map->p[i]->n_eq; ++k) {
308 for (l = 0; l < 2; ++l) {
309 int stat;
310 isl_seq_neg(map->p[i]->eq[k], map->p[i]->eq[k], 1+dim);
311 stat = status_in(map->p[i]->eq[k], tab);
312 if (stat != STATUS_VALID)
313 return 0;
317 for (k = 0; k < map->p[i]->n_ineq; ++k) {
318 int stat;
319 if (ineq_i[k] == STATUS_REDUNDANT)
320 continue;
321 stat = status_in(map->p[i]->ineq[k], tab);
322 if (stat != STATUS_VALID)
323 return 0;
325 return 1;
328 /* At least one of the basic maps has an equality that is adjacent
329 * to inequality. Make sure that only one of the basic maps has
330 * such an equality and that the other basic map has exactly one
331 * inequality adjacent to an equality.
332 * We call the basic map that has the inequality "i" and the basic
333 * map that has the equality "j".
334 * If "i" has any "cut" inequality, then relaxing the inequality
335 * by one would not result in a basic map that contains the other
336 * basic map.
337 * Otherwise, we relax the constraint, compute the corresponding
338 * facet and check whether it is included in the other basic map.
339 * If so, we know that relaxing the constraint extend the basic
340 * map with exactly the other basic map (we already know that this
341 * other basic map is included in the extension, because there
342 * were no "cut" inequalities in "i") and we can replace the
343 * two basic maps by thie extension.
344 * ____ _____
345 * / || / |
346 * / || / |
347 * \ || => \ |
348 * \ || \ |
349 * \___|| \____|
351 static int check_adj_eq(struct isl_map *map, int i, int j,
352 struct isl_tab **tabs, int *eq_i, int *ineq_i, int *eq_j, int *ineq_j)
354 int changed = 0;
355 int super;
356 int k;
357 struct isl_tab_undo *snap, *snap2;
358 unsigned n_eq = map->p[i]->n_eq;
360 if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_INEQ) &&
361 any(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_INEQ))
362 /* ADJ EQ TOO MANY */
363 return 0;
365 if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_INEQ))
366 return check_adj_eq(map, j, i, tabs,
367 eq_j, ineq_j, eq_i, ineq_i);
369 /* j has an equality adjacent to an inequality in i */
371 if (any(ineq_i, map->p[i]->n_ineq, STATUS_CUT))
372 /* ADJ EQ CUT */
373 return 0;
374 if (count(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_INEQ) != 1 ||
375 count(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_EQ) != 1 ||
376 any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_EQ) ||
377 any(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ) ||
378 any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ))
379 /* ADJ EQ TOO MANY */
380 return 0;
382 for (k = 0; k < map->p[i]->n_ineq ; ++k)
383 if (ineq_i[k] == STATUS_ADJ_EQ)
384 break;
386 snap = isl_tab_snap(tabs[i]);
387 tabs[i] = isl_tab_relax(tabs[i], n_eq + k);
388 snap2 = isl_tab_snap(tabs[i]);
389 tabs[i] = isl_tab_select_facet(tabs[i], n_eq + k);
390 super = contains(map, j, ineq_j, tabs[i]);
391 if (super) {
392 isl_tab_rollback(tabs[i], snap2);
393 map->p[i] = isl_basic_map_cow(map->p[i]);
394 if (!map->p[i])
395 return -1;
396 isl_int_add_ui(map->p[i]->ineq[k][0], map->p[i]->ineq[k][0], 1);
397 ISL_F_SET(map->p[i], ISL_BASIC_MAP_FINAL);
398 drop(map, j, tabs);
399 changed = 1;
400 } else
401 isl_tab_rollback(tabs[i], snap);
403 return changed;
406 /* Check if the union of the given pair of basic maps
407 * can be represented by a single basic map.
408 * If so, replace the pair by the single basic map and return 1.
409 * Otherwise, return 0;
411 * We first check the effect of each constraint of one basic map
412 * on the other basic map.
413 * The constraint may be
414 * redundant the constraint is redundant in its own
415 * basic map and should be ignore and removed
416 * in the end
417 * valid all (integer) points of the other basic map
418 * satisfy the constraint
419 * separate no (integer) point of the other basic map
420 * satisfies the constraint
421 * cut some but not all points of the other basic map
422 * satisfy the constraint
423 * adj_eq the given constraint is adjacent (on the outside)
424 * to an equality of the other basic map
425 * adj_ineq the given constraint is adjacent (on the outside)
426 * to an inequality of the other basic map
428 * We consider four cases in which we can replace the pair by a single
429 * basic map. We ignore all "redundant" constraints.
431 * 1. all constraints of one basic map are valid
432 * => the other basic map is a subset and can be removed
434 * 2. all constraints of both basic maps are either "valid" or "cut"
435 * and the facets corresponding to the "cut" constraints
436 * of one of the basic maps lies entirely inside the other basic map
437 * => the pair can be replaced by a basic map consisting
438 * of the valid constraints in both basic maps
440 * 3. there is a single pair of adjacent inequalities
441 * (all other constraints are "valid")
442 * => the pair can be replaced by a basic map consisting
443 * of the valid constraints in both basic maps
445 * 4. there is a single adjacent pair of an inequality and an equality,
446 * the other constraints of the basic map containing the inequality are
447 * "valid". Moreover, if the inequality the basic map is relaxed
448 * and then turned into an equality, then resulting facet lies
449 * entirely inside the other basic map
450 * => the pair can be replaced by the basic map containing
451 * the inequality, with the inequality relaxed.
453 * Throughout the computation, we maintain a collection of tableaus
454 * corresponding to the basic maps. When the basic maps are dropped
455 * or combined, the tableaus are modified accordingly.
457 static int coalesce_pair(struct isl_map *map, int i, int j,
458 struct isl_tab **tabs)
460 int changed = 0;
461 int *eq_i = NULL;
462 int *eq_j = NULL;
463 int *ineq_i = NULL;
464 int *ineq_j = NULL;
466 eq_i = eq_status_in(map, i, j, tabs);
467 if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ERROR))
468 goto error;
469 if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_SEPARATE))
470 goto done;
472 eq_j = eq_status_in(map, j, i, tabs);
473 if (any(eq_j, 2 * map->p[j]->n_eq, STATUS_ERROR))
474 goto error;
475 if (any(eq_j, 2 * map->p[j]->n_eq, STATUS_SEPARATE))
476 goto done;
478 ineq_i = ineq_status_in(map, i, j, tabs);
479 if (any(ineq_i, map->p[i]->n_ineq, STATUS_ERROR))
480 goto error;
481 if (any(ineq_i, map->p[i]->n_ineq, STATUS_SEPARATE))
482 goto done;
484 ineq_j = ineq_status_in(map, j, i, tabs);
485 if (any(ineq_j, map->p[j]->n_ineq, STATUS_ERROR))
486 goto error;
487 if (any(ineq_j, map->p[j]->n_ineq, STATUS_SEPARATE))
488 goto done;
490 if (all(eq_i, 2 * map->p[i]->n_eq, STATUS_VALID) &&
491 all(ineq_i, map->p[i]->n_ineq, STATUS_VALID)) {
492 drop(map, j, tabs);
493 changed = 1;
494 } else if (all(eq_j, 2 * map->p[j]->n_eq, STATUS_VALID) &&
495 all(ineq_j, map->p[j]->n_ineq, STATUS_VALID)) {
496 drop(map, i, tabs);
497 changed = 1;
498 } else if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_CUT) ||
499 any(eq_j, 2 * map->p[j]->n_eq, STATUS_CUT)) {
500 /* BAD CUT */
501 } else if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_EQ) ||
502 any(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_EQ)) {
503 /* ADJ EQ PAIR */
504 } else if (any(eq_i, 2 * map->p[i]->n_eq, STATUS_ADJ_INEQ) ||
505 any(eq_j, 2 * map->p[j]->n_eq, STATUS_ADJ_INEQ)) {
506 changed = check_adj_eq(map, i, j, tabs,
507 eq_i, ineq_i, eq_j, ineq_j);
508 } else if (any(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_EQ) ||
509 any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_EQ)) {
510 /* Can't happen */
511 /* BAD ADJ INEQ */
512 } else if (any(ineq_i, map->p[i]->n_ineq, STATUS_ADJ_INEQ) ||
513 any(ineq_j, map->p[j]->n_ineq, STATUS_ADJ_INEQ)) {
514 changed = check_adj_ineq(map, i, j, tabs, ineq_i, ineq_j);
515 } else
516 changed = check_facets(map, i, j, tabs, ineq_i, ineq_j);
518 done:
519 free(eq_i);
520 free(eq_j);
521 free(ineq_i);
522 free(ineq_j);
523 return changed;
524 error:
525 free(eq_i);
526 free(eq_j);
527 free(ineq_i);
528 free(ineq_j);
529 return -1;
532 static struct isl_map *coalesce(struct isl_map *map, struct isl_tab **tabs)
534 int i, j;
536 for (i = 0; i < map->n - 1; ++i)
537 for (j = i + 1; j < map->n; ++j) {
538 int changed;
539 changed = coalesce_pair(map, i, j, tabs);
540 if (changed < 0)
541 goto error;
542 if (changed)
543 return coalesce(map, tabs);
545 return map;
546 error:
547 isl_map_free(map);
548 return NULL;
551 /* For each pair of basic maps in the map, check if the union of the two
552 * can be represented by a single basic map.
553 * If so, replace the pair by the single basic map and start over.
555 struct isl_map *isl_map_coalesce(struct isl_map *map)
557 int i;
558 unsigned n;
559 struct isl_tab **tabs = NULL;
561 if (!map)
562 return NULL;
564 if (map->n <= 1)
565 return map;
567 map = isl_map_align_divs(map);
569 tabs = isl_calloc_array(map->ctx, struct isl_tab *, map->n);
570 if (!tabs)
571 goto error;
573 n = map->n;
574 for (i = 0; i < map->n; ++i) {
575 tabs[i] = isl_tab_from_basic_map(map->p[i]);
576 if (!tabs[i])
577 goto error;
578 if (!ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_NO_IMPLICIT))
579 tabs[i] = isl_tab_detect_equalities(tabs[i]);
580 if (!ISL_F_ISSET(map->p[i], ISL_BASIC_MAP_NO_REDUNDANT))
581 tabs[i] = isl_tab_detect_redundant(tabs[i]);
583 for (i = map->n - 1; i >= 0; --i)
584 if (tabs[i]->empty)
585 drop(map, i, tabs);
587 map = coalesce(map, tabs);
589 if (map)
590 for (i = 0; i < map->n; ++i) {
591 map->p[i] = isl_basic_map_update_from_tab(map->p[i],
592 tabs[i]);
593 map->p[i] = isl_basic_map_finalize(map->p[i]);
594 if (!map->p[i])
595 goto error;
596 ISL_F_SET(map->p[i], ISL_BASIC_MAP_NO_IMPLICIT);
597 ISL_F_SET(map->p[i], ISL_BASIC_MAP_NO_REDUNDANT);
600 for (i = 0; i < n; ++i)
601 isl_tab_free(tabs[i]);
603 free(tabs);
605 return map;
606 error:
607 if (tabs)
608 for (i = 0; i < n; ++i)
609 isl_tab_free(tabs[i]);
610 free(tabs);
611 return NULL;
614 /* For each pair of basic sets in the set, check if the union of the two
615 * can be represented by a single basic set.
616 * If so, replace the pair by the single basic set and start over.
618 struct isl_set *isl_set_coalesce(struct isl_set *set)
620 (struct isl_set *)isl_map_coalesce((struct isl_map *)set);