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