isl_vertices.c

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