isl_mat_parameter_compression: add memory management annotations
[isl.git] / isl_vertices.c
blobf70e410daac03d817ce5078d58624caa4d7a0c12
1 /*
2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the MIT license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8 * 91893 Orsay, France
9 */
11 #include <isl_map_private.h>
12 #include <isl_aff_private.h>
13 #include <isl/set.h>
14 #include <isl_seq.h>
15 #include <isl_tab.h>
16 #include <isl_space_private.h>
17 #include <isl_morph.h>
18 #include <isl_vertices_private.h>
19 #include <isl_mat_private.h>
20 #include <isl_vec_private.h>
22 #define SELECTED 1
23 #define DESELECTED -1
24 #define UNSELECTED 0
26 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
27 __isl_take isl_vertices *vertices);
29 __isl_give isl_vertices *isl_vertices_copy(__isl_keep isl_vertices *vertices)
31 if (!vertices)
32 return NULL;
34 vertices->ref++;
35 return vertices;
38 __isl_null isl_vertices *isl_vertices_free(__isl_take isl_vertices *vertices)
40 int i;
42 if (!vertices)
43 return NULL;
45 if (--vertices->ref > 0)
46 return NULL;
48 for (i = 0; i < vertices->n_vertices; ++i) {
49 isl_basic_set_free(vertices->v[i].vertex);
50 isl_basic_set_free(vertices->v[i].dom);
52 free(vertices->v);
54 for (i = 0; i < vertices->n_chambers; ++i) {
55 free(vertices->c[i].vertices);
56 isl_basic_set_free(vertices->c[i].dom);
58 free(vertices->c);
60 isl_basic_set_free(vertices->bset);
61 free(vertices);
63 return NULL;
66 struct isl_vertex_list {
67 struct isl_vertex v;
68 struct isl_vertex_list *next;
71 static struct isl_vertex_list *free_vertex_list(struct isl_vertex_list *list)
73 struct isl_vertex_list *next;
75 for (; list; list = next) {
76 next = list->next;
77 isl_basic_set_free(list->v.vertex);
78 isl_basic_set_free(list->v.dom);
79 free(list);
82 return NULL;
85 static __isl_give isl_vertices *vertices_from_list(__isl_keep isl_basic_set *bset,
86 int n_vertices, struct isl_vertex_list *list)
88 int i;
89 struct isl_vertex_list *next;
90 isl_vertices *vertices;
92 vertices = isl_calloc_type(bset->ctx, isl_vertices);
93 if (!vertices)
94 goto error;
95 vertices->ref = 1;
96 vertices->bset = isl_basic_set_copy(bset);
97 vertices->v = isl_alloc_array(bset->ctx, struct isl_vertex, n_vertices);
98 if (n_vertices && !vertices->v)
99 goto error;
100 vertices->n_vertices = n_vertices;
102 for (i = 0; list; list = next, i++) {
103 next = list->next;
104 vertices->v[i] = list->v;
105 free(list);
108 return vertices;
109 error:
110 isl_vertices_free(vertices);
111 free_vertex_list(list);
112 return NULL;
115 /* Prepend a vertex to the linked list "list" based on the equalities in "tab".
116 * Return isl_bool_true if the vertex was actually added and
117 * isl_bool_false otherwise.
118 * In particular, vertices with a lower-dimensional activity domain are
119 * not added to the list because they would not be included in any chamber.
120 * Return isl_bool_error on error.
122 static isl_bool add_vertex(struct isl_vertex_list **list,
123 __isl_keep isl_basic_set *bset, struct isl_tab *tab)
125 unsigned nvar;
126 struct isl_vertex_list *v = NULL;
128 if (isl_tab_detect_implicit_equalities(tab) < 0)
129 return isl_bool_error;
131 nvar = isl_basic_set_dim(bset, isl_dim_set);
133 v = isl_calloc_type(tab->mat->ctx, struct isl_vertex_list);
134 if (!v)
135 goto error;
137 v->v.vertex = isl_basic_set_copy(bset);
138 v->v.vertex = isl_basic_set_cow(v->v.vertex);
139 v->v.vertex = isl_basic_set_update_from_tab(v->v.vertex, tab);
140 v->v.vertex = isl_basic_set_simplify(v->v.vertex);
141 v->v.vertex = isl_basic_set_finalize(v->v.vertex);
142 if (!v->v.vertex)
143 goto error;
144 isl_assert(bset->ctx, v->v.vertex->n_eq >= nvar, goto error);
145 v->v.dom = isl_basic_set_copy(v->v.vertex);
146 v->v.dom = isl_basic_set_params(v->v.dom);
147 if (!v->v.dom)
148 goto error;
150 if (v->v.dom->n_eq > 0) {
151 free_vertex_list(v);
152 return isl_bool_false;
155 v->next = *list;
156 *list = v;
158 return isl_bool_true;
159 error:
160 free_vertex_list(v);
161 return isl_bool_error;
164 /* Compute the parametric vertices and the chamber decomposition
165 * of an empty parametric polytope.
167 static __isl_give isl_vertices *vertices_empty(__isl_keep isl_basic_set *bset)
169 isl_vertices *vertices;
171 if (!bset)
172 return NULL;
174 vertices = isl_calloc_type(bset->ctx, isl_vertices);
175 if (!vertices)
176 return NULL;
177 vertices->bset = isl_basic_set_copy(bset);
178 vertices->ref = 1;
180 vertices->n_vertices = 0;
181 vertices->n_chambers = 0;
183 return vertices;
186 /* Compute the parametric vertices and the chamber decomposition
187 * of the parametric polytope defined using the same constraints
188 * as "bset" in the 0D case.
189 * There is exactly one 0D vertex and a single chamber containing
190 * the vertex.
192 static __isl_give isl_vertices *vertices_0D(__isl_keep isl_basic_set *bset)
194 isl_vertices *vertices;
196 if (!bset)
197 return NULL;
199 vertices = isl_calloc_type(bset->ctx, isl_vertices);
200 if (!vertices)
201 return NULL;
202 vertices->ref = 1;
203 vertices->bset = isl_basic_set_copy(bset);
205 vertices->v = isl_calloc_array(bset->ctx, struct isl_vertex, 1);
206 if (!vertices->v)
207 goto error;
208 vertices->n_vertices = 1;
209 vertices->v[0].vertex = isl_basic_set_copy(bset);
210 vertices->v[0].dom = isl_basic_set_params(isl_basic_set_copy(bset));
211 if (!vertices->v[0].vertex || !vertices->v[0].dom)
212 goto error;
214 vertices->c = isl_calloc_array(bset->ctx, struct isl_chamber, 1);
215 if (!vertices->c)
216 goto error;
217 vertices->n_chambers = 1;
218 vertices->c[0].n_vertices = 1;
219 vertices->c[0].vertices = isl_calloc_array(bset->ctx, int, 1);
220 if (!vertices->c[0].vertices)
221 goto error;
222 vertices->c[0].dom = isl_basic_set_copy(vertices->v[0].dom);
223 if (!vertices->c[0].dom)
224 goto error;
226 return vertices;
227 error:
228 isl_vertices_free(vertices);
229 return NULL;
232 static int isl_mat_rank(__isl_keep isl_mat *mat)
234 int row, col;
235 isl_mat *H;
237 H = isl_mat_left_hermite(isl_mat_copy(mat), 0, NULL, NULL);
238 if (!H)
239 return -1;
241 for (col = 0; col < H->n_col; ++col) {
242 for (row = 0; row < H->n_row; ++row)
243 if (!isl_int_is_zero(H->row[row][col]))
244 break;
245 if (row == H->n_row)
246 break;
249 isl_mat_free(H);
251 return col;
254 /* Is the row pointed to by "f" linearly independent of the "n" first
255 * rows in "facets"?
257 static int is_independent(__isl_keep isl_mat *facets, int n, isl_int *f)
259 int rank;
261 if (isl_seq_first_non_zero(f, facets->n_col) < 0)
262 return 0;
264 isl_seq_cpy(facets->row[n], f, facets->n_col);
265 facets->n_row = n + 1;
266 rank = isl_mat_rank(facets);
267 if (rank < 0)
268 return -1;
270 return rank == n + 1;
273 /* Check whether we can select constraint "level", given the current selection
274 * reflected by facets in "tab", the rows of "facets" and the earlier
275 * "selected" elements of "selection".
277 * If the constraint is (strictly) redundant in the tableau, selecting it would
278 * result in an empty tableau, so it can't be selected.
279 * If the set variable part of the constraint is not linearly independent
280 * of the set variable parts of the already selected constraints,
281 * the constraint cannot be selected.
282 * If selecting the constraint results in an empty tableau, the constraint
283 * cannot be selected.
284 * Finally, if selecting the constraint results in some explicitly
285 * deselected constraints turning into equalities, then the corresponding
286 * vertices have already been generated, so the constraint cannot be selected.
288 static int can_select(__isl_keep isl_basic_set *bset, int level,
289 struct isl_tab *tab, __isl_keep isl_mat *facets, int selected,
290 int *selection)
292 int i;
293 int indep;
294 unsigned ovar;
295 struct isl_tab_undo *snap;
297 if (isl_tab_is_redundant(tab, level))
298 return 0;
300 ovar = isl_space_offset(bset->dim, isl_dim_set);
302 indep = is_independent(facets, selected, bset->ineq[level] + 1 + ovar);
303 if (indep < 0)
304 return -1;
305 if (!indep)
306 return 0;
308 snap = isl_tab_snap(tab);
309 if (isl_tab_select_facet(tab, level) < 0)
310 return -1;
312 if (tab->empty) {
313 if (isl_tab_rollback(tab, snap) < 0)
314 return -1;
315 return 0;
318 for (i = 0; i < level; ++i) {
319 int sgn;
321 if (selection[i] != DESELECTED)
322 continue;
324 if (isl_tab_is_equality(tab, i))
325 sgn = 0;
326 else if (isl_tab_is_redundant(tab, i))
327 sgn = 1;
328 else
329 sgn = isl_tab_sign_of_max(tab, i);
330 if (sgn < -1)
331 return -1;
332 if (sgn <= 0) {
333 if (isl_tab_rollback(tab, snap) < 0)
334 return -1;
335 return 0;
339 return 1;
342 /* Compute the parametric vertices and the chamber decomposition
343 * of a parametric polytope that is not full-dimensional.
345 * Simply map the parametric polytope to a lower dimensional space
346 * and map the resulting vertices back.
348 static __isl_give isl_vertices *lower_dim_vertices(
349 __isl_keep isl_basic_set *bset)
351 isl_morph *morph;
352 isl_vertices *vertices;
354 bset = isl_basic_set_copy(bset);
355 morph = isl_basic_set_full_compression(bset);
356 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
358 vertices = isl_basic_set_compute_vertices(bset);
359 isl_basic_set_free(bset);
361 morph = isl_morph_inverse(morph);
363 vertices = isl_morph_vertices(morph, vertices);
365 return vertices;
368 /* Compute the parametric vertices and the chamber decomposition
369 * of the parametric polytope defined using the same constraints
370 * as "bset". "bset" is assumed to have no existentially quantified
371 * variables.
373 * The vertices themselves are computed in a fairly simplistic way.
374 * We simply run through all combinations of d constraints,
375 * with d the number of set variables, and check if those d constraints
376 * define a vertex. To avoid the generation of duplicate vertices,
377 * which we may happen if a vertex is defined by more that d constraints,
378 * we make sure we only generate the vertex for the d constraints with
379 * smallest index.
381 * We set up a tableau and keep track of which facets have been
382 * selected. The tableau is marked strict_redundant so that we can be
383 * sure that any constraint that is marked redundant (and that is not
384 * also marked zero) is not an equality.
385 * If a constraint is marked DESELECTED, it means the constraint was
386 * SELECTED before (in combination with the same selection of earlier
387 * constraints). If such a deselected constraint turns out to be an
388 * equality, then any vertex that may still be found with the current
389 * selection has already been generated when the constraint was selected.
390 * A constraint is marked UNSELECTED when there is no way selecting
391 * the constraint could lead to a vertex (in combination with the current
392 * selection of earlier constraints).
394 * The set variable coefficients of the selected constraints are stored
395 * in the facets matrix.
397 __isl_give isl_vertices *isl_basic_set_compute_vertices(
398 __isl_keep isl_basic_set *bset)
400 struct isl_tab *tab;
401 int level;
402 int init;
403 unsigned nvar;
404 int *selection = NULL;
405 int selected;
406 struct isl_tab_undo **snap = NULL;
407 isl_mat *facets = NULL;
408 struct isl_vertex_list *list = NULL;
409 int n_vertices = 0;
410 isl_vertices *vertices;
412 if (!bset)
413 return NULL;
415 if (isl_basic_set_plain_is_empty(bset))
416 return vertices_empty(bset);
418 if (bset->n_eq != 0)
419 return lower_dim_vertices(bset);
421 isl_assert(bset->ctx, isl_basic_set_dim(bset, isl_dim_div) == 0,
422 return NULL);
424 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
425 return vertices_0D(bset);
427 nvar = isl_basic_set_dim(bset, isl_dim_set);
429 bset = isl_basic_set_copy(bset);
430 bset = isl_basic_set_set_rational(bset);
431 if (!bset)
432 return NULL;
434 tab = isl_tab_from_basic_set(bset, 0);
435 if (!tab)
436 goto error;
437 tab->strict_redundant = 1;
439 if (tab->empty) {
440 vertices = vertices_empty(bset);
441 isl_basic_set_free(bset);
442 isl_tab_free(tab);
443 return vertices;
446 selection = isl_alloc_array(bset->ctx, int, bset->n_ineq);
447 snap = isl_alloc_array(bset->ctx, struct isl_tab_undo *, bset->n_ineq);
448 facets = isl_mat_alloc(bset->ctx, nvar, nvar);
449 if ((bset->n_ineq && (!selection || !snap)) || !facets)
450 goto error;
452 level = 0;
453 init = 1;
454 selected = 0;
456 while (level >= 0) {
457 if (level >= bset->n_ineq ||
458 (!init && selection[level] != SELECTED)) {
459 --level;
460 init = 0;
461 continue;
463 if (init) {
464 int ok;
465 snap[level] = isl_tab_snap(tab);
466 ok = can_select(bset, level, tab, facets, selected,
467 selection);
468 if (ok < 0)
469 goto error;
470 if (ok) {
471 selection[level] = SELECTED;
472 selected++;
473 } else
474 selection[level] = UNSELECTED;
475 } else {
476 selection[level] = DESELECTED;
477 selected--;
478 if (isl_tab_rollback(tab, snap[level]) < 0)
479 goto error;
481 if (selected == nvar) {
482 if (tab->n_dead == nvar) {
483 isl_bool added = add_vertex(&list, bset, tab);
484 if (added < 0)
485 goto error;
486 if (added)
487 n_vertices++;
489 init = 0;
490 continue;
492 ++level;
493 init = 1;
496 isl_mat_free(facets);
497 free(selection);
498 free(snap);
500 isl_tab_free(tab);
502 vertices = vertices_from_list(bset, n_vertices, list);
504 vertices = compute_chambers(bset, vertices);
506 return vertices;
507 error:
508 free_vertex_list(list);
509 isl_mat_free(facets);
510 free(selection);
511 free(snap);
512 isl_tab_free(tab);
513 isl_basic_set_free(bset);
514 return NULL;
517 struct isl_chamber_list {
518 struct isl_chamber c;
519 struct isl_chamber_list *next;
522 static void free_chamber_list(struct isl_chamber_list *list)
524 struct isl_chamber_list *next;
526 for (; list; list = next) {
527 next = list->next;
528 isl_basic_set_free(list->c.dom);
529 free(list->c.vertices);
530 free(list);
534 /* Check whether the basic set "bset" is a superset of the basic set described
535 * by "tab", i.e., check whether all constraints of "bset" are redundant.
537 static isl_bool bset_covers_tab(__isl_keep isl_basic_set *bset,
538 struct isl_tab *tab)
540 int i;
542 if (!bset || !tab)
543 return isl_bool_error;
545 for (i = 0; i < bset->n_ineq; ++i) {
546 enum isl_ineq_type type = isl_tab_ineq_type(tab, bset->ineq[i]);
547 switch (type) {
548 case isl_ineq_error: return isl_bool_error;
549 case isl_ineq_redundant: continue;
550 default: return isl_bool_false;
554 return isl_bool_true;
557 static __isl_give isl_vertices *vertices_add_chambers(
558 __isl_take isl_vertices *vertices, int n_chambers,
559 struct isl_chamber_list *list)
561 int i;
562 isl_ctx *ctx;
563 struct isl_chamber_list *next;
565 ctx = isl_vertices_get_ctx(vertices);
566 vertices->c = isl_alloc_array(ctx, struct isl_chamber, n_chambers);
567 if (!vertices->c)
568 goto error;
569 vertices->n_chambers = n_chambers;
571 for (i = 0; list; list = next, i++) {
572 next = list->next;
573 vertices->c[i] = list->c;
574 free(list);
577 return vertices;
578 error:
579 isl_vertices_free(vertices);
580 free_chamber_list(list);
581 return NULL;
584 /* Can "tab" be intersected with "bset" without resulting in
585 * a lower-dimensional set.
586 * "bset" itself is assumed to be full-dimensional.
588 static isl_bool can_intersect(struct isl_tab *tab,
589 __isl_keep isl_basic_set *bset)
591 int i;
592 struct isl_tab_undo *snap;
594 if (bset->n_eq > 0)
595 isl_die(isl_basic_set_get_ctx(bset), isl_error_internal,
596 "expecting full-dimensional input",
597 return isl_bool_error);
599 if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
600 return isl_bool_error;
602 snap = isl_tab_snap(tab);
604 for (i = 0; i < bset->n_ineq; ++i) {
605 if (isl_tab_ineq_type(tab, bset->ineq[i]) == isl_ineq_redundant)
606 continue;
607 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
608 return isl_bool_error;
611 if (isl_tab_detect_implicit_equalities(tab) < 0)
612 return isl_bool_error;
613 if (tab->n_dead) {
614 if (isl_tab_rollback(tab, snap) < 0)
615 return isl_bool_error;
616 return isl_bool_false;
619 return isl_bool_true;
622 static int add_chamber(struct isl_chamber_list **list,
623 __isl_keep isl_vertices *vertices, struct isl_tab *tab, int *selection)
625 int n_frozen;
626 int i, j;
627 int n_vertices = 0;
628 struct isl_tab_undo *snap;
629 struct isl_chamber_list *c = NULL;
631 for (i = 0; i < vertices->n_vertices; ++i)
632 if (selection[i])
633 n_vertices++;
635 snap = isl_tab_snap(tab);
637 for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
638 tab->con[i].frozen = 0;
639 n_frozen = i;
641 if (isl_tab_detect_redundant(tab) < 0)
642 return -1;
644 c = isl_calloc_type(tab->mat->ctx, struct isl_chamber_list);
645 if (!c)
646 goto error;
647 c->c.vertices = isl_alloc_array(tab->mat->ctx, int, n_vertices);
648 if (n_vertices && !c->c.vertices)
649 goto error;
650 c->c.dom = isl_basic_set_copy(isl_tab_peek_bset(tab));
651 c->c.dom = isl_basic_set_set_rational(c->c.dom);
652 c->c.dom = isl_basic_set_cow(c->c.dom);
653 c->c.dom = isl_basic_set_update_from_tab(c->c.dom, tab);
654 c->c.dom = isl_basic_set_simplify(c->c.dom);
655 c->c.dom = isl_basic_set_finalize(c->c.dom);
656 if (!c->c.dom)
657 goto error;
659 c->c.n_vertices = n_vertices;
661 for (i = 0, j = 0; i < vertices->n_vertices; ++i)
662 if (selection[i]) {
663 c->c.vertices[j] = i;
664 j++;
667 c->next = *list;
668 *list = c;
670 for (i = 0; i < n_frozen; ++i)
671 tab->con[i].frozen = 1;
673 if (isl_tab_rollback(tab, snap) < 0)
674 return -1;
676 return 0;
677 error:
678 free_chamber_list(c);
679 return -1;
682 struct isl_facet_todo {
683 struct isl_tab *tab; /* A tableau representation of the facet */
684 isl_basic_set *bset; /* A normalized basic set representation */
685 isl_vec *constraint; /* Constraint pointing to the other side */
686 struct isl_facet_todo *next;
689 static void free_todo(struct isl_facet_todo *todo)
691 while (todo) {
692 struct isl_facet_todo *next = todo->next;
694 isl_tab_free(todo->tab);
695 isl_basic_set_free(todo->bset);
696 isl_vec_free(todo->constraint);
697 free(todo);
699 todo = next;
703 static struct isl_facet_todo *create_todo(struct isl_tab *tab, int con)
705 int i;
706 int n_frozen;
707 struct isl_tab_undo *snap;
708 struct isl_facet_todo *todo;
710 snap = isl_tab_snap(tab);
712 for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
713 tab->con[i].frozen = 0;
714 n_frozen = i;
716 if (isl_tab_detect_redundant(tab) < 0)
717 return NULL;
719 todo = isl_calloc_type(tab->mat->ctx, struct isl_facet_todo);
720 if (!todo)
721 return NULL;
723 todo->constraint = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
724 if (!todo->constraint)
725 goto error;
726 isl_seq_neg(todo->constraint->el, tab->bmap->ineq[con], 1 + tab->n_var);
727 todo->bset = isl_basic_set_copy(isl_tab_peek_bset(tab));
728 todo->bset = isl_basic_set_set_rational(todo->bset);
729 todo->bset = isl_basic_set_cow(todo->bset);
730 todo->bset = isl_basic_set_update_from_tab(todo->bset, tab);
731 todo->bset = isl_basic_set_simplify(todo->bset);
732 todo->bset = isl_basic_set_sort_constraints(todo->bset);
733 if (!todo->bset)
734 goto error;
735 ISL_F_SET(todo->bset, ISL_BASIC_SET_NORMALIZED);
736 todo->tab = isl_tab_dup(tab);
737 if (!todo->tab)
738 goto error;
740 for (i = 0; i < n_frozen; ++i)
741 tab->con[i].frozen = 1;
743 if (isl_tab_rollback(tab, snap) < 0)
744 goto error;
746 return todo;
747 error:
748 free_todo(todo);
749 return NULL;
752 /* Create todo items for all interior facets of the chamber represented
753 * by "tab" and collect them in "next".
755 static int init_todo(struct isl_facet_todo **next, struct isl_tab *tab)
757 int i;
758 struct isl_tab_undo *snap;
759 struct isl_facet_todo *todo;
761 snap = isl_tab_snap(tab);
763 for (i = 0; i < tab->n_con; ++i) {
764 if (tab->con[i].frozen)
765 continue;
766 if (tab->con[i].is_redundant)
767 continue;
769 if (isl_tab_select_facet(tab, i) < 0)
770 return -1;
772 todo = create_todo(tab, i);
773 if (!todo)
774 return -1;
776 todo->next = *next;
777 *next = todo;
779 if (isl_tab_rollback(tab, snap) < 0)
780 return -1;
783 return 0;
786 /* Does the linked list contain a todo item that is the opposite of "todo".
787 * If so, return 1 and remove the opposite todo item.
789 static int has_opposite(struct isl_facet_todo *todo,
790 struct isl_facet_todo **list)
792 for (; *list; list = &(*list)->next) {
793 int eq;
794 eq = isl_basic_set_plain_is_equal(todo->bset, (*list)->bset);
795 if (eq < 0)
796 return -1;
797 if (!eq)
798 continue;
799 todo = *list;
800 *list = todo->next;
801 todo->next = NULL;
802 free_todo(todo);
803 return 1;
806 return 0;
809 /* Create todo items for all interior facets of the chamber represented
810 * by "tab" and collect them in first->next, taking care to cancel
811 * opposite todo items.
813 static int update_todo(struct isl_facet_todo *first, struct isl_tab *tab)
815 int i;
816 struct isl_tab_undo *snap;
817 struct isl_facet_todo *todo;
819 snap = isl_tab_snap(tab);
821 for (i = 0; i < tab->n_con; ++i) {
822 int drop;
824 if (tab->con[i].frozen)
825 continue;
826 if (tab->con[i].is_redundant)
827 continue;
829 if (isl_tab_select_facet(tab, i) < 0)
830 return -1;
832 todo = create_todo(tab, i);
833 if (!todo)
834 return -1;
836 drop = has_opposite(todo, &first->next);
837 if (drop < 0)
838 return -1;
840 if (drop)
841 free_todo(todo);
842 else {
843 todo->next = first->next;
844 first->next = todo;
847 if (isl_tab_rollback(tab, snap) < 0)
848 return -1;
851 return 0;
854 /* Compute the chamber decomposition of the parametric polytope respresented
855 * by "bset" given the parametric vertices and their activity domains.
857 * We are only interested in full-dimensional chambers.
858 * Each of these chambers is the intersection of the activity domains of
859 * one or more vertices and the union of all chambers is equal to the
860 * projection of the entire parametric polytope onto the parameter space.
862 * We first create an initial chamber by intersecting as many activity
863 * domains as possible without ending up with an empty or lower-dimensional
864 * set. As a minor optimization, we only consider those activity domains
865 * that contain some arbitrary point.
867 * For each of the interior facets of the chamber, we construct a todo item,
868 * containing the facet and a constraint containing the other side of the facet,
869 * for constructing the chamber on the other side.
870 * While their are any todo items left, we pick a todo item and
871 * create the required chamber by intersecting all activity domains
872 * that contain the facet and have a full-dimensional intersection with
873 * the other side of the facet. For each of the interior facets, we
874 * again create todo items, taking care to cancel opposite todo items.
876 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
877 __isl_take isl_vertices *vertices)
879 int i;
880 isl_ctx *ctx;
881 isl_vec *sample = NULL;
882 struct isl_tab *tab = NULL;
883 struct isl_tab_undo *snap;
884 int *selection = NULL;
885 int n_chambers = 0;
886 struct isl_chamber_list *list = NULL;
887 struct isl_facet_todo *todo = NULL;
889 if (!bset || !vertices)
890 goto error;
892 ctx = isl_vertices_get_ctx(vertices);
893 selection = isl_alloc_array(ctx, int, vertices->n_vertices);
894 if (vertices->n_vertices && !selection)
895 goto error;
897 bset = isl_basic_set_params(bset);
899 tab = isl_tab_from_basic_set(bset, 1);
900 if (!tab)
901 goto error;
902 for (i = 0; i < bset->n_ineq; ++i)
903 if (isl_tab_freeze_constraint(tab, i) < 0)
904 goto error;
905 isl_basic_set_free(bset);
907 snap = isl_tab_snap(tab);
909 sample = isl_tab_get_sample_value(tab);
911 for (i = 0; i < vertices->n_vertices; ++i) {
912 selection[i] = isl_basic_set_contains(vertices->v[i].dom, sample);
913 if (selection[i] < 0)
914 goto error;
915 if (!selection[i])
916 continue;
917 selection[i] = can_intersect(tab, vertices->v[i].dom);
918 if (selection[i] < 0)
919 goto error;
922 if (isl_tab_detect_redundant(tab) < 0)
923 goto error;
925 if (add_chamber(&list, vertices, tab, selection) < 0)
926 goto error;
927 n_chambers++;
929 if (init_todo(&todo, tab) < 0)
930 goto error;
932 while (todo) {
933 struct isl_facet_todo *next;
935 if (isl_tab_rollback(tab, snap) < 0)
936 goto error;
938 if (isl_tab_add_ineq(tab, todo->constraint->el) < 0)
939 goto error;
940 if (isl_tab_freeze_constraint(tab, tab->n_con - 1) < 0)
941 goto error;
943 for (i = 0; i < vertices->n_vertices; ++i) {
944 selection[i] = bset_covers_tab(vertices->v[i].dom,
945 todo->tab);
946 if (selection[i] < 0)
947 goto error;
948 if (!selection[i])
949 continue;
950 selection[i] = can_intersect(tab, vertices->v[i].dom);
951 if (selection[i] < 0)
952 goto error;
955 if (isl_tab_detect_redundant(tab) < 0)
956 goto error;
958 if (add_chamber(&list, vertices, tab, selection) < 0)
959 goto error;
960 n_chambers++;
962 if (update_todo(todo, tab) < 0)
963 goto error;
965 next = todo->next;
966 todo->next = NULL;
967 free_todo(todo);
968 todo = next;
971 isl_vec_free(sample);
973 isl_tab_free(tab);
974 free(selection);
976 vertices = vertices_add_chambers(vertices, n_chambers, list);
978 for (i = 0; vertices && i < vertices->n_vertices; ++i) {
979 isl_basic_set_free(vertices->v[i].dom);
980 vertices->v[i].dom = NULL;
983 return vertices;
984 error:
985 free_chamber_list(list);
986 free_todo(todo);
987 isl_vec_free(sample);
988 isl_tab_free(tab);
989 free(selection);
990 if (!tab)
991 isl_basic_set_free(bset);
992 isl_vertices_free(vertices);
993 return NULL;
996 isl_ctx *isl_vertex_get_ctx(__isl_keep isl_vertex *vertex)
998 return vertex ? isl_vertices_get_ctx(vertex->vertices) : NULL;
1001 int isl_vertex_get_id(__isl_keep isl_vertex *vertex)
1003 return vertex ? vertex->id : -1;
1006 /* Return the activity domain of the vertex "vertex".
1008 __isl_give isl_basic_set *isl_vertex_get_domain(__isl_keep isl_vertex *vertex)
1010 struct isl_vertex *v;
1012 if (!vertex)
1013 return NULL;
1015 v = &vertex->vertices->v[vertex->id];
1016 if (!v->dom) {
1017 v->dom = isl_basic_set_copy(v->vertex);
1018 v->dom = isl_basic_set_params(v->dom);
1019 v->dom = isl_basic_set_set_integral(v->dom);
1022 return isl_basic_set_copy(v->dom);
1025 /* Return a multiple quasi-affine expression describing the vertex "vertex"
1026 * in terms of the parameters,
1028 __isl_give isl_multi_aff *isl_vertex_get_expr(__isl_keep isl_vertex *vertex)
1030 struct isl_vertex *v;
1031 isl_basic_set *bset;
1033 if (!vertex)
1034 return NULL;
1036 v = &vertex->vertices->v[vertex->id];
1038 bset = isl_basic_set_copy(v->vertex);
1039 return isl_multi_aff_from_basic_set_equalities(bset);
1042 static __isl_give isl_vertex *isl_vertex_alloc(__isl_take isl_vertices *vertices,
1043 int id)
1045 isl_ctx *ctx;
1046 isl_vertex *vertex;
1048 if (!vertices)
1049 return NULL;
1051 ctx = isl_vertices_get_ctx(vertices);
1052 vertex = isl_alloc_type(ctx, isl_vertex);
1053 if (!vertex)
1054 goto error;
1056 vertex->vertices = vertices;
1057 vertex->id = id;
1059 return vertex;
1060 error:
1061 isl_vertices_free(vertices);
1062 return NULL;
1065 void isl_vertex_free(__isl_take isl_vertex *vertex)
1067 if (!vertex)
1068 return;
1069 isl_vertices_free(vertex->vertices);
1070 free(vertex);
1073 isl_ctx *isl_cell_get_ctx(__isl_keep isl_cell *cell)
1075 return cell ? cell->dom->ctx : NULL;
1078 __isl_give isl_basic_set *isl_cell_get_domain(__isl_keep isl_cell *cell)
1080 return cell ? isl_basic_set_copy(cell->dom) : NULL;
1083 static __isl_give isl_cell *isl_cell_alloc(__isl_take isl_vertices *vertices,
1084 __isl_take isl_basic_set *dom, int id)
1086 int i;
1087 isl_cell *cell = NULL;
1089 if (!vertices || !dom)
1090 goto error;
1092 cell = isl_calloc_type(dom->ctx, isl_cell);
1093 if (!cell)
1094 goto error;
1096 cell->n_vertices = vertices->c[id].n_vertices;
1097 cell->ids = isl_alloc_array(dom->ctx, int, cell->n_vertices);
1098 if (cell->n_vertices && !cell->ids)
1099 goto error;
1100 for (i = 0; i < cell->n_vertices; ++i)
1101 cell->ids[i] = vertices->c[id].vertices[i];
1102 cell->vertices = vertices;
1103 cell->dom = dom;
1105 return cell;
1106 error:
1107 isl_cell_free(cell);
1108 isl_vertices_free(vertices);
1109 isl_basic_set_free(dom);
1110 return NULL;
1113 void isl_cell_free(__isl_take isl_cell *cell)
1115 if (!cell)
1116 return;
1118 isl_vertices_free(cell->vertices);
1119 free(cell->ids);
1120 isl_basic_set_free(cell->dom);
1121 free(cell);
1124 /* Create a tableau of the cone obtained by first homogenizing the given
1125 * polytope and then making all inequalities strict by setting the
1126 * constant term to -1.
1128 static struct isl_tab *tab_for_shifted_cone(__isl_keep isl_basic_set *bset)
1130 int i;
1131 isl_vec *c = NULL;
1132 struct isl_tab *tab;
1134 if (!bset)
1135 return NULL;
1136 tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq + 1,
1137 1 + isl_basic_set_total_dim(bset), 0);
1138 if (!tab)
1139 return NULL;
1140 tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1141 if (ISL_F_ISSET(bset, ISL_BASIC_MAP_EMPTY)) {
1142 if (isl_tab_mark_empty(tab) < 0)
1143 goto error;
1144 return tab;
1147 c = isl_vec_alloc(bset->ctx, 1 + 1 + isl_basic_set_total_dim(bset));
1148 if (!c)
1149 goto error;
1151 isl_int_set_si(c->el[0], 0);
1152 for (i = 0; i < bset->n_eq; ++i) {
1153 isl_seq_cpy(c->el + 1, bset->eq[i], c->size - 1);
1154 if (isl_tab_add_eq(tab, c->el) < 0)
1155 goto error;
1158 isl_int_set_si(c->el[0], -1);
1159 for (i = 0; i < bset->n_ineq; ++i) {
1160 isl_seq_cpy(c->el + 1, bset->ineq[i], c->size - 1);
1161 if (isl_tab_add_ineq(tab, c->el) < 0)
1162 goto error;
1163 if (tab->empty) {
1164 isl_vec_free(c);
1165 return tab;
1169 isl_seq_clr(c->el + 1, c->size - 1);
1170 isl_int_set_si(c->el[1], 1);
1171 if (isl_tab_add_ineq(tab, c->el) < 0)
1172 goto error;
1174 isl_vec_free(c);
1175 return tab;
1176 error:
1177 isl_vec_free(c);
1178 isl_tab_free(tab);
1179 return NULL;
1182 /* Compute an interior point of "bset" by selecting an interior
1183 * point in homogeneous space and projecting the point back down.
1185 static __isl_give isl_vec *isl_basic_set_interior_point(
1186 __isl_keep isl_basic_set *bset)
1188 isl_vec *vec;
1189 struct isl_tab *tab;
1191 tab = tab_for_shifted_cone(bset);
1192 vec = isl_tab_get_sample_value(tab);
1193 isl_tab_free(tab);
1194 if (!vec)
1195 return NULL;
1197 isl_seq_cpy(vec->el, vec->el + 1, vec->size - 1);
1198 vec->size--;
1200 return vec;
1203 /* Call "fn" on all chambers of the parametric polytope with the shared
1204 * facets of neighboring chambers only appearing in one of the chambers.
1206 * We pick an interior point from one of the chambers and then make
1207 * all constraints that do not satisfy this point strict.
1208 * For constraints that saturate the interior point, the sign
1209 * of the first non-zero coefficient is used to determine which
1210 * of the two (internal) constraints should be tightened.
1212 isl_stat isl_vertices_foreach_disjoint_cell(__isl_keep isl_vertices *vertices,
1213 isl_stat (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1215 int i;
1216 isl_vec *vec;
1217 isl_cell *cell;
1219 if (!vertices)
1220 return isl_stat_error;
1222 if (vertices->n_chambers == 0)
1223 return isl_stat_ok;
1225 if (vertices->n_chambers == 1) {
1226 isl_basic_set *dom = isl_basic_set_copy(vertices->c[0].dom);
1227 dom = isl_basic_set_set_integral(dom);
1228 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, 0);
1229 if (!cell)
1230 return isl_stat_error;
1231 return fn(cell, user);
1234 vec = isl_basic_set_interior_point(vertices->c[0].dom);
1235 if (!vec)
1236 return isl_stat_error;
1238 for (i = 0; i < vertices->n_chambers; ++i) {
1239 int r;
1240 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1241 if (i)
1242 dom = isl_basic_set_tighten_outward(dom, vec);
1243 dom = isl_basic_set_set_integral(dom);
1244 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1245 if (!cell)
1246 goto error;
1247 r = fn(cell, user);
1248 if (r < 0)
1249 goto error;
1252 isl_vec_free(vec);
1254 return isl_stat_ok;
1255 error:
1256 isl_vec_free(vec);
1257 return isl_stat_error;
1260 isl_stat isl_vertices_foreach_cell(__isl_keep isl_vertices *vertices,
1261 isl_stat (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1263 int i;
1264 isl_cell *cell;
1266 if (!vertices)
1267 return isl_stat_error;
1269 if (vertices->n_chambers == 0)
1270 return isl_stat_ok;
1272 for (i = 0; i < vertices->n_chambers; ++i) {
1273 isl_stat r;
1274 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1276 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1277 if (!cell)
1278 return isl_stat_error;
1280 r = fn(cell, user);
1281 if (r < 0)
1282 return isl_stat_error;
1285 return isl_stat_ok;
1288 isl_stat isl_vertices_foreach_vertex(__isl_keep isl_vertices *vertices,
1289 isl_stat (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1291 int i;
1292 isl_vertex *vertex;
1294 if (!vertices)
1295 return isl_stat_error;
1297 if (vertices->n_vertices == 0)
1298 return isl_stat_ok;
1300 for (i = 0; i < vertices->n_vertices; ++i) {
1301 isl_stat r;
1303 vertex = isl_vertex_alloc(isl_vertices_copy(vertices), i);
1304 if (!vertex)
1305 return isl_stat_error;
1307 r = fn(vertex, user);
1308 if (r < 0)
1309 return isl_stat_error;
1312 return isl_stat_ok;
1315 isl_stat isl_cell_foreach_vertex(__isl_keep isl_cell *cell,
1316 isl_stat (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1318 int i;
1319 isl_vertex *vertex;
1321 if (!cell)
1322 return isl_stat_error;
1324 if (cell->n_vertices == 0)
1325 return isl_stat_ok;
1327 for (i = 0; i < cell->n_vertices; ++i) {
1328 isl_stat r;
1330 vertex = isl_vertex_alloc(isl_vertices_copy(cell->vertices),
1331 cell->ids[i]);
1332 if (!vertex)
1333 return isl_stat_error;
1335 r = fn(vertex, user);
1336 if (r < 0)
1337 return isl_stat_error;
1340 return isl_stat_ok;
1343 isl_ctx *isl_vertices_get_ctx(__isl_keep isl_vertices *vertices)
1345 return vertices ? vertices->bset->ctx : NULL;
1348 int isl_vertices_get_n_vertices(__isl_keep isl_vertices *vertices)
1350 return vertices ? vertices->n_vertices : -1;
1353 __isl_give isl_vertices *isl_morph_vertices(__isl_take isl_morph *morph,
1354 __isl_take isl_vertices *vertices)
1356 int i;
1357 isl_morph *param_morph = NULL;
1359 if (!morph || !vertices)
1360 goto error;
1362 isl_assert(vertices->bset->ctx, vertices->ref == 1, goto error);
1364 param_morph = isl_morph_copy(morph);
1365 param_morph = isl_morph_dom_params(param_morph);
1366 param_morph = isl_morph_ran_params(param_morph);
1368 for (i = 0; i < vertices->n_vertices; ++i) {
1369 vertices->v[i].dom = isl_morph_basic_set(
1370 isl_morph_copy(param_morph), vertices->v[i].dom);
1371 vertices->v[i].vertex = isl_morph_basic_set(
1372 isl_morph_copy(morph), vertices->v[i].vertex);
1373 if (!vertices->v[i].vertex)
1374 goto error;
1377 for (i = 0; i < vertices->n_chambers; ++i) {
1378 vertices->c[i].dom = isl_morph_basic_set(
1379 isl_morph_copy(param_morph), vertices->c[i].dom);
1380 if (!vertices->c[i].dom)
1381 goto error;
1384 isl_morph_free(param_morph);
1385 isl_morph_free(morph);
1386 return vertices;
1387 error:
1388 isl_morph_free(param_morph);
1389 isl_morph_free(morph);
1390 isl_vertices_free(vertices);
1391 return NULL;
1394 /* Construct a simplex isl_cell spanned by the vertices with indices in
1395 * "simplex_ids" and "other_ids" and call "fn" on this isl_cell.
1397 static isl_stat call_on_simplex(__isl_keep isl_cell *cell,
1398 int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1399 isl_stat (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1401 int i;
1402 isl_ctx *ctx;
1403 struct isl_cell *simplex;
1405 ctx = isl_cell_get_ctx(cell);
1407 simplex = isl_calloc_type(ctx, struct isl_cell);
1408 if (!simplex)
1409 return isl_stat_error;
1410 simplex->vertices = isl_vertices_copy(cell->vertices);
1411 if (!simplex->vertices)
1412 goto error;
1413 simplex->dom = isl_basic_set_copy(cell->dom);
1414 if (!simplex->dom)
1415 goto error;
1416 simplex->n_vertices = n_simplex + n_other;
1417 simplex->ids = isl_alloc_array(ctx, int, simplex->n_vertices);
1418 if (!simplex->ids)
1419 goto error;
1421 for (i = 0; i < n_simplex; ++i)
1422 simplex->ids[i] = simplex_ids[i];
1423 for (i = 0; i < n_other; ++i)
1424 simplex->ids[n_simplex + i] = other_ids[i];
1426 return fn(simplex, user);
1427 error:
1428 isl_cell_free(simplex);
1429 return isl_stat_error;
1432 /* Check whether the parametric vertex described by "vertex"
1433 * lies on the facet corresponding to constraint "facet" of "bset".
1434 * The isl_vec "v" is a temporary vector than can be used by this function.
1436 * We eliminate the variables from the facet constraint using the
1437 * equalities defining the vertex and check if the result is identical
1438 * to zero.
1440 * It would probably be better to keep track of the constraints defining
1441 * a vertex during the vertex construction so that we could simply look
1442 * it up here.
1444 static int vertex_on_facet(__isl_keep isl_basic_set *vertex,
1445 __isl_keep isl_basic_set *bset, int facet, __isl_keep isl_vec *v)
1447 int i;
1448 isl_int m;
1450 isl_seq_cpy(v->el, bset->ineq[facet], v->size);
1452 isl_int_init(m);
1453 for (i = 0; i < vertex->n_eq; ++i) {
1454 int k = isl_seq_last_non_zero(vertex->eq[i], v->size);
1455 isl_seq_elim(v->el, vertex->eq[i], k, v->size, &m);
1457 isl_int_clear(m);
1459 return isl_seq_first_non_zero(v->el, v->size) == -1;
1462 /* Triangulate the polytope spanned by the vertices with ids
1463 * in "simplex_ids" and "other_ids" and call "fn" on each of
1464 * the resulting simplices.
1465 * If the input polytope is already a simplex, we simply call "fn".
1466 * Otherwise, we pick a point from "other_ids" and add it to "simplex_ids".
1467 * Then we consider each facet of "bset" that does not contain the point
1468 * we just picked, but does contain some of the other points in "other_ids"
1469 * and call ourselves recursively on the polytope spanned by the new
1470 * "simplex_ids" and those points in "other_ids" that lie on the facet.
1472 static isl_stat triangulate(__isl_keep isl_cell *cell, __isl_keep isl_vec *v,
1473 int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1474 isl_stat (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1476 int i, j, k;
1477 int d, nparam;
1478 int *ids;
1479 isl_ctx *ctx;
1480 isl_basic_set *vertex;
1481 isl_basic_set *bset;
1483 ctx = isl_cell_get_ctx(cell);
1484 d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1485 nparam = isl_basic_set_dim(cell->vertices->bset, isl_dim_param);
1487 if (n_simplex + n_other == d + 1)
1488 return call_on_simplex(cell, simplex_ids, n_simplex,
1489 other_ids, n_other, fn, user);
1491 simplex_ids[n_simplex] = other_ids[0];
1492 vertex = cell->vertices->v[other_ids[0]].vertex;
1493 bset = cell->vertices->bset;
1495 ids = isl_alloc_array(ctx, int, n_other - 1);
1496 if (!ids)
1497 goto error;
1498 for (i = 0; i < bset->n_ineq; ++i) {
1499 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + nparam, d) == -1)
1500 continue;
1501 if (vertex_on_facet(vertex, bset, i, v))
1502 continue;
1504 for (j = 1, k = 0; j < n_other; ++j) {
1505 isl_basic_set *ov;
1506 ov = cell->vertices->v[other_ids[j]].vertex;
1507 if (vertex_on_facet(ov, bset, i, v))
1508 ids[k++] = other_ids[j];
1510 if (k == 0)
1511 continue;
1513 if (triangulate(cell, v, simplex_ids, n_simplex + 1,
1514 ids, k, fn, user) < 0)
1515 goto error;
1517 free(ids);
1519 return isl_stat_ok;
1520 error:
1521 free(ids);
1522 return isl_stat_error;
1525 /* Triangulate the given cell and call "fn" on each of the resulting
1526 * simplices.
1528 isl_stat isl_cell_foreach_simplex(__isl_take isl_cell *cell,
1529 isl_stat (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1531 int d, total;
1532 int r;
1533 isl_ctx *ctx;
1534 isl_vec *v = NULL;
1535 int *simplex_ids = NULL;
1537 if (!cell)
1538 return isl_stat_error;
1540 d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1541 total = isl_basic_set_total_dim(cell->vertices->bset);
1543 if (cell->n_vertices == d + 1)
1544 return fn(cell, user);
1546 ctx = isl_cell_get_ctx(cell);
1547 simplex_ids = isl_alloc_array(ctx, int, d + 1);
1548 if (!simplex_ids)
1549 goto error;
1551 v = isl_vec_alloc(ctx, 1 + total);
1552 if (!v)
1553 goto error;
1555 r = triangulate(cell, v, simplex_ids, 0,
1556 cell->ids, cell->n_vertices, fn, user);
1558 isl_vec_free(v);
1559 free(simplex_ids);
1561 isl_cell_free(cell);
1563 return r;
1564 error:
1565 free(simplex_ids);
1566 isl_vec_free(v);
1567 isl_cell_free(cell);
1568 return isl_stat_error;