polly/lib/External/isl/isl_vertices.c

   1 /*
   2  * Copyright 2010      INRIA Saclay
   3  *
   4  * Use of this software is governed by the MIT license
   5  *
   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  */
  10
  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>
  21
  22 #define SELECTED        1
  23 #define DESELECTED      -1
  24 #define UNSELECTED      0
  25
  26 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
  27         __isl_take isl_vertices *vertices);
  28
  29 __isl_give isl_vertices *isl_vertices_copy(__isl_keep isl_vertices *vertices)
  30 {
  31         if (!vertices)
  32                 return NULL;
  33
  34         vertices->ref++;
  35         return vertices;
  36 }
  37
  38 __isl_null isl_vertices *isl_vertices_free(__isl_take isl_vertices *vertices)
  39 {
  40         int i;
  41
  42         if (!vertices)
  43                 return NULL;
  44
  45         if (--vertices->ref > 0)
  46                 return NULL;
  47
  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);
  51         }
  52         free(vertices->v);
  53
  54         for (i = 0; i < vertices->n_chambers; ++i) {
  55                 free(vertices->c[i].vertices);
  56                 isl_basic_set_free(vertices->c[i].dom);
  57         }
  58         free(vertices->c);
  59
  60         isl_basic_set_free(vertices->bset);
  61         free(vertices);
  62
  63         return NULL;
  64 }
  65
  66 struct isl_vertex_list {
  67         struct isl_vertex v;
  68         struct isl_vertex_list *next;
  69 };
  70
  71 static struct isl_vertex_list *free_vertex_list(struct isl_vertex_list *list)
  72 {
  73         struct isl_vertex_list *next;
  74
  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);
  80         }
  81
  82         return NULL;
  83 }
  84
  85 static __isl_give isl_vertices *vertices_from_list(__isl_keep isl_basic_set *bset,
  86         int n_vertices, struct isl_vertex_list *list)
  87 {
  88         int i;
  89         struct isl_vertex_list *next;
  90         isl_vertices *vertices;
  91
  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;
 101
 102         for (i = 0; list; list = next, i++) {
 103                 next = list->next;
 104                 vertices->v[i] = list->v;
 105                 free(list);
 106         }
 107
 108         return vertices;
 109 error:
 110         isl_vertices_free(vertices);
 111         free_vertex_list(list);
 112         return NULL;
 113 }
 114
 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.
 121  */
 122 static isl_bool add_vertex(struct isl_vertex_list **list,
 123         __isl_keep isl_basic_set *bset, struct isl_tab *tab)
 124 {
 125         isl_size nvar;
 126         struct isl_vertex_list *v = NULL;
 127
 128         if (isl_tab_detect_implicit_equalities(tab) < 0)
 129                 return isl_bool_error;
 130
 131         nvar = isl_basic_set_dim(bset, isl_dim_set);
 132         if (nvar < 0)
 133                 return isl_bool_error;
 134
 135         v = isl_calloc_type(tab->mat->ctx, struct isl_vertex_list);
 136         if (!v)
 137                 goto error;
 138
 139         v->v.vertex = isl_basic_set_copy(bset);
 140         v->v.vertex = isl_basic_set_cow(v->v.vertex);
 141         v->v.vertex = isl_basic_set_update_from_tab(v->v.vertex, tab);
 142         v->v.vertex = isl_basic_set_simplify(v->v.vertex);
 143         v->v.vertex = isl_basic_set_finalize(v->v.vertex);
 144         if (!v->v.vertex)
 145                 goto error;
 146         isl_assert(bset->ctx, v->v.vertex->n_eq >= nvar, goto error);
 147         v->v.dom = isl_basic_set_copy(v->v.vertex);
 148         v->v.dom = isl_basic_set_params(v->v.dom);
 149         if (!v->v.dom)
 150                 goto error;
 151
 152         if (v->v.dom->n_eq > 0) {
 153                 free_vertex_list(v);
 154                 return isl_bool_false;
 155         }
 156
 157         v->next = *list;
 158         *list = v;
 159
 160         return isl_bool_true;
 161 error:
 162         free_vertex_list(v);
 163         return isl_bool_error;
 164 }
 165
 166 /* Compute the parametric vertices and the chamber decomposition
 167  * of an empty parametric polytope.
 168  */
 169 static __isl_give isl_vertices *vertices_empty(__isl_keep isl_basic_set *bset)
 170 {
 171         isl_vertices *vertices;
 172
 173         if (!bset)
 174                 return NULL;
 175
 176         vertices = isl_calloc_type(bset->ctx, isl_vertices);
 177         if (!vertices)
 178                 return NULL;
 179         vertices->bset = isl_basic_set_copy(bset);
 180         vertices->ref = 1;
 181
 182         vertices->n_vertices = 0;
 183         vertices->n_chambers = 0;
 184
 185         return vertices;
 186 }
 187
 188 /* Compute the parametric vertices and the chamber decomposition
 189  * of the parametric polytope defined using the same constraints
 190  * as "bset" in the 0D case.
 191  * There is exactly one 0D vertex and a single chamber containing
 192  * the vertex.
 193  */
 194 static __isl_give isl_vertices *vertices_0D(__isl_keep isl_basic_set *bset)
 195 {
 196         isl_vertices *vertices;
 197
 198         if (!bset)
 199                 return NULL;
 200
 201         vertices = isl_calloc_type(bset->ctx, isl_vertices);
 202         if (!vertices)
 203                 return NULL;
 204         vertices->ref = 1;
 205         vertices->bset = isl_basic_set_copy(bset);
 206
 207         vertices->v = isl_calloc_array(bset->ctx, struct isl_vertex, 1);
 208         if (!vertices->v)
 209                 goto error;
 210         vertices->n_vertices = 1;
 211         vertices->v[0].vertex = isl_basic_set_copy(bset);
 212         vertices->v[0].dom = isl_basic_set_params(isl_basic_set_copy(bset));
 213         if (!vertices->v[0].vertex || !vertices->v[0].dom)
 214                 goto error;
 215
 216         vertices->c = isl_calloc_array(bset->ctx, struct isl_chamber, 1);
 217         if (!vertices->c)
 218                 goto error;
 219         vertices->n_chambers = 1;
 220         vertices->c[0].n_vertices = 1;
 221         vertices->c[0].vertices = isl_calloc_array(bset->ctx, int, 1);
 222         if (!vertices->c[0].vertices)
 223                 goto error;
 224         vertices->c[0].dom = isl_basic_set_copy(vertices->v[0].dom);
 225         if (!vertices->c[0].dom)
 226                 goto error;
 227
 228         return vertices;
 229 error:
 230         isl_vertices_free(vertices);
 231         return NULL;
 232 }
 233
 234 /* Is the row pointed to by "f" linearly independent of the "n" first
 235  * rows in "facets"?
 236  */
 237 static isl_bool is_independent(__isl_keep isl_mat *facets, int n, isl_int *f)
 238 {
 239         isl_size rank;
 240
 241         if (isl_seq_first_non_zero(f, facets->n_col) < 0)
 242                 return isl_bool_false;
 243
 244         isl_seq_cpy(facets->row[n], f, facets->n_col);
 245         facets->n_row = n + 1;
 246         rank = isl_mat_rank(facets);
 247         if (rank < 0)
 248                 return isl_bool_error;
 249
 250         return isl_bool_ok(rank == n + 1);
 251 }
 252
 253 /* Check whether we can select constraint "level", given the current selection
 254  * reflected by facets in "tab", the rows of "facets" and the earlier
 255  * "selected" elements of "selection".
 256  *
 257  * If the constraint is (strictly) redundant in the tableau, selecting it would
 258  * result in an empty tableau, so it can't be selected.
 259  * If the set variable part of the constraint is not linearly independent
 260  * of the set variable parts of the already selected constraints,
 261  * the constraint cannot be selected.
 262  * If selecting the constraint results in an empty tableau, the constraint
 263  * cannot be selected.
 264  * Finally, if selecting the constraint results in some explicitly
 265  * deselected constraints turning into equalities, then the corresponding
 266  * vertices have already been generated, so the constraint cannot be selected.
 267  */
 268 static isl_bool can_select(__isl_keep isl_basic_set *bset, int level,
 269         struct isl_tab *tab, __isl_keep isl_mat *facets, int selected,
 270         int *selection)
 271 {
 272         int i;
 273         isl_bool indep;
 274         unsigned ovar;
 275         struct isl_tab_undo *snap;
 276
 277         if (isl_tab_is_redundant(tab, level))
 278                 return isl_bool_false;
 279
 280         ovar = isl_space_offset(bset->dim, isl_dim_set);
 281
 282         indep = is_independent(facets, selected, bset->ineq[level] + 1 + ovar);
 283         if (indep < 0 || !indep)
 284                 return indep;
 285
 286         snap = isl_tab_snap(tab);
 287         if (isl_tab_select_facet(tab, level) < 0)
 288                 return isl_bool_error;
 289
 290         if (tab->empty) {
 291                 if (isl_tab_rollback(tab, snap) < 0)
 292                         return isl_bool_error;
 293                 return isl_bool_false;
 294         }
 295
 296         for (i = 0; i < level; ++i) {
 297                 int sgn;
 298
 299                 if (selection[i] != DESELECTED)
 300                         continue;
 301
 302                 if (isl_tab_is_equality(tab, i))
 303                         sgn = 0;
 304                 else if (isl_tab_is_redundant(tab, i))
 305                         sgn = 1;
 306                 else
 307                         sgn = isl_tab_sign_of_max(tab, i);
 308                 if (sgn < -1)
 309                         return isl_bool_error;
 310                 if (sgn <= 0) {
 311                         if (isl_tab_rollback(tab, snap) < 0)
 312                                 return isl_bool_error;
 313                         return isl_bool_false;
 314                 }
 315         }
 316
 317         return isl_bool_true;
 318 }
 319
 320 /* Compute the parametric vertices and the chamber decomposition
 321  * of a parametric polytope that is not full-dimensional.
 322  *
 323  * Simply map the parametric polytope to a lower dimensional space
 324  * and map the resulting vertices back.
 325  */
 326 static __isl_give isl_vertices *lower_dim_vertices(
 327         __isl_take isl_basic_set *bset)
 328 {
 329         isl_morph *morph;
 330         isl_vertices *vertices;
 331
 332         morph = isl_basic_set_full_compression(bset);
 333         bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
 334
 335         vertices = isl_basic_set_compute_vertices(bset);
 336         isl_basic_set_free(bset);
 337
 338         morph = isl_morph_inverse(morph);
 339
 340         vertices = isl_morph_vertices(morph, vertices);
 341
 342         return vertices;
 343 }
 344
 345 /* Compute the parametric vertices and the chamber decomposition
 346  * of a parametric polytope "bset" that is not full-dimensional.
 347  * Additionally, free both "copy" and "tab".
 348  */
 349 static __isl_give isl_vertices *lower_dim_vertices_free(
 350         __isl_take isl_basic_set *bset, __isl_take isl_basic_set *copy,
 351         struct isl_tab *tab)
 352 {
 353         isl_basic_set_free(copy);
 354         isl_tab_free(tab);
 355         return lower_dim_vertices(bset);
 356 }
 357
 358 /* Detect implicit equality constraints in "bset" using the tableau
 359  * representation "tab".
 360  * Return a copy of "bset" with the implicit equality constraints
 361  * made explicit, leaving the original "bset" unmodified.
 362  */
 363 static __isl_give isl_basic_set *detect_implicit_equality_constraints(
 364         __isl_keep isl_basic_set *bset, struct isl_tab *tab)
 365 {
 366         if (isl_tab_detect_implicit_equalities(tab) < 0)
 367                 return NULL;
 368
 369         bset = isl_basic_set_copy(bset);
 370         bset = isl_basic_set_cow(bset);
 371         bset = isl_basic_set_update_from_tab(bset, tab);
 372
 373         return bset;
 374 }
 375
 376 /* Compute the parametric vertices and the chamber decomposition
 377  * of the parametric polytope defined using the same constraints
 378  * as "bset".  "bset" is assumed to have no existentially quantified
 379  * variables.
 380  *
 381  * The vertices themselves are computed in a fairly simplistic way.
 382  * We simply run through all combinations of d constraints,
 383  * with d the number of set variables, and check if those d constraints
 384  * define a vertex.  To avoid the generation of duplicate vertices,
 385  * which may happen if a vertex is defined by more than d constraints,
 386  * we make sure we only generate the vertex for the d constraints with
 387  * smallest index.
 388  *
 389  * Only potential vertices with a full-dimensional activity domain
 390  * are considered.  However, if the input has (implicit) equality
 391  * constraints among the parameters, then activity domain
 392  * should be considered full-dimensional if it does not satisfy
 393  * any extra equality constraints beyond those of the input.
 394  * The implicit equality constraints of the input are therefore first detected.
 395  * If there are any, then the input is mapped to a lower dimensional space
 396  * such that the check for full-dimensional activity domains
 397  * can be performed with respect to a full-dimensional space.
 398  * Note that it is important to leave "bset" unmodified while detecting
 399  * equality constraints since the inequality constraints of "bset"
 400  * are assumed to correspond to those of the tableau.
 401  *
 402  * We set up a tableau and keep track of which facets have been
 403  * selected.  The tableau is marked strict_redundant so that we can be
 404  * sure that any constraint that is marked redundant (and that is not
 405  * also marked zero) is not an equality.
 406  * If a constraint is marked DESELECTED, it means the constraint was
 407  * SELECTED before (in combination with the same selection of earlier
 408  * constraints).  If such a deselected constraint turns out to be an
 409  * equality, then any vertex that may still be found with the current
 410  * selection has already been generated when the constraint was selected.
 411  * A constraint is marked UNSELECTED when there is no way selecting
 412  * the constraint could lead to a vertex (in combination with the current
 413  * selection of earlier constraints).
 414  *
 415  * The set variable coefficients of the selected constraints are stored
 416  * in the facets matrix.
 417  */
 418 __isl_give isl_vertices *isl_basic_set_compute_vertices(
 419         __isl_keep isl_basic_set *bset)
 420 {
 421         struct isl_tab *tab;
 422         int level;
 423         int init;
 424         isl_size n_eq;
 425         isl_size nvar;
 426         int *selection = NULL;
 427         int selected;
 428         struct isl_tab_undo **snap = NULL;
 429         isl_mat *facets = NULL;
 430         struct isl_vertex_list *list = NULL;
 431         int n_vertices = 0;
 432         isl_vertices *vertices;
 433         isl_basic_set *copy;
 434         isl_basic_set *test;
 435
 436         if (!bset)
 437                 return NULL;
 438
 439         if (isl_basic_set_plain_is_empty(bset))
 440                 return vertices_empty(bset);
 441
 442         if (bset->n_eq != 0)
 443                 return lower_dim_vertices(isl_basic_set_copy(bset));
 444
 445         if (isl_basic_set_check_no_locals(bset) < 0)
 446                 return NULL;
 447
 448         nvar = isl_basic_set_dim(bset, isl_dim_set);
 449         if (nvar < 0)
 450                 return NULL;
 451         if (nvar == 0)
 452                 return vertices_0D(bset);
 453
 454         copy = isl_basic_set_copy(bset);
 455         copy = isl_basic_set_set_rational(copy);
 456         if (!copy)
 457                 return NULL;
 458
 459         tab = isl_tab_from_basic_set(copy, 0);
 460         if (!tab)
 461                 goto error;
 462         tab->strict_redundant = 1;
 463
 464         if (tab->empty) {
 465                 vertices = vertices_empty(copy);
 466                 isl_basic_set_free(copy);
 467                 isl_tab_free(tab);
 468                 return vertices;
 469         }
 470
 471         test = detect_implicit_equality_constraints(bset, tab);
 472         n_eq = isl_basic_set_n_equality(test);
 473         if (n_eq < 0)
 474                 test = isl_basic_set_free(test);
 475         if (n_eq < 0 || n_eq > 0)
 476                 return lower_dim_vertices_free(test, copy, tab);
 477         isl_basic_set_free(test);
 478
 479         selection = isl_alloc_array(copy->ctx, int, copy->n_ineq);
 480         snap = isl_alloc_array(copy->ctx, struct isl_tab_undo *, copy->n_ineq);
 481         facets = isl_mat_alloc(copy->ctx, nvar, nvar);
 482         if ((copy->n_ineq && (!selection || !snap)) || !facets)
 483                 goto error;
 484
 485         level = 0;
 486         init = 1;
 487         selected = 0;
 488
 489         while (level >= 0) {
 490                 if (level >= copy->n_ineq ||
 491                     (!init && selection[level] != SELECTED)) {
 492                         --level;
 493                         init = 0;
 494                         continue;
 495                 }
 496                 if (init) {
 497                         isl_bool ok;
 498                         snap[level] = isl_tab_snap(tab);
 499                         ok = can_select(copy, level, tab, facets, selected,
 500                                         selection);
 501                         if (ok < 0)
 502                                 goto error;
 503                         if (ok) {
 504                                 selection[level] = SELECTED;
 505                                 selected++;
 506                         } else
 507                                 selection[level] = UNSELECTED;
 508                 } else {
 509                         selection[level] = DESELECTED;
 510                         selected--;
 511                         if (isl_tab_rollback(tab, snap[level]) < 0)
 512                                 goto error;
 513                 }
 514                 if (selected == nvar) {
 515                         if (tab->n_dead == nvar) {
 516                                 isl_bool added = add_vertex(&list, copy, tab);
 517                                 if (added < 0)
 518                                         goto error;
 519                                 if (added)
 520                                         n_vertices++;
 521                         }
 522                         init = 0;
 523                         continue;
 524                 }
 525                 ++level;
 526                 init = 1;
 527         }
 528
 529         isl_mat_free(facets);
 530         free(selection);
 531         free(snap);
 532
 533         isl_tab_free(tab);
 534
 535         vertices = vertices_from_list(copy, n_vertices, list);
 536
 537         vertices = compute_chambers(copy, vertices);
 538
 539         return vertices;
 540 error:
 541         free_vertex_list(list);
 542         isl_mat_free(facets);
 543         free(selection);
 544         free(snap);
 545         isl_tab_free(tab);
 546         isl_basic_set_free(copy);
 547         return NULL;
 548 }
 549
 550 struct isl_chamber_list {
 551         struct isl_chamber c;
 552         struct isl_chamber_list *next;
 553 };
 554
 555 static void free_chamber_list(struct isl_chamber_list *list)
 556 {
 557         struct isl_chamber_list *next;
 558
 559         for (; list; list = next) {
 560                 next = list->next;
 561                 isl_basic_set_free(list->c.dom);
 562                 free(list->c.vertices);
 563                 free(list);
 564         }
 565 }
 566
 567 /* Check whether the basic set "bset" is a superset of the basic set described
 568  * by "tab", i.e., check whether all constraints of "bset" are redundant.
 569  */
 570 static isl_bool bset_covers_tab(__isl_keep isl_basic_set *bset,
 571         struct isl_tab *tab)
 572 {
 573         int i;
 574
 575         if (!bset || !tab)
 576                 return isl_bool_error;
 577
 578         for (i = 0; i < bset->n_ineq; ++i) {
 579                 enum isl_ineq_type type = isl_tab_ineq_type(tab, bset->ineq[i]);
 580                 switch (type) {
 581                 case isl_ineq_error:            return isl_bool_error;
 582                 case isl_ineq_redundant:        continue;
 583                 default:                        return isl_bool_false;
 584                 }
 585         }
 586
 587         return isl_bool_true;
 588 }
 589
 590 static __isl_give isl_vertices *vertices_add_chambers(
 591         __isl_take isl_vertices *vertices, int n_chambers,
 592         struct isl_chamber_list *list)
 593 {
 594         int i;
 595         isl_ctx *ctx;
 596         struct isl_chamber_list *next;
 597
 598         ctx = isl_vertices_get_ctx(vertices);
 599         vertices->c = isl_alloc_array(ctx, struct isl_chamber, n_chambers);
 600         if (!vertices->c)
 601                 goto error;
 602         vertices->n_chambers = n_chambers;
 603
 604         for (i = 0; list; list = next, i++) {
 605                 next = list->next;
 606                 vertices->c[i] = list->c;
 607                 free(list);
 608         }
 609
 610         return vertices;
 611 error:
 612         isl_vertices_free(vertices);
 613         free_chamber_list(list);
 614         return NULL;
 615 }
 616
 617 /* Can "tab" be intersected with "bset" without resulting in
 618  * a lower-dimensional set.
 619  * "bset" itself is assumed to be full-dimensional.
 620  */
 621 static isl_bool can_intersect(struct isl_tab *tab,
 622         __isl_keep isl_basic_set *bset)
 623 {
 624         int i;
 625         struct isl_tab_undo *snap;
 626
 627         if (bset->n_eq > 0)
 628                 isl_die(isl_basic_set_get_ctx(bset), isl_error_internal,
 629                         "expecting full-dimensional input",
 630                         return isl_bool_error);
 631
 632         if (isl_tab_extend_cons(tab, bset->n_ineq) < 0)
 633                 return isl_bool_error;
 634
 635         snap = isl_tab_snap(tab);
 636
 637         for (i = 0; i < bset->n_ineq; ++i) {
 638                 enum isl_ineq_type type;
 639
 640                 type = isl_tab_ineq_type(tab, bset->ineq[i]);
 641                 if (type < 0)
 642                         return isl_bool_error;
 643                 if (type == isl_ineq_redundant)
 644                         continue;
 645                 if (isl_tab_add_ineq(tab, bset->ineq[i]) < 0)
 646                         return isl_bool_error;
 647         }
 648
 649         if (isl_tab_detect_implicit_equalities(tab) < 0)
 650                 return isl_bool_error;
 651         if (tab->n_dead) {
 652                 if (isl_tab_rollback(tab, snap) < 0)
 653                         return isl_bool_error;
 654                 return isl_bool_false;
 655         }
 656
 657         return isl_bool_true;
 658 }
 659
 660 static int add_chamber(struct isl_chamber_list **list,
 661         __isl_keep isl_vertices *vertices, struct isl_tab *tab, int *selection)
 662 {
 663         int n_frozen;
 664         int i, j;
 665         int n_vertices = 0;
 666         struct isl_tab_undo *snap;
 667         struct isl_chamber_list *c = NULL;
 668
 669         for (i = 0; i < vertices->n_vertices; ++i)
 670                 if (selection[i])
 671                         n_vertices++;
 672
 673         snap = isl_tab_snap(tab);
 674
 675         for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
 676                 tab->con[i].frozen = 0;
 677         n_frozen = i;
 678
 679         if (isl_tab_detect_redundant(tab) < 0)
 680                 return -1;
 681
 682         c = isl_calloc_type(tab->mat->ctx, struct isl_chamber_list);
 683         if (!c)
 684                 goto error;
 685         c->c.vertices = isl_alloc_array(tab->mat->ctx, int, n_vertices);
 686         if (n_vertices && !c->c.vertices)
 687                 goto error;
 688         c->c.dom = isl_basic_set_copy(isl_tab_peek_bset(tab));
 689         c->c.dom = isl_basic_set_set_rational(c->c.dom);
 690         c->c.dom = isl_basic_set_cow(c->c.dom);
 691         c->c.dom = isl_basic_set_update_from_tab(c->c.dom, tab);
 692         c->c.dom = isl_basic_set_simplify(c->c.dom);
 693         c->c.dom = isl_basic_set_finalize(c->c.dom);
 694         if (!c->c.dom)
 695                 goto error;
 696
 697         c->c.n_vertices = n_vertices;
 698
 699         for (i = 0, j = 0; i < vertices->n_vertices; ++i)
 700                 if (selection[i]) {
 701                         c->c.vertices[j] = i;
 702                         j++;
 703                 }
 704
 705         c->next = *list;
 706         *list = c;
 707
 708         for (i = 0; i < n_frozen; ++i)
 709                 tab->con[i].frozen = 1;
 710
 711         if (isl_tab_rollback(tab, snap) < 0)
 712                 return -1;
 713
 714         return 0;
 715 error:
 716         free_chamber_list(c);
 717         return -1;
 718 }
 719
 720 struct isl_facet_todo {
 721         struct isl_tab *tab;    /* A tableau representation of the facet */
 722         isl_basic_set *bset;    /* A normalized basic set representation */
 723         isl_vec *constraint;    /* Constraint pointing to the other side */
 724         struct isl_facet_todo *next;
 725 };
 726
 727 static void free_todo(struct isl_facet_todo *todo)
 728 {
 729         while (todo) {
 730                 struct isl_facet_todo *next = todo->next;
 731
 732                 isl_tab_free(todo->tab);
 733                 isl_basic_set_free(todo->bset);
 734                 isl_vec_free(todo->constraint);
 735                 free(todo);
 736
 737                 todo = next;
 738         }
 739 }
 740
 741 static struct isl_facet_todo *create_todo(struct isl_tab *tab, int con)
 742 {
 743         int i;
 744         int n_frozen;
 745         struct isl_tab_undo *snap;
 746         struct isl_facet_todo *todo;
 747
 748         snap = isl_tab_snap(tab);
 749
 750         for (i = 0; i < tab->n_con && tab->con[i].frozen; ++i)
 751                 tab->con[i].frozen = 0;
 752         n_frozen = i;
 753
 754         if (isl_tab_detect_redundant(tab) < 0)
 755                 return NULL;
 756
 757         todo = isl_calloc_type(tab->mat->ctx, struct isl_facet_todo);
 758         if (!todo)
 759                 return NULL;
 760
 761         todo->constraint = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
 762         if (!todo->constraint)
 763                 goto error;
 764         isl_seq_neg(todo->constraint->el, tab->bmap->ineq[con], 1 + tab->n_var);
 765         todo->bset = isl_basic_set_copy(isl_tab_peek_bset(tab));
 766         todo->bset = isl_basic_set_set_rational(todo->bset);
 767         todo->bset = isl_basic_set_cow(todo->bset);
 768         todo->bset = isl_basic_set_update_from_tab(todo->bset, tab);
 769         todo->bset = isl_basic_set_simplify(todo->bset);
 770         todo->bset = isl_basic_set_sort_constraints(todo->bset);
 771         if (!todo->bset)
 772                 goto error;
 773         ISL_F_SET(todo->bset, ISL_BASIC_SET_NO_REDUNDANT);
 774         todo->tab = isl_tab_dup(tab);
 775         if (!todo->tab)
 776                 goto error;
 777
 778         for (i = 0; i < n_frozen; ++i)
 779                 tab->con[i].frozen = 1;
 780
 781         if (isl_tab_rollback(tab, snap) < 0)
 782                 goto error;
 783
 784         return todo;
 785 error:
 786         free_todo(todo);
 787         return NULL;
 788 }
 789
 790 /* Create todo items for all interior facets of the chamber represented
 791  * by "tab" and collect them in "next".
 792  */
 793 static int init_todo(struct isl_facet_todo **next, struct isl_tab *tab)
 794 {
 795         int i;
 796         struct isl_tab_undo *snap;
 797         struct isl_facet_todo *todo;
 798
 799         snap = isl_tab_snap(tab);
 800
 801         for (i = 0; i < tab->n_con; ++i) {
 802                 if (tab->con[i].frozen)
 803                         continue;
 804                 if (tab->con[i].is_redundant)
 805                         continue;
 806
 807                 if (isl_tab_select_facet(tab, i) < 0)
 808                         return -1;
 809
 810                 todo = create_todo(tab, i);
 811                 if (!todo)
 812                         return -1;
 813
 814                 todo->next = *next;
 815                 *next = todo;
 816
 817                 if (isl_tab_rollback(tab, snap) < 0)
 818                         return -1;
 819         }
 820
 821         return 0;
 822 }
 823
 824 /* Does the linked list contain a todo item that is the opposite of "todo".
 825  * If so, return 1 and remove the opposite todo item.
 826  */
 827 static int has_opposite(struct isl_facet_todo *todo,
 828         struct isl_facet_todo **list)
 829 {
 830         for (; *list; list = &(*list)->next) {
 831                 int eq;
 832                 eq = isl_basic_set_plain_is_equal(todo->bset, (*list)->bset);
 833                 if (eq < 0)
 834                         return -1;
 835                 if (!eq)
 836                         continue;
 837                 todo = *list;
 838                 *list = todo->next;
 839                 todo->next = NULL;
 840                 free_todo(todo);
 841                 return 1;
 842         }
 843
 844         return 0;
 845 }
 846
 847 /* Create todo items for all interior facets of the chamber represented
 848  * by "tab" and collect them in first->next, taking care to cancel
 849  * opposite todo items.
 850  */
 851 static int update_todo(struct isl_facet_todo *first, struct isl_tab *tab)
 852 {
 853         int i;
 854         struct isl_tab_undo *snap;
 855         struct isl_facet_todo *todo;
 856
 857         snap = isl_tab_snap(tab);
 858
 859         for (i = 0; i < tab->n_con; ++i) {
 860                 int drop;
 861
 862                 if (tab->con[i].frozen)
 863                         continue;
 864                 if (tab->con[i].is_redundant)
 865                         continue;
 866
 867                 if (isl_tab_select_facet(tab, i) < 0)
 868                         return -1;
 869
 870                 todo = create_todo(tab, i);
 871                 if (!todo)
 872                         return -1;
 873
 874                 drop = has_opposite(todo, &first->next);
 875                 if (drop < 0)
 876                         return -1;
 877
 878                 if (drop)
 879                         free_todo(todo);
 880                 else {
 881                         todo->next = first->next;
 882                         first->next = todo;
 883                 }
 884
 885                 if (isl_tab_rollback(tab, snap) < 0)
 886                         return -1;
 887         }
 888
 889         return 0;
 890 }
 891
 892 /* Compute the chamber decomposition of the parametric polytope respresented
 893  * by "bset" given the parametric vertices and their activity domains.
 894  *
 895  * We are only interested in full-dimensional chambers.
 896  * Each of these chambers is the intersection of the activity domains of
 897  * one or more vertices and the union of all chambers is equal to the
 898  * projection of the entire parametric polytope onto the parameter space.
 899  *
 900  * We first create an initial chamber by intersecting as many activity
 901  * domains as possible without ending up with an empty or lower-dimensional
 902  * set.  As a minor optimization, we only consider those activity domains
 903  * that contain some arbitrary point.
 904  *
 905  * For each of the interior facets of the chamber, we construct a todo item,
 906  * containing the facet and a constraint containing the other side of the facet,
 907  * for constructing the chamber on the other side.
 908  * While their are any todo items left, we pick a todo item and
 909  * create the required chamber by intersecting all activity domains
 910  * that contain the facet and have a full-dimensional intersection with
 911  * the other side of the facet.  For each of the interior facets, we
 912  * again create todo items, taking care to cancel opposite todo items.
 913  */
 914 static __isl_give isl_vertices *compute_chambers(__isl_take isl_basic_set *bset,
 915         __isl_take isl_vertices *vertices)
 916 {
 917         int i;
 918         isl_ctx *ctx;
 919         isl_size n_eq;
 920         isl_vec *sample = NULL;
 921         struct isl_tab *tab = NULL;
 922         struct isl_tab_undo *snap;
 923         int *selection = NULL;
 924         int n_chambers = 0;
 925         struct isl_chamber_list *list = NULL;
 926         struct isl_facet_todo *todo = NULL;
 927
 928         if (!bset || !vertices)
 929                 goto error;
 930
 931         ctx = isl_vertices_get_ctx(vertices);
 932         selection = isl_alloc_array(ctx, int, vertices->n_vertices);
 933         if (vertices->n_vertices && !selection)
 934                 goto error;
 935
 936         bset = isl_basic_set_params(bset);
 937         n_eq = isl_basic_set_n_equality(bset);
 938         if (n_eq < 0)
 939                 goto error;
 940         if (n_eq > 0)
 941                 isl_die(isl_basic_set_get_ctx(bset), isl_error_internal,
 942                         "expecting full-dimensional input", goto error);
 943
 944         tab = isl_tab_from_basic_set(bset, 1);
 945         if (!tab)
 946                 goto error;
 947         for (i = 0; i < bset->n_ineq; ++i)
 948                 if (isl_tab_freeze_constraint(tab, i) < 0)
 949                         goto error;
 950         isl_basic_set_free(bset);
 951
 952         snap = isl_tab_snap(tab);
 953
 954         sample = isl_tab_get_sample_value(tab);
 955
 956         for (i = 0; i < vertices->n_vertices; ++i) {
 957                 selection[i] = isl_basic_set_contains(vertices->v[i].dom, sample);
 958                 if (selection[i] < 0)
 959                         goto error;
 960                 if (!selection[i])
 961                         continue;
 962                 selection[i] = can_intersect(tab, vertices->v[i].dom);
 963                 if (selection[i] < 0)
 964                         goto error;
 965         }
 966
 967         if (isl_tab_detect_redundant(tab) < 0)
 968                 goto error;
 969
 970         if (add_chamber(&list, vertices, tab, selection) < 0)
 971                 goto error;
 972         n_chambers++;
 973
 974         if (init_todo(&todo, tab) < 0)
 975                 goto error;
 976
 977         while (todo) {
 978                 struct isl_facet_todo *next;
 979
 980                 if (isl_tab_rollback(tab, snap) < 0)
 981                         goto error;
 982
 983                 if (isl_tab_add_ineq(tab, todo->constraint->el) < 0)
 984                         goto error;
 985                 if (isl_tab_freeze_constraint(tab, tab->n_con - 1) < 0)
 986                         goto error;
 987
 988                 for (i = 0; i < vertices->n_vertices; ++i) {
 989                         selection[i] = bset_covers_tab(vertices->v[i].dom,
 990                                                         todo->tab);
 991                         if (selection[i] < 0)
 992                                 goto error;
 993                         if (!selection[i])
 994                                 continue;
 995                         selection[i] = can_intersect(tab, vertices->v[i].dom);
 996                         if (selection[i] < 0)
 997                                 goto error;
 998                 }
 999
1000                 if (isl_tab_detect_redundant(tab) < 0)
1001                         goto error;
1002
1003                 if (add_chamber(&list, vertices, tab, selection) < 0)
1004                         goto error;
1005                 n_chambers++;
1006
1007                 if (update_todo(todo, tab) < 0)
1008                         goto error;
1009
1010                 next = todo->next;
1011                 todo->next = NULL;
1012                 free_todo(todo);
1013                 todo = next;
1014         }
1015
1016         isl_vec_free(sample);
1017
1018         isl_tab_free(tab);
1019         free(selection);
1020
1021         vertices = vertices_add_chambers(vertices, n_chambers, list);
1022
1023         for (i = 0; vertices && i < vertices->n_vertices; ++i) {
1024                 isl_basic_set_free(vertices->v[i].dom);
1025                 vertices->v[i].dom = NULL;
1026         }
1027
1028         return vertices;
1029 error:
1030         free_chamber_list(list);
1031         free_todo(todo);
1032         isl_vec_free(sample);
1033         isl_tab_free(tab);
1034         free(selection);
1035         if (!tab)
1036                 isl_basic_set_free(bset);
1037         isl_vertices_free(vertices);
1038         return NULL;
1039 }
1040
1041 isl_ctx *isl_vertex_get_ctx(__isl_keep isl_vertex *vertex)
1042 {
1043         return vertex ? isl_vertices_get_ctx(vertex->vertices) : NULL;
1044 }
1045
1046 isl_size isl_vertex_get_id(__isl_keep isl_vertex *vertex)
1047 {
1048         return vertex ? vertex->id : isl_size_error;
1049 }
1050
1051 /* Return the activity domain of the vertex "vertex".
1052  */
1053 __isl_give isl_basic_set *isl_vertex_get_domain(__isl_keep isl_vertex *vertex)
1054 {
1055         struct isl_vertex *v;
1056
1057         if (!vertex)
1058                 return NULL;
1059
1060         v = &vertex->vertices->v[vertex->id];
1061         if (!v->dom) {
1062                 v->dom = isl_basic_set_copy(v->vertex);
1063                 v->dom = isl_basic_set_params(v->dom);
1064                 v->dom = isl_basic_set_set_integral(v->dom);
1065         }
1066
1067         return isl_basic_set_copy(v->dom);
1068 }
1069
1070 /* Return a multiple quasi-affine expression describing the vertex "vertex"
1071  * in terms of the parameters,
1072  */
1073 __isl_give isl_multi_aff *isl_vertex_get_expr(__isl_keep isl_vertex *vertex)
1074 {
1075         struct isl_vertex *v;
1076         isl_basic_set *bset;
1077
1078         if (!vertex)
1079                 return NULL;
1080
1081         v = &vertex->vertices->v[vertex->id];
1082
1083         bset = isl_basic_set_copy(v->vertex);
1084         return isl_multi_aff_from_basic_set_equalities(bset);
1085 }
1086
1087 static __isl_give isl_vertex *isl_vertex_alloc(__isl_take isl_vertices *vertices,
1088         int id)
1089 {
1090         isl_ctx *ctx;
1091         isl_vertex *vertex;
1092
1093         if (!vertices)
1094                 return NULL;
1095
1096         ctx = isl_vertices_get_ctx(vertices);
1097         vertex = isl_alloc_type(ctx, isl_vertex);
1098         if (!vertex)
1099                 goto error;
1100
1101         vertex->vertices = vertices;
1102         vertex->id = id;
1103
1104         return vertex;
1105 error:
1106         isl_vertices_free(vertices);
1107         return NULL;
1108 }
1109
1110 __isl_null isl_vertex *isl_vertex_free(__isl_take isl_vertex *vertex)
1111 {
1112         if (!vertex)
1113                 return NULL;
1114         isl_vertices_free(vertex->vertices);
1115         free(vertex);
1116
1117         return NULL;
1118 }
1119
1120 isl_ctx *isl_cell_get_ctx(__isl_keep isl_cell *cell)
1121 {
1122         return cell ? cell->dom->ctx : NULL;
1123 }
1124
1125 __isl_give isl_basic_set *isl_cell_get_domain(__isl_keep isl_cell *cell)
1126 {
1127         return cell ? isl_basic_set_copy(cell->dom) : NULL;
1128 }
1129
1130 static __isl_give isl_cell *isl_cell_alloc(__isl_take isl_vertices *vertices,
1131         __isl_take isl_basic_set *dom, int id)
1132 {
1133         int i;
1134         isl_cell *cell = NULL;
1135
1136         if (!vertices || !dom)
1137                 goto error;
1138
1139         cell = isl_calloc_type(dom->ctx, isl_cell);
1140         if (!cell)
1141                 goto error;
1142
1143         cell->n_vertices = vertices->c[id].n_vertices;
1144         cell->ids = isl_alloc_array(dom->ctx, int, cell->n_vertices);
1145         if (cell->n_vertices && !cell->ids)
1146                 goto error;
1147         for (i = 0; i < cell->n_vertices; ++i)
1148                 cell->ids[i] = vertices->c[id].vertices[i];
1149         cell->vertices = vertices;
1150         cell->dom = dom;
1151
1152         return cell;
1153 error:
1154         isl_cell_free(cell);
1155         isl_vertices_free(vertices);
1156         isl_basic_set_free(dom);
1157         return NULL;
1158 }
1159
1160 __isl_null isl_cell *isl_cell_free(__isl_take isl_cell *cell)
1161 {
1162         if (!cell)
1163                 return NULL;
1164
1165         isl_vertices_free(cell->vertices);
1166         free(cell->ids);
1167         isl_basic_set_free(cell->dom);
1168         free(cell);
1169
1170         return NULL;
1171 }
1172
1173 /* Create a tableau of the cone obtained by first homogenizing the given
1174  * polytope and then making all inequalities strict by setting the
1175  * constant term to -1.
1176  */
1177 static struct isl_tab *tab_for_shifted_cone(__isl_keep isl_basic_set *bset)
1178 {
1179         int i;
1180         isl_vec *c = NULL;
1181         struct isl_tab *tab;
1182         isl_size total;
1183
1184         total = isl_basic_set_dim(bset, isl_dim_all);
1185         if (total < 0)
1186                 return NULL;
1187         tab = isl_tab_alloc(bset->ctx, bset->n_eq + bset->n_ineq + 1,
1188                             1 + total, 0);
1189         if (!tab)
1190                 return NULL;
1191         tab->rational = ISL_F_ISSET(bset, ISL_BASIC_SET_RATIONAL);
1192         if (ISL_F_ISSET(bset, ISL_BASIC_MAP_EMPTY)) {
1193                 if (isl_tab_mark_empty(tab) < 0)
1194                         goto error;
1195                 return tab;
1196         }
1197
1198         c = isl_vec_alloc(bset->ctx, 1 + 1 + total);
1199         if (!c)
1200                 goto error;
1201
1202         isl_int_set_si(c->el[0], 0);
1203         for (i = 0; i < bset->n_eq; ++i) {
1204                 isl_seq_cpy(c->el + 1, bset->eq[i], c->size - 1);
1205                 if (isl_tab_add_eq(tab, c->el) < 0)
1206                         goto error;
1207         }
1208
1209         isl_int_set_si(c->el[0], -1);
1210         for (i = 0; i < bset->n_ineq; ++i) {
1211                 isl_seq_cpy(c->el + 1, bset->ineq[i], c->size - 1);
1212                 if (isl_tab_add_ineq(tab, c->el) < 0)
1213                         goto error;
1214                 if (tab->empty) {
1215                         isl_vec_free(c);
1216                         return tab;
1217                 }
1218         }
1219
1220         isl_seq_clr(c->el + 1, c->size - 1);
1221         isl_int_set_si(c->el[1], 1);
1222         if (isl_tab_add_ineq(tab, c->el) < 0)
1223                 goto error;
1224
1225         isl_vec_free(c);
1226         return tab;
1227 error:
1228         isl_vec_free(c);
1229         isl_tab_free(tab);
1230         return NULL;
1231 }
1232
1233 /* Compute an interior point of "bset" by selecting an interior
1234  * point in homogeneous space and projecting the point back down.
1235  */
1236 static __isl_give isl_vec *isl_basic_set_interior_point(
1237         __isl_keep isl_basic_set *bset)
1238 {
1239         isl_vec *vec;
1240         struct isl_tab *tab;
1241
1242         tab = tab_for_shifted_cone(bset);
1243         vec = isl_tab_get_sample_value(tab);
1244         isl_tab_free(tab);
1245         if (!vec)
1246                 return NULL;
1247
1248         isl_seq_cpy(vec->el, vec->el + 1, vec->size - 1);
1249         vec->size--;
1250
1251         return vec;
1252 }
1253
1254 /* Call "fn" on all chambers of the parametric polytope with the shared
1255  * facets of neighboring chambers only appearing in one of the chambers.
1256  *
1257  * We pick an interior point from one of the chambers and then make
1258  * all constraints that do not satisfy this point strict.
1259  * For constraints that saturate the interior point, the sign
1260  * of the first non-zero coefficient is used to determine which
1261  * of the two (internal) constraints should be tightened.
1262  */
1263 isl_stat isl_vertices_foreach_disjoint_cell(__isl_keep isl_vertices *vertices,
1264         isl_stat (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1265 {
1266         int i;
1267         isl_vec *vec;
1268         isl_cell *cell;
1269
1270         if (!vertices)
1271                 return isl_stat_error;
1272
1273         if (vertices->n_chambers == 0)
1274                 return isl_stat_ok;
1275
1276         if (vertices->n_chambers == 1) {
1277                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[0].dom);
1278                 dom = isl_basic_set_set_integral(dom);
1279                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, 0);
1280                 if (!cell)
1281                         return isl_stat_error;
1282                 return fn(cell, user);
1283         }
1284
1285         vec = isl_basic_set_interior_point(vertices->c[0].dom);
1286         if (!vec)
1287                 return isl_stat_error;
1288
1289         for (i = 0; i < vertices->n_chambers; ++i) {
1290                 int r;
1291                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1292                 if (i)
1293                         dom = isl_basic_set_tighten_outward(dom, vec);
1294                 dom = isl_basic_set_set_integral(dom);
1295                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1296                 if (!cell)
1297                         goto error;
1298                 r = fn(cell, user);
1299                 if (r < 0)
1300                         goto error;
1301         }
1302
1303         isl_vec_free(vec);
1304
1305         return isl_stat_ok;
1306 error:
1307         isl_vec_free(vec);
1308         return isl_stat_error;
1309 }
1310
1311 isl_stat isl_vertices_foreach_cell(__isl_keep isl_vertices *vertices,
1312         isl_stat (*fn)(__isl_take isl_cell *cell, void *user), void *user)
1313 {
1314         int i;
1315         isl_cell *cell;
1316
1317         if (!vertices)
1318                 return isl_stat_error;
1319
1320         if (vertices->n_chambers == 0)
1321                 return isl_stat_ok;
1322
1323         for (i = 0; i < vertices->n_chambers; ++i) {
1324                 isl_stat r;
1325                 isl_basic_set *dom = isl_basic_set_copy(vertices->c[i].dom);
1326
1327                 cell = isl_cell_alloc(isl_vertices_copy(vertices), dom, i);
1328                 if (!cell)
1329                         return isl_stat_error;
1330
1331                 r = fn(cell, user);
1332                 if (r < 0)
1333                         return isl_stat_error;
1334         }
1335
1336         return isl_stat_ok;
1337 }
1338
1339 isl_stat isl_vertices_foreach_vertex(__isl_keep isl_vertices *vertices,
1340         isl_stat (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1341 {
1342         int i;
1343         isl_vertex *vertex;
1344
1345         if (!vertices)
1346                 return isl_stat_error;
1347
1348         if (vertices->n_vertices == 0)
1349                 return isl_stat_ok;
1350
1351         for (i = 0; i < vertices->n_vertices; ++i) {
1352                 isl_stat r;
1353
1354                 vertex = isl_vertex_alloc(isl_vertices_copy(vertices), i);
1355                 if (!vertex)
1356                         return isl_stat_error;
1357
1358                 r = fn(vertex, user);
1359                 if (r < 0)
1360                         return isl_stat_error;
1361         }
1362
1363         return isl_stat_ok;
1364 }
1365
1366 isl_stat isl_cell_foreach_vertex(__isl_keep isl_cell *cell,
1367         isl_stat (*fn)(__isl_take isl_vertex *vertex, void *user), void *user)
1368 {
1369         int i;
1370         isl_vertex *vertex;
1371
1372         if (!cell)
1373                 return isl_stat_error;
1374
1375         if (cell->n_vertices == 0)
1376                 return isl_stat_ok;
1377
1378         for (i = 0; i < cell->n_vertices; ++i) {
1379                 isl_stat r;
1380
1381                 vertex = isl_vertex_alloc(isl_vertices_copy(cell->vertices),
1382                                           cell->ids[i]);
1383                 if (!vertex)
1384                         return isl_stat_error;
1385
1386                 r = fn(vertex, user);
1387                 if (r < 0)
1388                         return isl_stat_error;
1389         }
1390
1391         return isl_stat_ok;
1392 }
1393
1394 isl_ctx *isl_vertices_get_ctx(__isl_keep isl_vertices *vertices)
1395 {
1396         return vertices ? vertices->bset->ctx : NULL;
1397 }
1398
1399 isl_size isl_vertices_get_n_vertices(__isl_keep isl_vertices *vertices)
1400 {
1401         return vertices ? vertices->n_vertices : isl_size_error;
1402 }
1403
1404 __isl_give isl_vertices *isl_morph_vertices(__isl_take isl_morph *morph,
1405         __isl_take isl_vertices *vertices)
1406 {
1407         int i;
1408         isl_morph *param_morph = NULL;
1409
1410         if (!morph || !vertices)
1411                 goto error;
1412
1413         isl_assert(vertices->bset->ctx, vertices->ref == 1, goto error);
1414
1415         param_morph = isl_morph_copy(morph);
1416         param_morph = isl_morph_dom_params(param_morph);
1417         param_morph = isl_morph_ran_params(param_morph);
1418
1419         for (i = 0; i < vertices->n_vertices; ++i) {
1420                 vertices->v[i].dom = isl_morph_basic_set(
1421                         isl_morph_copy(param_morph), vertices->v[i].dom);
1422                 vertices->v[i].vertex = isl_morph_basic_set(
1423                         isl_morph_copy(morph), vertices->v[i].vertex);
1424                 if (!vertices->v[i].vertex)
1425                         goto error;
1426         }
1427
1428         for (i = 0; i < vertices->n_chambers; ++i) {
1429                 vertices->c[i].dom = isl_morph_basic_set(
1430                         isl_morph_copy(param_morph), vertices->c[i].dom);
1431                 if (!vertices->c[i].dom)
1432                         goto error;
1433         }
1434
1435         isl_morph_free(param_morph);
1436         isl_morph_free(morph);
1437         return vertices;
1438 error:
1439         isl_morph_free(param_morph);
1440         isl_morph_free(morph);
1441         isl_vertices_free(vertices);
1442         return NULL;
1443 }
1444
1445 /* Construct a simplex isl_cell spanned by the vertices with indices in
1446  * "simplex_ids" and "other_ids" and call "fn" on this isl_cell.
1447  */
1448 static isl_stat call_on_simplex(__isl_keep isl_cell *cell,
1449         int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1450         isl_stat (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1451 {
1452         int i;
1453         isl_ctx *ctx;
1454         struct isl_cell *simplex;
1455
1456         ctx = isl_cell_get_ctx(cell);
1457
1458         simplex = isl_calloc_type(ctx, struct isl_cell);
1459         if (!simplex)
1460                 return isl_stat_error;
1461         simplex->vertices = isl_vertices_copy(cell->vertices);
1462         if (!simplex->vertices)
1463                 goto error;
1464         simplex->dom = isl_basic_set_copy(cell->dom);
1465         if (!simplex->dom)
1466                 goto error;
1467         simplex->n_vertices = n_simplex + n_other;
1468         simplex->ids = isl_alloc_array(ctx, int, simplex->n_vertices);
1469         if (!simplex->ids)
1470                 goto error;
1471
1472         for (i = 0; i < n_simplex; ++i)
1473                 simplex->ids[i] = simplex_ids[i];
1474         for (i = 0; i < n_other; ++i)
1475                 simplex->ids[n_simplex + i] = other_ids[i];
1476
1477         return fn(simplex, user);
1478 error:
1479         isl_cell_free(simplex);
1480         return isl_stat_error;
1481 }
1482
1483 /* Check whether the parametric vertex described by "vertex"
1484  * lies on the facet corresponding to constraint "facet" of "bset".
1485  * The isl_vec "v" is a temporary vector than can be used by this function.
1486  *
1487  * We eliminate the variables from the facet constraint using the
1488  * equalities defining the vertex and check if the result is identical
1489  * to zero.
1490  *
1491  * It would probably be better to keep track of the constraints defining
1492  * a vertex during the vertex construction so that we could simply look
1493  * it up here.
1494  */
1495 static int vertex_on_facet(__isl_keep isl_basic_set *vertex,
1496         __isl_keep isl_basic_set *bset, int facet, __isl_keep isl_vec *v)
1497 {
1498         int i;
1499         isl_int m;
1500
1501         isl_seq_cpy(v->el, bset->ineq[facet], v->size);
1502
1503         isl_int_init(m);
1504         for (i = 0; i < vertex->n_eq; ++i) {
1505                 int k = isl_seq_last_non_zero(vertex->eq[i], v->size);
1506                 isl_seq_elim(v->el, vertex->eq[i], k, v->size, &m);
1507         }
1508         isl_int_clear(m);
1509
1510         return isl_seq_first_non_zero(v->el, v->size) == -1;
1511 }
1512
1513 /* Triangulate the polytope spanned by the vertices with ids
1514  * in "simplex_ids" and "other_ids" and call "fn" on each of
1515  * the resulting simplices.
1516  * If the input polytope is already a simplex, we simply call "fn".
1517  * Otherwise, we pick a point from "other_ids" and add it to "simplex_ids".
1518  * Then we consider each facet of "bset" that does not contain the point
1519  * we just picked, but does contain some of the other points in "other_ids"
1520  * and call ourselves recursively on the polytope spanned by the new
1521  * "simplex_ids" and those points in "other_ids" that lie on the facet.
1522  */
1523 static isl_stat triangulate(__isl_keep isl_cell *cell, __isl_keep isl_vec *v,
1524         int *simplex_ids, int n_simplex, int *other_ids, int n_other,
1525         isl_stat (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1526 {
1527         int i, j, k;
1528         isl_size d, nparam;
1529         int *ids;
1530         isl_ctx *ctx;
1531         isl_basic_set *vertex;
1532         isl_basic_set *bset;
1533
1534         ctx = isl_cell_get_ctx(cell);
1535         d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1536         nparam = isl_basic_set_dim(cell->vertices->bset, isl_dim_param);
1537         if (d < 0 || nparam < 0)
1538                 return isl_stat_error;
1539
1540         if (n_simplex + n_other == d + 1)
1541                 return call_on_simplex(cell, simplex_ids, n_simplex,
1542                                        other_ids, n_other, fn, user);
1543
1544         simplex_ids[n_simplex] = other_ids[0];
1545         vertex = cell->vertices->v[other_ids[0]].vertex;
1546         bset = cell->vertices->bset;
1547
1548         ids = isl_alloc_array(ctx, int, n_other - 1);
1549         if (!ids)
1550                 goto error;
1551         for (i = 0; i < bset->n_ineq; ++i) {
1552                 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + nparam, d) == -1)
1553                         continue;
1554                 if (vertex_on_facet(vertex, bset, i, v))
1555                         continue;
1556
1557                 for (j = 1, k = 0; j < n_other; ++j) {
1558                         isl_basic_set *ov;
1559                         ov = cell->vertices->v[other_ids[j]].vertex;
1560                         if (vertex_on_facet(ov, bset, i, v))
1561                                 ids[k++] = other_ids[j];
1562                 }
1563                 if (k == 0)
1564                         continue;
1565
1566                 if (triangulate(cell, v, simplex_ids, n_simplex + 1,
1567                                 ids, k, fn, user) < 0)
1568                         goto error;
1569         }
1570         free(ids);
1571
1572         return isl_stat_ok;
1573 error:
1574         free(ids);
1575         return isl_stat_error;
1576 }
1577
1578 /* Triangulate the given cell and call "fn" on each of the resulting
1579  * simplices.
1580  */
1581 isl_stat isl_cell_foreach_simplex(__isl_take isl_cell *cell,
1582         isl_stat (*fn)(__isl_take isl_cell *simplex, void *user), void *user)
1583 {
1584         isl_size d, total;
1585         isl_stat r;
1586         isl_ctx *ctx;
1587         isl_vec *v = NULL;
1588         int *simplex_ids = NULL;
1589
1590         if (!cell)
1591                 return isl_stat_error;
1592
1593         d = isl_basic_set_dim(cell->vertices->bset, isl_dim_set);
1594         total = isl_basic_set_dim(cell->vertices->bset, isl_dim_all);
1595         if (d < 0 || total < 0)
1596                 return isl_stat_error;
1597
1598         if (cell->n_vertices == d + 1)
1599                 return fn(cell, user);
1600
1601         ctx = isl_cell_get_ctx(cell);
1602         simplex_ids = isl_alloc_array(ctx, int, d + 1);
1603         if (!simplex_ids)
1604                 goto error;
1605
1606         v = isl_vec_alloc(ctx, 1 + total);
1607         if (!v)
1608                 goto error;
1609
1610         r = triangulate(cell, v, simplex_ids, 0,
1611                         cell->ids, cell->n_vertices, fn, user);
1612
1613         isl_vec_free(v);
1614         free(simplex_ids);
1615
1616         isl_cell_free(cell);
1617
1618         return r;
1619 error:
1620         free(simplex_ids);
1621         isl_vec_free(v);
1622         isl_cell_free(cell);
1623         return isl_stat_error;
1624 }