add isl_vec_sort
[isl.git] / isl_vertices.c
blob1dd0f03bdb04714a504950c5754994295d4ea4ce
1 /*
2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 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/set.h>
13 #include <isl/seq.h>
14 #include <isl_tab.h>
15 #include <isl_dim_private.h>
16 #include <isl_morph.h>
17 #include <isl_vertices_private.h>
18 #include <isl_mat_private.h>
20 #define SELECTED 1
21 #define DESELECTED -1
22 #define UNSELECTED 0
24 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
25 __isl_take isl_vertices *vertices);
27 __isl_give isl_vertices *isl_vertices_copy(__isl_keep isl_vertices *vertices)
29 if (!vertices)
30 return NULL;
32 vertices->ref++;
33 return vertices;
36 void isl_vertices_free(__isl_take isl_vertices *vertices)
38 int i;
40 if (!vertices)
41 return;
43 if (--vertices->ref > 0)
44 return;
46 for (i = 0; i < vertices->n_vertices; ++i) {
47 isl_basic_set_free(vertices->v[i].vertex);
48 isl_basic_set_free(vertices->v[i].dom);
50 free(vertices->v);
52 for (i = 0; i < vertices->n_chambers; ++i) {
53 free(vertices->c[i].vertices);
54 isl_basic_set_free(vertices->c[i].dom);
56 free(vertices->c);
58 isl_basic_set_free(vertices->bset);
59 free(vertices);
62 struct isl_vertex_list {
63 struct isl_vertex v;
64 struct isl_vertex_list *next;
67 static void free_vertex_list(struct isl_vertex_list *list)
69 struct isl_vertex_list *next;
71 for (; list; list = next) {
72 next = list->next;
73 isl_basic_set_free(list->v.vertex);
74 isl_basic_set_free(list->v.dom);
75 free(list);
79 static __isl_give isl_vertices *vertices_from_list(__isl_keep isl_basic_set *bset,
80 int n_vertices, struct isl_vertex_list *list)
82 int i;
83 struct isl_vertex_list *next;
84 isl_vertices *vertices;
86 vertices = isl_calloc_type(bset->ctx, isl_vertices);
87 if (!vertices)
88 goto error;
89 vertices->ref = 1;
90 vertices->bset = isl_basic_set_copy(bset);
91 vertices->v = isl_alloc_array(bset->ctx, struct isl_vertex, n_vertices);
92 if (!vertices->v)
93 goto error;
94 vertices->n_vertices = n_vertices;
96 for (i = 0; list; list = next, i++) {
97 next = list->next;
98 vertices->v[i] = list->v;
99 free(list);
102 return vertices;
103 error:
104 free(vertices);
105 free_vertex_list(list);
106 return NULL;
109 /* Prepend a vertex to the linked list "list" based on the equalities in "tab".
111 static int add_vertex(struct isl_vertex_list **list,
112 __isl_keep isl_basic_set *bset, struct isl_tab *tab)
114 unsigned nvar;
115 unsigned nparam;
116 struct isl_vertex_list *v = NULL;
118 if (isl_tab_detect_implicit_equalities(tab) < 0)
119 return -1;
121 nvar = isl_basic_set_dim(bset, isl_dim_set);
122 nparam = isl_basic_set_dim(bset, isl_dim_param);
124 v = isl_calloc_type(tab->mat->ctx, struct isl_vertex_list);
125 if (!v)
126 goto error;
128 v->v.vertex = isl_basic_set_copy(bset);
129 v->v.vertex = isl_basic_set_cow(v->v.vertex);
130 v->v.vertex = isl_basic_set_update_from_tab(v->v.vertex, tab);
131 v->v.vertex = isl_basic_set_simplify(v->v.vertex);
132 v->v.vertex = isl_basic_set_finalize(v->v.vertex);
133 if (!v->v.vertex)
134 goto error;
135 isl_assert(bset->ctx, v->v.vertex->n_eq >= nvar, goto error);
136 v->v.dom = isl_basic_set_copy(v->v.vertex);
137 v->v.dom = isl_basic_set_project_out(v->v.dom, isl_dim_set, 0, nvar);
138 if (!v->v.dom)
139 goto error;
141 v->next = *list;
142 *list = v;
144 return 0;
145 error:
146 free_vertex_list(v);
147 return -1;
150 /* Compute the parametric vertices and the chamber decomposition
151 * of an empty parametric polytope.
153 static __isl_give isl_vertices *vertices_empty(__isl_keep isl_basic_set *bset)
155 isl_vertices *vertices;
156 unsigned nparam;
158 if (!bset)
159 return NULL;
161 nparam = isl_basic_set_dim(bset, isl_dim_param);
163 vertices = isl_calloc_type(bset->ctx, isl_vertices);
164 if (!vertices)
165 return NULL;
166 vertices->bset = isl_basic_set_copy(bset);
167 vertices->ref = 1;
169 vertices->n_vertices = 0;
170 vertices->n_chambers = 0;
172 return vertices;
175 /* Compute the parametric vertices and the chamber decomposition
176 * of the parametric polytope defined using the same constraints
177 * as "bset" in the 0D case.
178 * There is exactly one 0D vertex and a single chamber containing
179 * the vertex.
181 static __isl_give isl_vertices *vertices_0D(__isl_keep isl_basic_set *bset)
183 isl_vertices *vertices;
184 unsigned nparam;
186 if (!bset)
187 return NULL;
189 nparam = isl_basic_set_dim(bset, isl_dim_param);
191 vertices = isl_calloc_type(bset->ctx, isl_vertices);
192 if (!vertices)
193 return NULL;
194 vertices->ref = 1;
195 vertices->bset = isl_basic_set_copy(bset);
197 vertices->v = isl_calloc_array(bset->ctx, struct isl_vertex, 1);
198 if (!vertices->v)
199 goto error;
200 vertices->n_vertices = 1;
201 vertices->v[0].vertex = isl_basic_set_copy(bset);
202 if (!vertices->v[0].vertex)
203 goto error;
205 vertices->c = isl_calloc_array(bset->ctx, struct isl_chamber, 1);
206 if (!vertices->c)
207 goto error;
208 vertices->n_chambers = 1;
209 vertices->c[0].n_vertices = 1;
210 vertices->c[0].vertices = isl_calloc_array(bset->ctx, int, 1);
211 if (!vertices->c[0].vertices)
212 goto error;
213 vertices->c[0].dom = isl_basic_set_copy(bset);
214 if (!vertices->c[0].dom)
215 goto error;
217 return vertices;
218 error:
219 isl_vertices_free(vertices);
220 return NULL;
223 static int isl_mat_rank(__isl_keep isl_mat *mat)
225 int row, col;
226 isl_mat *H;
228 H = isl_mat_left_hermite(isl_mat_copy(mat), 0, NULL, NULL);
229 if (!H)
230 return -1;
232 for (col = 0; col < H->n_col; ++col) {
233 for (row = 0; row < H->n_row; ++row)
234 if (!isl_int_is_zero(H->row[row][col]))
235 break;
236 if (row == H->n_row)
237 break;
240 isl_mat_free(H);
242 return col;
245 /* Is the row pointed to by "f" linearly independent of the "n" first
246 * rows in "facets"?
248 static int is_independent(__isl_keep isl_mat *facets, int n, isl_int *f)
250 int rank;
252 if (isl_seq_first_non_zero(f, facets->n_col) < 0)
253 return 0;
255 isl_seq_cpy(facets->row[n], f, facets->n_col);
256 facets->n_row = n + 1;
257 rank = isl_mat_rank(facets);
258 if (rank < 0)
259 return -1;
261 return rank == n + 1;
264 /* Check whether we can select constraint "level", given the current selection
265 * reflected by facets in "tab", the rows of "facets" and the earlier
266 * "selected" elements of "selection".
268 * If the constraint is (strictly) redundant in the tableau, selecting it would
269 * result in an empty tableau, so it can't be selected.
270 * If the set variable part of the constraint is not linearly indepedent
271 * of the set variable parts of the already selected constraints,
272 * the constraint cannot be selected.
273 * If selecting the constraint results in an empty tableau, the constraint
274 * cannot be selected.
275 * Finally, if selecting the constraint results in some explicitly
276 * deselected constraints turning into equalities, then the corresponding
277 * vertices have already been generated, so the constraint cannot be selected.
279 static int can_select(__isl_keep isl_basic_set *bset, int level,
280 struct isl_tab *tab, __isl_keep isl_mat *facets, int selected,
281 int *selection)
283 int i;
284 int indep;
285 unsigned ovar;
286 struct isl_tab_undo *snap;
288 if (isl_tab_is_redundant(tab, level))
289 return 0;
291 ovar = isl_dim_offset(bset->dim, isl_dim_set);
293 indep = is_independent(facets, selected, bset->ineq[level] + 1 + ovar);
294 if (indep < 0)
295 return -1;
296 if (!indep)
297 return 0;
299 snap = isl_tab_snap(tab);
300 if (isl_tab_select_facet(tab, level) < 0)
301 return -1;
303 if (tab->empty) {
304 if (isl_tab_rollback(tab, snap) < 0)
305 return -1;
306 return 0;
309 for (i = 0; i < level; ++i) {
310 int sgn;
312 if (selection[i] != DESELECTED)
313 continue;
315 if (isl_tab_is_equality(tab, i))
316 sgn = 0;
317 else if (isl_tab_is_redundant(tab, i))
318 sgn = 1;
319 else
320 sgn = isl_tab_sign_of_max(tab, i);
321 if (sgn < -1)
322 return -1;
323 if (sgn <= 0) {
324 if (isl_tab_rollback(tab, snap) < 0)
325 return -1;
326 return 0;
330 return 1;
333 /* Compute the parametric vertices and the chamber decomposition
334 * of a parametric polytope that is not full-dimensional.
336 * Simply map the parametric polytope to a lower dimensional space
337 * and map the resulting vertices back.
339 static __isl_give isl_vertices *lower_dim_vertices(
340 __isl_keep isl_basic_set *bset)
342 isl_morph *morph;
343 isl_vertices *vertices;
345 bset = isl_basic_set_copy(bset);
346 morph = isl_basic_set_full_compression(bset);
347 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
349 vertices = isl_basic_set_compute_vertices(bset);
350 isl_basic_set_free(bset);
352 morph = isl_morph_inverse(morph);
354 vertices = isl_morph_vertices(morph, vertices);
356 return vertices;
359 /* Compute the parametric vertices and the chamber decomposition
360 * of the parametric polytope defined using the same constraints
361 * as "bset". "bset" is assumed to have no existentially quantified
362 * variables.
364 * The vertices themselves are computed in a fairly simplistic way.
365 * We simply run through all combinations of d constraints,
366 * with d the number of set variables, and check if those d constraints
367 * define a vertex. To avoid the generation of duplicate vertices,
368 * which we may happen if a vertex is defined by more that d constraints,
369 * we make sure we only generate the vertex for the d constraints with
370 * smallest index.
372 * We set up a tableau and keep track of which facets have been
373 * selected. The tableau is marked strict_redundant so that we can be
374 * sure that any constraint that is marked redundant (and that is not
375 * also marked zero) is not an equality.
376 * If a constraint is marked DESELECTED, it means the constraint was
377 * SELECTED before (in combination with the same selection of earlier
378 * constraints). If such a deselected constraint turns out to be an
379 * equality, then any vertex that may still be found with the current
380 * selection has already been generated when the constraint was selected.
381 * A constraint is marked UNSELECTED when there is no way selecting
382 * the constraint could lead to a vertex (in combination with the current
383 * selection of earlier constraints).
385 * The set variable coefficients of the selected constraints are stored
386 * in the facets matrix.
388 __isl_give isl_vertices *isl_basic_set_compute_vertices(
389 __isl_keep isl_basic_set *bset)
391 struct isl_tab *tab;
392 int level;
393 int init;
394 unsigned nvar;
395 int *selection;
396 int selected;
397 struct isl_tab_undo **snap;
398 isl_mat *facets;
399 struct isl_vertex_list *list = NULL;
400 int n_vertices = 0;
401 isl_vertices *vertices;
403 if (!bset)
404 return NULL;
406 if (isl_basic_set_plain_is_empty(bset))
407 return vertices_empty(bset);
409 if (bset->n_eq != 0)
410 return lower_dim_vertices(bset);
412 isl_assert(bset->ctx, isl_basic_set_dim(bset, isl_dim_div) == 0,
413 return NULL);
415 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
416 return vertices_0D(bset);
418 nvar = isl_basic_set_dim(bset, isl_dim_set);
420 bset = isl_basic_set_copy(bset);
421 bset = isl_basic_set_set_rational(bset);
422 if (!bset)
423 return NULL;
425 tab = isl_tab_from_basic_set(bset);
426 if (!tab)
427 goto error;
428 tab->strict_redundant = 1;
430 if (tab->empty) {
431 vertices = vertices_empty(bset);
432 isl_basic_set_free(bset);
433 isl_tab_free(tab);
434 return vertices;
437 selection = isl_alloc_array(bset->ctx, int, bset->n_ineq);
438 snap = isl_alloc_array(bset->ctx, struct isl_tab_undo *, bset->n_ineq);
439 facets = isl_mat_alloc(bset->ctx, nvar, nvar);
440 if (!selection || !snap || !facets)
441 goto error;
443 level = 0;
444 init = 1;
445 selected = 0;
447 while (level >= 0) {
448 if (level >= bset->n_ineq ||
449 (!init && selection[level] != SELECTED)) {
450 --level;
451 init = 0;
452 continue;
454 if (init) {
455 int ok;
456 snap[level] = isl_tab_snap(tab);
457 ok = can_select(bset, level, tab, facets, selected,
458 selection);
459 if (ok < 0)
460 goto error;
461 if (ok) {
462 selection[level] = SELECTED;
463 selected++;
464 } else
465 selection[level] = UNSELECTED;
466 } else {
467 selection[level] = DESELECTED;
468 selected--;
469 if (isl_tab_rollback(tab, snap[level]) < 0)
470 goto error;
472 if (selected == nvar) {
473 if (tab->n_dead == nvar) {
474 if (add_vertex(&list, bset, tab) < 0)
475 goto error;
476 n_vertices++;
478 init = 0;
479 continue;
481 ++level;
482 init = 1;
485 isl_mat_free(facets);
486 free(selection);
487 free(snap);
489 isl_tab_free(tab);
491 vertices = vertices_from_list(bset, n_vertices, list);
493 vertices = compute_chambers(bset, vertices);
495 return vertices;
496 error:
497 isl_mat_free(facets);
498 free(selection);
499 free(snap);
500 isl_tab_free(tab);
501 isl_basic_set_free(bset);
502 return NULL;
505 struct isl_chamber_list {
506 struct isl_chamber c;
507 struct isl_chamber_list *next;
510 static void free_chamber_list(struct isl_chamber_list *list)
512 struct isl_chamber_list *next;
514 for (; list; list = next) {
515 next = list->next;
516 isl_basic_set_free(list->c.dom);
517 free(list->c.vertices);
518 free(list);
522 /* Check whether the basic set "bset" is a superset of the basic set described
523 * by "tab", i.e., check whether all constraints of "bset" are redundant.
525 static int bset_covers_tab(__isl_keep isl_basic_set *bset, struct isl_tab *tab)
527 int i;
529 if (!bset || !tab)
530 return -1;
532 for (i = 0; i < bset->n_ineq; ++i) {
533 enum isl_ineq_type type = isl_tab_ineq_type(tab, bset->ineq[i]);
534 switch (type) {
535 case isl_ineq_error: return -1;
536 case isl_ineq_redundant: continue;
537 default: return 0;
541 return 1;
544 static __isl_give isl_vertices *vertices_add_chambers(
545 __isl_take isl_vertices *vertices, int n_chambers,
546 struct isl_chamber_list *list)
548 int i;
549 isl_ctx *ctx;
550 struct isl_chamber_list *next;
552 ctx = isl_vertices_get_ctx(vertices);
553 vertices->c = isl_alloc_array(ctx, struct isl_chamber, n_chambers);
554 if (!vertices->c)
555 goto error;
556 vertices->n_chambers = n_chambers;
558 for (i = 0; list; list = next, i++) {
559 next = list->next;
560 vertices->c[i] = list->c;
561 free(list);
564 return vertices;
565 error:
566 isl_vertices_free(vertices);
567 free_chamber_list(list);
568 return NULL;
571 /* Can "tab" be intersected with "bset" without resulting in
572 * a lower-dimensional set.
574 static int can_intersect(struct isl_tab *tab, __isl_keep isl_basic_set *bset)
576 int i;
577 struct isl_tab_undo *snap;
579 if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
580 return -1;
582 snap = isl_tab_snap(tab);
584 for (i = 0; i < bset->n_ineq; ++i) {
585 if (isl_tab_ineq_type(tab, bset->ineq[i]) == isl_ineq_redundant)
586 continue;
587 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
588 return -1;
591 if (isl_tab_detect_implicit_equalities(tab) < 0)
592 return -1;
593 if (tab->n_dead) {
594 if (isl_tab_rollback(tab, snap) < 0)
595 return -1;
596 return 0;
599 return 1;
602 static int add_chamber(struct isl_chamber_list **list,
603 __isl_keep isl_vertices *vertices, struct isl_tab *tab, int *selection)
605 int n_frozen;
606 int i, j;
607 int n_vertices = 0;
608 struct isl_tab_undo *snap;
609 struct isl_chamber_list *c = NULL;
611 for (i = 0; i < vertices->n_vertices; ++i)
612 if (selection[i])
613 n_vertices++;
615 snap = isl_tab_snap(tab);
617 for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
618 tab->con[i].frozen = 0;
619 n_frozen = i;
621 if (isl_tab_detect_redundant(tab) < 0)
622 return -1;
624 c = isl_calloc_type(tab->mat->ctx, struct isl_chamber_list);
625 if (!c)
626 goto error;
627 c->c.vertices = isl_alloc_array(tab->mat->ctx, int, n_vertices);
628 if (!c->c.vertices)
629 goto error;
630 c->c.dom = isl_basic_set_from_basic_map(isl_basic_map_copy(tab->bmap));
631 c->c.dom = isl_basic_set_set_rational(c->c.dom);
632 c->c.dom = isl_basic_set_cow(c->c.dom);
633 c->c.dom = isl_basic_set_update_from_tab(c->c.dom, tab);
634 c->c.dom = isl_basic_set_simplify(c->c.dom);
635 c->c.dom = isl_basic_set_finalize(c->c.dom);
636 if (!c->c.dom)
637 goto error;
639 c->c.n_vertices = n_vertices;
641 for (i = 0, j = 0; i < vertices->n_vertices; ++i)
642 if (selection[i]) {
643 c->c.vertices[j] = i;
644 j++;
647 c->next = *list;
648 *list = c;
650 for (i = 0; i < n_frozen; ++i)
651 tab->con[i].frozen = 1;
653 if (isl_tab_rollback(tab, snap) < 0)
654 return -1;
656 return 0;
657 error:
658 free_chamber_list(c);
659 return -1;
662 struct isl_facet_todo {
663 struct isl_tab *tab; /* A tableau representation of the facet */
664 isl_basic_set *bset; /* A normalized basic set representation */
665 isl_vec *constraint; /* Constraint pointing to the other side */
666 struct isl_facet_todo *next;
669 static void free_todo(struct isl_facet_todo *todo)
671 while (todo) {
672 struct isl_facet_todo *next = todo->next;
674 isl_tab_free(todo->tab);
675 isl_basic_set_free(todo->bset);
676 isl_vec_free(todo->constraint);
677 free(todo);
679 todo = next;
683 static struct isl_facet_todo *create_todo(struct isl_tab *tab, int con)
685 int i;
686 int n_frozen;
687 struct isl_tab_undo *snap;
688 struct isl_facet_todo *todo;
690 snap = isl_tab_snap(tab);
692 for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
693 tab->con[i].frozen = 0;
694 n_frozen = i;
696 if (isl_tab_detect_redundant(tab) < 0)
697 return NULL;
699 todo = isl_calloc_type(tab->mat->ctx, struct isl_facet_todo);
700 if (!todo)
701 return NULL;
703 todo->constraint = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
704 if (!todo->constraint)
705 goto error;
706 isl_seq_neg(todo->constraint->el, tab->bmap->ineq[con], 1 + tab->n_var);
707 todo->bset = isl_basic_set_from_basic_map(isl_basic_map_copy(tab->bmap));
708 todo->bset = isl_basic_set_set_rational(todo->bset);
709 todo->bset = isl_basic_set_cow(todo->bset);
710 todo->bset = isl_basic_set_update_from_tab(todo->bset, tab);
711 todo->bset = isl_basic_set_simplify(todo->bset);
712 todo->bset = isl_basic_set_sort_constraints(todo->bset);
713 if (!todo->bset)
714 goto error;
715 ISL_F_SET(todo->bset, ISL_BASIC_SET_NORMALIZED);
716 todo->tab = isl_tab_dup(tab);
717 if (!todo->tab)
718 goto error;
720 for (i = 0; i < n_frozen; ++i)
721 tab->con[i].frozen = 1;
723 if (isl_tab_rollback(tab, snap) < 0)
724 goto error;
726 return todo;
727 error:
728 free_todo(todo);
729 return NULL;
732 /* Create todo items for all interior facets of the chamber represented
733 * by "tab" and collect them in "next".
735 static int init_todo(struct isl_facet_todo **next, struct isl_tab *tab)
737 int i;
738 struct isl_tab_undo *snap;
739 struct isl_facet_todo *todo;
741 snap = isl_tab_snap(tab);
743 for (i = 0; i < tab->n_con; ++i) {
744 if (tab->con[i].frozen)
745 continue;
746 if (tab->con[i].is_redundant)
747 continue;
749 if (isl_tab_select_facet(tab, i) < 0)
750 return -1;
752 todo = create_todo(tab, i);
753 if (!todo)
754 return -1;
756 todo->next = *next;
757 *next = todo;
759 if (isl_tab_rollback(tab, snap) < 0)
760 return -1;
763 return 0;
766 /* Does the linked list contain a todo item that is the opposite of "todo".
767 * If so, return 1 and remove the opposite todo item.
769 static int has_opposite(struct isl_facet_todo *todo,
770 struct isl_facet_todo **list)
772 for (; *list; list = &(*list)->next) {
773 int eq;
774 eq = isl_basic_set_plain_is_equal(todo->bset, (*list)->bset);
775 if (eq < 0)
776 return -1;
777 if (!eq)
778 continue;
779 todo = *list;
780 *list = todo->next;
781 todo->next = NULL;
782 free_todo(todo);
783 return 1;
786 return 0;
789 /* Create todo items for all interior facets of the chamber represented
790 * by "tab" and collect them in first->next, taking care to cancel
791 * opposite todo items.
793 static int update_todo(struct isl_facet_todo *first, struct isl_tab *tab)
795 int i;
796 struct isl_tab_undo *snap;
797 struct isl_facet_todo *todo;
799 snap = isl_tab_snap(tab);
801 for (i = 0; i < tab->n_con; ++i) {
802 int drop;
804 if (tab->con[i].frozen)
805 continue;
806 if (tab->con[i].is_redundant)
807 continue;
809 if (isl_tab_select_facet(tab, i) < 0)
810 return -1;
812 todo = create_todo(tab, i);
813 if (!todo)
814 return -1;
816 drop = has_opposite(todo, &first->next);
817 if (drop < 0)
818 return -1;
820 if (drop)
821 free_todo(todo);
822 else {
823 todo->next = first->next;
824 first->next = todo;
827 if (isl_tab_rollback(tab, snap) < 0)
828 return -1;
831 return 0;
834 /* Compute the chamber decomposition of the parametric polytope respresented
835 * by "bset" given the parametric vertices and their activity domains.
837 * We are only interested in full-dimensional chambers.
838 * Each of these chambers is the intersection of the activity domains of
839 * one or more vertices and the union of all chambers is equal to the
840 * projection of the entire parametric polytope onto the parameter space.
842 * We first create an initial chamber by intersecting as many activity
843 * domains as possible without ending up with an empty or lower-dimensional
844 * set. As a minor optimization, we only consider those activity domains
845 * that contain some arbitrary point.
847 * For each of interior facets of the chamber, we construct a todo item,
848 * containing the facet and a constraint containing the other side of the facet,
849 * for constructing the chamber on the other side.
850 * While their are any todo items left, we pick a todo item and
851 * create the required chamber by intersecting all activity domains
852 * that contain the facet and have a full-dimensional intersection with
853 * the other side of the facet. For each of the interior facets, we
854 * again create todo items, taking care to cancel opposite todo items.
856 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
857 __isl_take isl_vertices *vertices)
859 int i;
860 isl_ctx *ctx;
861 isl_vec *sample = NULL;
862 struct isl_tab *tab = NULL;
863 struct isl_tab_undo *snap;
864 unsigned nvar;
865 int *selection = NULL;
866 int n_chambers = 0;
867 struct isl_chamber_list *list = NULL;
868 struct isl_facet_todo *todo = NULL;
870 if (!bset || !vertices)
871 goto error;
873 ctx = isl_vertices_get_ctx(vertices);
874 selection = isl_alloc_array(ctx, int, vertices->n_vertices);
875 if (!selection)
876 goto error;
878 nvar = isl_basic_set_dim(bset, isl_dim_set);
879 bset = isl_basic_set_project_out(bset, isl_dim_set, 0, nvar);
881 tab = isl_tab_from_basic_set(bset);
882 for (i = 0; i < bset->n_ineq; ++i)
883 if (isl_tab_freeze_constraint(tab, i) < 0)
884 goto error;
885 if (isl_tab_track_bset(tab, bset) < 0)
886 goto error;
888 snap = isl_tab_snap(tab);
890 sample = isl_tab_get_sample_value(tab);
892 for (i = 0; i < vertices->n_vertices; ++i) {
893 selection[i] = isl_basic_set_contains(vertices->v[i].dom, sample);
894 if (selection[i] < 0)
895 goto error;
896 if (!selection[i])
897 continue;
898 selection[i] = can_intersect(tab, vertices->v[i].dom);
899 if (selection[i] < 0)
900 goto error;
903 if (isl_tab_detect_redundant(tab) < 0)
904 goto error;
906 if (add_chamber(&list, vertices, tab, selection) < 0)
907 goto error;
908 n_chambers++;
910 if (init_todo(&todo, tab) < 0)
911 goto error;
913 while (todo) {
914 struct isl_facet_todo *next;
916 if (isl_tab_rollback(tab, snap) < 0)
917 goto error;
919 if (isl_tab_add_ineq(tab, todo->constraint->el) < 0)
920 goto error;
921 if (isl_tab_freeze_constraint(tab, tab->n_con - 1) < 0)
922 goto error;
924 for (i = 0; i < vertices->n_vertices; ++i) {
925 selection[i] = bset_covers_tab(vertices->v[i].dom,
926 todo->tab);
927 if (selection[i] < 0)
928 goto error;
929 if (!selection[i])
930 continue;
931 selection[i] = can_intersect(tab, vertices->v[i].dom);
932 if (selection[i] < 0)
933 goto error;
936 if (isl_tab_detect_redundant(tab) < 0)
937 goto error;
939 if (add_chamber(&list, vertices, tab, selection) < 0)
940 goto error;
941 n_chambers++;
943 if (update_todo(todo, tab) < 0)
944 goto error;
946 next = todo->next;
947 todo->next = NULL;
948 free_todo(todo);
949 todo = next;
952 isl_vec_free(sample);
954 isl_tab_free(tab);
955 free(selection);
957 vertices = vertices_add_chambers(vertices, n_chambers, list);
959 for (i = 0; vertices && i < vertices->n_vertices; ++i) {
960 isl_basic_set_free(vertices->v[i].dom);
961 vertices->v[i].dom = NULL;
964 return vertices;
965 error:
966 free_chamber_list(list);
967 free_todo(todo);
968 isl_vec_free(sample);
969 isl_tab_free(tab);
970 free(selection);
971 if (!tab)
972 isl_basic_set_free(bset);
973 isl_vertices_free(vertices);
974 return NULL;
977 isl_ctx *isl_vertex_get_ctx(__isl_keep isl_vertex *vertex)
979 return vertex ? isl_vertices_get_ctx(vertex->vertices) : NULL;
982 int isl_vertex_get_id(__isl_keep isl_vertex *vertex)
984 return vertex ? vertex->id : -1;
987 __isl_give isl_basic_set *isl_vertex_get_domain(__isl_keep isl_vertex *vertex)
989 struct isl_vertex *v;
991 if (!vertex)
992 return NULL;
994 v = &vertex->vertices->v[vertex->id];
995 if (!v->dom) {
996 unsigned nvar;
997 nvar = isl_basic_set_dim(v->vertex, isl_dim_set);
998 v->dom = isl_basic_set_copy(v->vertex);
999 v->dom = isl_basic_set_project_out(v->dom, isl_dim_set, 0, nvar);
1002 return isl_basic_set_copy(v->dom);
1005 __isl_give isl_basic_set *isl_vertex_get_expr(__isl_keep isl_vertex *vertex)
1007 struct isl_vertex *v;
1009 if (!vertex)
1010 return NULL;
1012 v = &vertex->vertices->v[vertex->id];
1014 return isl_basic_set_copy(v->vertex);
1017 static __isl_give isl_vertex *isl_vertex_alloc(__isl_take isl_vertices *vertices,
1018 int id)
1020 isl_ctx *ctx;
1021 isl_vertex *vertex;
1023 if (!vertices)
1024 return NULL;
1026 ctx = isl_vertices_get_ctx(vertices);
1027 vertex = isl_alloc_type(ctx, isl_vertex);
1028 if (!vertex)
1029 goto error;
1031 vertex->vertices = vertices;
1032 vertex->id = id;
1034 return vertex;
1035 error:
1036 isl_vertices_free(vertices);
1037 return NULL;
1040 void isl_vertex_free(__isl_take isl_vertex *vertex)
1042 if (!vertex)
1043 return;
1044 isl_vertices_free(vertex->vertices);
1045 free(vertex);
1048 __isl_give isl_basic_set *isl_basic_set_set_integral(__isl_take isl_basic_set *bset)
1050 if (!bset)
1051 return NULL;
1053 if (!ISL_F_ISSET(bset, ISL_BASIC_MAP_RATIONAL))
1054 return bset;
1056 bset = isl_basic_set_cow(bset);
1057 if (!bset)
1058 return NULL;
1060 ISL_F_CLR(bset, ISL_BASIC_MAP_RATIONAL);
1062 return isl_basic_set_finalize(bset);
1065 isl_ctx *isl_cell_get_ctx(__isl_keep isl_cell *cell)
1067 return cell ? cell->dom->ctx : NULL;
1070 __isl_give isl_basic_set *isl_cell_get_domain(__isl_keep isl_cell *cell)
1072 return cell ? isl_basic_set_copy(cell->dom) : NULL;
1075 static __isl_give isl_cell *isl_cell_alloc(__isl_take isl_vertices *vertices,
1076 __isl_take isl_basic_set *dom, int id)
1078 int i;
1079 isl_cell *cell = NULL;
1081 if (!vertices || !dom)
1082 goto error;
1084 cell = isl_calloc_type(dom->ctx, isl_cell);
1085 if (!cell)
1086 goto error;
1088 cell->n_vertices = vertices->c[id].n_vertices;
1089 cell->ids = isl_alloc_array(dom->ctx, int, cell->n_vertices);
1090 if (!cell->ids)
1091 goto error;
1092 for (i = 0; i < cell->n_vertices; ++i)
1093 cell->ids[i] = vertices->c[id].vertices[i];
1094 cell->vertices = vertices;
1095 cell->dom = dom;
1097 return cell;
1098 error:
1099 isl_cell_free(cell);
1100 isl_vertices_free(vertices);
1101 isl_basic_set_free(dom);
1102 return NULL;
1105 void isl_cell_free(__isl_take isl_cell *cell)
1107 if (!cell)
1108 return;
1110 isl_vertices_free(cell->vertices);
1111 free(cell->ids);
1112 isl_basic_set_free(cell->dom);
1113 free(cell);
1116 /* Create a tableau of the cone obtained by first homogenizing the given
1117 * polytope and then making all inequalities strict by setting the
1118 * constant term to -1.
1120 static struct isl_tab *tab_for_shifted_cone(__isl_keep isl_basic_set *bset)
1122 int i;
1123 isl_vec *c = NULL;
1124 struct isl_tab *tab;
1126 if (!bset)
1127 return NULL;
1128 tab = isl_tab_alloc(bset->ctx, bset->n_ineq + 1,
1129 1 + isl_basic_set_total_dim(bset), 0);
1130 if (!tab)
1131 return NULL;
1132 tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1133 if (ISL_F_ISSET(bset, ISL_BASIC_MAP_EMPTY)) {
1134 if (isl_tab_mark_empty(tab) < 0)
1135 goto error;
1136 return tab;
1139 c = isl_vec_alloc(bset->ctx, 1 + 1 + isl_basic_set_total_dim(bset));
1140 if (!c)
1141 goto error;
1143 isl_int_set_si(c->el[0], 0);
1144 for (i = 0; i < bset->n_eq; ++i) {
1145 isl_seq_cpy(c->el + 1, bset->eq[i], c->size - 1);
1146 if (isl_tab_add_eq(tab, c->el) < 0)
1147 goto error;
1150 isl_int_set_si(c->el[0], -1);
1151 for (i = 0; i < bset->n_ineq; ++i) {
1152 isl_seq_cpy(c->el + 1, bset->ineq[i], c->size - 1);
1153 if (isl_tab_add_ineq(tab, c->el) < 0)
1154 goto error;
1155 if (tab->empty) {
1156 isl_vec_free(c);
1157 return tab;
1161 isl_seq_clr(c->el + 1, c->size - 1);
1162 isl_int_set_si(c->el[1], 1);
1163 if (isl_tab_add_ineq(tab, c->el) < 0)
1164 goto error;
1166 isl_vec_free(c);
1167 return tab;
1168 error:
1169 isl_vec_free(c);
1170 isl_tab_free(tab);
1171 return NULL;
1174 /* Compute an interior point of "bset" by selecting an interior
1175 * point in homogeneous space and projecting the point back down.
1177 static __isl_give isl_vec *isl_basic_set_interior_point(
1178 __isl_keep isl_basic_set *bset)
1180 isl_vec *vec;
1181 struct isl_tab *tab;
1183 tab = tab_for_shifted_cone(bset);
1184 vec = isl_tab_get_sample_value(tab);
1185 isl_tab_free(tab);
1186 if (!vec)
1187 return NULL;
1189 isl_seq_cpy(vec->el, vec->el + 1, vec->size - 1);
1190 vec->size--;
1192 return vec;
1195 /* Call "fn" on all chambers of the parametric polytope with the shared
1196 * facets of neighboring chambers only appearing in one of the chambers.
1198 * We pick an interior point from one of the chambers and then make
1199 * all constraints that do not satisfy this point strict.
1201 int isl_vertices_foreach_disjoint_cell(__isl_keep isl_vertices *vertices,
1202 int (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1204 int i, j;
1205 isl_vec *vec;
1206 isl_int v;
1207 isl_cell *cell;
1209 if (!vertices)
1210 return -1;
1212 if (vertices->n_chambers == 0)
1213 return 0;
1215 if (vertices->n_chambers == 1) {
1216 isl_basic_set *dom = isl_basic_set_copy(vertices->c[0].dom);
1217 dom = isl_basic_set_set_integral(dom);
1218 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, 0);
1219 if (!cell)
1220 return -1;
1221 return fn(cell, user);
1224 vec = isl_basic_set_interior_point(vertices->c[0].dom);
1225 if (!vec)
1226 return -1;
1228 isl_int_init(v);
1230 for (i = 0; i < vertices->n_chambers; ++i) {
1231 int r;
1232 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1233 dom = isl_basic_set_cow(dom);
1234 if (!dom)
1235 goto error;
1236 for (j = 0; i && j < dom->n_ineq; ++j) {
1237 isl_seq_inner_product(vec->el, dom->ineq[j], vec->size,
1238 &v);
1239 if (!isl_int_is_neg(v))
1240 continue;
1241 isl_int_sub_ui(dom->ineq[j][0], dom->ineq[j][0], 1);
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_int_clear(v);
1253 isl_vec_free(vec);
1255 return 0;
1256 error:
1257 isl_int_clear(v);
1258 isl_vec_free(vec);
1259 return -1;
1262 int isl_vertices_foreach_cell(__isl_keep isl_vertices *vertices,
1263 int (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1265 int i;
1266 isl_cell *cell;
1268 if (!vertices)
1269 return -1;
1271 if (vertices->n_chambers == 0)
1272 return 0;
1274 for (i = 0; i < vertices->n_chambers; ++i) {
1275 int r;
1276 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1278 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1279 if (!cell)
1280 return -1;
1282 r = fn(cell, user);
1283 if (r < 0)
1284 return -1;
1287 return 0;
1290 int isl_vertices_foreach_vertex(__isl_keep isl_vertices *vertices,
1291 int (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1293 int i;
1294 isl_vertex *vertex;
1296 if (!vertices)
1297 return -1;
1299 if (vertices->n_vertices == 0)
1300 return 0;
1302 for (i = 0; i < vertices->n_vertices; ++i) {
1303 int r;
1305 vertex = isl_vertex_alloc(isl_vertices_copy(vertices), i);
1306 if (!vertex)
1307 return -1;
1309 r = fn(vertex, user);
1310 if (r < 0)
1311 return -1;
1314 return 0;
1317 int isl_cell_foreach_vertex(__isl_keep isl_cell *cell,
1318 int (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1320 int i;
1321 isl_vertex *vertex;
1323 if (!cell)
1324 return -1;
1326 if (cell->n_vertices == 0)
1327 return 0;
1329 for (i = 0; i < cell->n_vertices; ++i) {
1330 int r;
1332 vertex = isl_vertex_alloc(isl_vertices_copy(cell->vertices),
1333 cell->ids[i]);
1334 if (!vertex)
1335 return -1;
1337 r = fn(vertex, user);
1338 if (r < 0)
1339 return -1;
1342 return 0;
1345 isl_ctx *isl_vertices_get_ctx(__isl_keep isl_vertices *vertices)
1347 return vertices ? vertices->bset->ctx : NULL;
1350 int isl_vertices_get_n_vertices(__isl_keep isl_vertices *vertices)
1352 return vertices ? vertices->n_vertices : -1;
1355 __isl_give isl_vertices *isl_morph_vertices(__isl_take isl_morph *morph,
1356 __isl_take isl_vertices *vertices)
1358 int i;
1359 isl_morph *param_morph = NULL;
1361 if (!morph || !vertices)
1362 goto error;
1364 isl_assert(vertices->bset->ctx, vertices->ref == 1, goto error);
1366 param_morph = isl_morph_copy(morph);
1367 param_morph = isl_morph_remove_dom_dims(param_morph, isl_dim_set,
1368 0, isl_morph_dom_dim(morph, isl_dim_set));
1369 param_morph = isl_morph_remove_ran_dims(param_morph, isl_dim_set,
1370 0, isl_morph_ran_dim(morph, isl_dim_set));
1372 for (i = 0; i < vertices->n_vertices; ++i) {
1373 vertices->v[i].dom = isl_morph_basic_set(
1374 isl_morph_copy(param_morph), vertices->v[i].dom);
1375 vertices->v[i].vertex = isl_morph_basic_set(
1376 isl_morph_copy(morph), vertices->v[i].vertex);
1377 if (!vertices->v[i].vertex)
1378 goto error;
1381 for (i = 0; i < vertices->n_chambers; ++i) {
1382 vertices->c[i].dom = isl_morph_basic_set(
1383 isl_morph_copy(param_morph), vertices->c[i].dom);
1384 if (!vertices->c[i].dom)
1385 goto error;
1388 isl_morph_free(param_morph);
1389 isl_morph_free(morph);
1390 return vertices;
1391 error:
1392 isl_morph_free(param_morph);
1393 isl_morph_free(morph);
1394 isl_vertices_free(vertices);
1395 return NULL;
1398 /* Construct a simplex isl_cell spanned by the vertices with indices in
1399 * "simplex_ids" and "other_ids" and call "fn" on this isl_cell.
1401 static int call_on_simplex(__isl_keep isl_cell *cell,
1402 int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1403 int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1405 int i;
1406 isl_ctx *ctx;
1407 struct isl_cell *simplex;
1409 ctx = isl_cell_get_ctx(cell);
1411 simplex = isl_calloc_type(ctx, struct isl_cell);
1412 if (!simplex)
1413 return -1;
1414 simplex->vertices = isl_vertices_copy(cell->vertices);
1415 if (!simplex->vertices)
1416 goto error;
1417 simplex->dom = isl_basic_set_copy(cell->dom);
1418 if (!simplex->dom)
1419 goto error;
1420 simplex->n_vertices = n_simplex + n_other;
1421 simplex->ids = isl_alloc_array(ctx, int, simplex->n_vertices);
1422 if (!simplex->ids)
1423 goto error;
1425 for (i = 0; i < n_simplex; ++i)
1426 simplex->ids[i] = simplex_ids[i];
1427 for (i = 0; i < n_other; ++i)
1428 simplex->ids[n_simplex + i] = other_ids[i];
1430 return fn(simplex, user);
1431 error:
1432 isl_cell_free(simplex);
1433 return -1;
1436 /* Check whether the parametric vertex described by "vertex"
1437 * lies on the facet corresponding to constraint "facet" of "bset".
1438 * The isl_vec "v" is a temporary vector than can be used by this function.
1440 * We eliminate the variables from the facet constraint using the
1441 * equalities defining the vertex and check if the result is identical
1442 * to zero.
1444 * It would probably be better to keep track of the constraints defining
1445 * a vertex during the vertex construction so that we could simply look
1446 * it up here.
1448 static int vertex_on_facet(__isl_keep isl_basic_set *vertex,
1449 __isl_keep isl_basic_set *bset, int facet, __isl_keep isl_vec *v)
1451 int i;
1452 isl_int m;
1454 isl_seq_cpy(v->el, bset->ineq[facet], v->size);
1456 isl_int_init(m);
1457 for (i = 0; i < vertex->n_eq; ++i) {
1458 int k = isl_seq_last_non_zero(vertex->eq[i], v->size);
1459 isl_seq_elim(v->el, vertex->eq[i], k, v->size, &m);
1461 isl_int_clear(m);
1463 return isl_seq_first_non_zero(v->el, v->size) == -1;
1466 /* Triangulate the polytope spanned by the vertices with ids
1467 * in "simplex_ids" and "other_ids" and call "fn" on each of
1468 * the resulting simplices.
1469 * If the input polytope is already a simplex, we simply call "fn".
1470 * Otherwise, we pick a point from "other_ids" and add it to "simplex_ids".
1471 * Then we consider each facet of "bset" that does not contain the point
1472 * we just picked, but does contain some of the other points in "other_ids"
1473 * and call ourselves recursively on the polytope spanned by the new
1474 * "simplex_ids" and those points in "other_ids" that lie on the facet.
1476 static int triangulate(__isl_keep isl_cell *cell, __isl_keep isl_vec *v,
1477 int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1478 int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1480 int i, j, k;
1481 int d, nparam;
1482 int *ids;
1483 isl_ctx *ctx;
1484 isl_basic_set *vertex;
1485 isl_basic_set *bset;
1487 ctx = isl_cell_get_ctx(cell);
1488 d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1489 nparam = isl_basic_set_dim(cell->vertices->bset, isl_dim_param);
1491 if (n_simplex + n_other == d + 1)
1492 return call_on_simplex(cell, simplex_ids, n_simplex,
1493 other_ids, n_other, fn, user);
1495 simplex_ids[n_simplex] = other_ids[0];
1496 vertex = cell->vertices->v[other_ids[0]].vertex;
1497 bset = cell->vertices->bset;
1499 ids = isl_alloc_array(ctx, int, n_other - 1);
1500 for (i = 0; i < bset->n_ineq; ++i) {
1501 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + nparam, d) == -1)
1502 continue;
1503 if (vertex_on_facet(vertex, bset, i, v))
1504 continue;
1506 for (j = 1, k = 0; j < n_other; ++j) {
1507 isl_basic_set *ov;
1508 ov = cell->vertices->v[other_ids[j]].vertex;
1509 if (vertex_on_facet(ov, bset, i, v))
1510 ids[k++] = other_ids[j];
1512 if (k == 0)
1513 continue;
1515 if (triangulate(cell, v, simplex_ids, n_simplex + 1,
1516 ids, k, fn, user) < 0)
1517 goto error;
1519 free(ids);
1521 return 0;
1522 error:
1523 free(ids);
1524 return -1;
1527 /* Triangulate the given cell and call "fn" on each of the resulting
1528 * simplices.
1530 int isl_cell_foreach_simplex(__isl_take isl_cell *cell,
1531 int (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1533 int d, total;
1534 int r;
1535 isl_ctx *ctx;
1536 isl_vec *v = NULL;
1537 int *simplex_ids = NULL;
1539 if (!cell)
1540 return -1;
1542 d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1543 total = isl_basic_set_total_dim(cell->vertices->bset);
1545 if (cell->n_vertices == d + 1)
1546 return fn(cell, user);
1548 ctx = isl_cell_get_ctx(cell);
1549 simplex_ids = isl_alloc_array(ctx, int, d + 1);
1550 if (!simplex_ids)
1551 goto error;
1553 v = isl_vec_alloc(ctx, 1 + total);
1554 if (!v)
1555 goto error;
1557 r = triangulate(cell, v, simplex_ids, 0,
1558 cell->ids, cell->n_vertices, fn, user);
1560 isl_vec_free(v);
1561 free(simplex_ids);
1563 isl_cell_free(cell);
1565 return r;
1566 error:
1567 free(simplex_ids);
1568 isl_vec_free(v);
1569 isl_cell_free(cell);
1570 return -1;