isl_affine_hull.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include "isl_ctx.h"
  11 #include "isl_seq.h"
  12 #include "isl_set.h"
  13 #include "isl_lp.h"
  14 #include "isl_map.h"
  15 #include "isl_map_private.h"
  16 #include "isl_equalities.h"
  17 #include "isl_sample.h"
  18 #include "isl_tab.h"
  19
  20 struct isl_basic_map *isl_basic_map_implicit_equalities(
  21                                                 struct isl_basic_map *bmap)
  22 {
  23         struct isl_tab *tab;
  24
  25         if (!bmap)
  26                 return bmap;
  27
  28         bmap = isl_basic_map_gauss(bmap, NULL);
  29         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
  30                 return bmap;
  31         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_IMPLICIT))
  32                 return bmap;
  33         if (bmap->n_ineq <= 1)
  34                 return bmap;
  35
  36         tab = isl_tab_from_basic_map(bmap);
  37         tab = isl_tab_detect_implicit_equalities(tab);
  38         bmap = isl_basic_map_update_from_tab(bmap, tab);
  39         isl_tab_free(tab);
  40         bmap = isl_basic_map_gauss(bmap, NULL);
  41         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
  42         return bmap;
  43 }
  44
  45 struct isl_basic_set *isl_basic_set_implicit_equalities(
  46                                                 struct isl_basic_set *bset)
  47 {
  48         return (struct isl_basic_set *)
  49                 isl_basic_map_implicit_equalities((struct isl_basic_map*)bset);
  50 }
  51
  52 struct isl_map *isl_map_implicit_equalities(struct isl_map *map)
  53 {
  54         int i;
  55
  56         if (!map)
  57                 return map;
  58
  59         for (i = 0; i < map->n; ++i) {
  60                 map->p[i] = isl_basic_map_implicit_equalities(map->p[i]);
  61                 if (!map->p[i])
  62                         goto error;
  63         }
  64
  65         return map;
  66 error:
  67         isl_map_free(map);
  68         return NULL;
  69 }
  70
  71 /* Make eq[row][col] of both bmaps equal so we can add the row
  72  * add the column to the common matrix.
  73  * Note that because of the echelon form, the columns of row row
  74  * after column col are zero.
  75  */
  76 static void set_common_multiple(
  77         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
  78         unsigned row, unsigned col)
  79 {
  80         isl_int m, c;
  81
  82         if (isl_int_eq(bset1->eq[row][col], bset2->eq[row][col]))
  83                 return;
  84
  85         isl_int_init(c);
  86         isl_int_init(m);
  87         isl_int_lcm(m, bset1->eq[row][col], bset2->eq[row][col]);
  88         isl_int_divexact(c, m, bset1->eq[row][col]);
  89         isl_seq_scale(bset1->eq[row], bset1->eq[row], c, col+1);
  90         isl_int_divexact(c, m, bset2->eq[row][col]);
  91         isl_seq_scale(bset2->eq[row], bset2->eq[row], c, col+1);
  92         isl_int_clear(c);
  93         isl_int_clear(m);
  94 }
  95
  96 /* Delete a given equality, moving all the following equalities one up.
  97  */
  98 static void delete_row(struct isl_basic_set *bset, unsigned row)
  99 {
 100         isl_int *t;
 101         int r;
 102
 103         t = bset->eq[row];
 104         bset->n_eq--;
 105         for (r = row; r < bset->n_eq; ++r)
 106                 bset->eq[r] = bset->eq[r+1];
 107         bset->eq[bset->n_eq] = t;
 108 }
 109
 110 /* Make first row entries in column col of bset1 identical to
 111  * those of bset2, using the fact that entry bset1->eq[row][col]=a
 112  * is non-zero.  Initially, these elements of bset1 are all zero.
 113  * For each row i < row, we set
 114  *              A[i] = a * A[i] + B[i][col] * A[row]
 115  *              B[i] = a * B[i]
 116  * so that
 117  *              A[i][col] = B[i][col] = a * old(B[i][col])
 118  */
 119 static void construct_column(
 120         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 121         unsigned row, unsigned col)
 122 {
 123         int r;
 124         isl_int a;
 125         isl_int b;
 126         unsigned total;
 127
 128         isl_int_init(a);
 129         isl_int_init(b);
 130         total = 1 + isl_basic_set_n_dim(bset1);
 131         for (r = 0; r < row; ++r) {
 132                 if (isl_int_is_zero(bset2->eq[r][col]))
 133                         continue;
 134                 isl_int_gcd(b, bset2->eq[r][col], bset1->eq[row][col]);
 135                 isl_int_divexact(a, bset1->eq[row][col], b);
 136                 isl_int_divexact(b, bset2->eq[r][col], b);
 137                 isl_seq_combine(bset1->eq[r], a, bset1->eq[r],
 138                                               b, bset1->eq[row], total);
 139                 isl_seq_scale(bset2->eq[r], bset2->eq[r], a, total);
 140         }
 141         isl_int_clear(a);
 142         isl_int_clear(b);
 143         delete_row(bset1, row);
 144 }
 145
 146 /* Make first row entries in column col of bset1 identical to
 147  * those of bset2, using only these entries of the two matrices.
 148  * Let t be the last row with different entries.
 149  * For each row i < t, we set
 150  *      A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
 151  *      B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
 152  * so that
 153  *      A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
 154  */
 155 static int transform_column(
 156         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 157         unsigned row, unsigned col)
 158 {
 159         int i, t;
 160         isl_int a, b, g;
 161         unsigned total;
 162
 163         for (t = row-1; t >= 0; --t)
 164                 if (isl_int_ne(bset1->eq[t][col], bset2->eq[t][col]))
 165                         break;
 166         if (t < 0)
 167                 return 0;
 168
 169         total = 1 + isl_basic_set_n_dim(bset1);
 170         isl_int_init(a);
 171         isl_int_init(b);
 172         isl_int_init(g);
 173         isl_int_sub(b, bset1->eq[t][col], bset2->eq[t][col]);
 174         for (i = 0; i < t; ++i) {
 175                 isl_int_sub(a, bset2->eq[i][col], bset1->eq[i][col]);
 176                 isl_int_gcd(g, a, b);
 177                 isl_int_divexact(a, a, g);
 178                 isl_int_divexact(g, b, g);
 179                 isl_seq_combine(bset1->eq[i], g, bset1->eq[i], a, bset1->eq[t],
 180                                 total);
 181                 isl_seq_combine(bset2->eq[i], g, bset2->eq[i], a, bset2->eq[t],
 182                                 total);
 183         }
 184         isl_int_clear(a);
 185         isl_int_clear(b);
 186         isl_int_clear(g);
 187         delete_row(bset1, t);
 188         delete_row(bset2, t);
 189         return 1;
 190 }
 191
 192 /* The implementation is based on Section 5.2 of Michael Karr,
 193  * "Affine Relationships Among Variables of a Program",
 194  * except that the echelon form we use starts from the last column
 195  * and that we are dealing with integer coefficients.
 196  */
 197 static struct isl_basic_set *affine_hull(
 198         struct isl_basic_set *bset1, struct isl_basic_set *bset2)
 199 {
 200         unsigned total;
 201         int col;
 202         int row;
 203
 204         total = 1 + isl_basic_set_n_dim(bset1);
 205
 206         row = 0;
 207         for (col = total-1; col >= 0; --col) {
 208                 int is_zero1 = row >= bset1->n_eq ||
 209                         isl_int_is_zero(bset1->eq[row][col]);
 210                 int is_zero2 = row >= bset2->n_eq ||
 211                         isl_int_is_zero(bset2->eq[row][col]);
 212                 if (!is_zero1 && !is_zero2) {
 213                         set_common_multiple(bset1, bset2, row, col);
 214                         ++row;
 215                 } else if (!is_zero1 && is_zero2) {
 216                         construct_column(bset1, bset2, row, col);
 217                 } else if (is_zero1 && !is_zero2) {
 218                         construct_column(bset2, bset1, row, col);
 219                 } else {
 220                         if (transform_column(bset1, bset2, row, col))
 221                                 --row;
 222                 }
 223         }
 224         isl_basic_set_free(bset2);
 225         isl_assert(bset1->ctx, row == bset1->n_eq, goto error);
 226         bset1 = isl_basic_set_normalize_constraints(bset1);
 227         return bset1;
 228 error:
 229         isl_basic_set_free(bset1);
 230         return NULL;
 231 }
 232
 233 /* Find an integer point in the set represented by "tab"
 234  * that lies outside of the equality "eq" e(x) = 0.
 235  * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
 236  * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
 237  * The point, if found, is returned.
 238  * If no point can be found, a zero-length vector is returned.
 239  *
 240  * Before solving an ILP problem, we first check if simply
 241  * adding the normal of the constraint to one of the known
 242  * integer points in the basic set represented by "tab"
 243  * yields another point inside the basic set.
 244  *
 245  * The caller of this function ensures that the tableau is bounded or
 246  * that tab->basis and tab->n_unbounded have been set appropriately.
 247  */
 248 static struct isl_vec *outside_point(struct isl_tab *tab, isl_int *eq, int up)
 249 {
 250         struct isl_ctx *ctx;
 251         struct isl_vec *sample = NULL;
 252         struct isl_tab_undo *snap;
 253         unsigned dim;
 254         int k;
 255
 256         if (!tab)
 257                 return NULL;
 258         ctx = tab->mat->ctx;
 259
 260         dim = tab->n_var;
 261         sample = isl_vec_alloc(ctx, 1 + dim);
 262         if (!sample)
 263                 return NULL;
 264         isl_int_set_si(sample->el[0], 1);
 265         isl_seq_combine(sample->el + 1,
 266                 ctx->one, tab->bmap->sample->el + 1,
 267                 up ? ctx->one : ctx->negone, eq + 1, dim);
 268         if (isl_basic_map_contains(tab->bmap, sample))
 269                 return sample;
 270         isl_vec_free(sample);
 271         sample = NULL;
 272
 273         snap = isl_tab_snap(tab);
 274
 275         if (!up)
 276                 isl_seq_neg(eq, eq, 1 + dim);
 277         isl_int_sub_ui(eq[0], eq[0], 1);
 278
 279         if (isl_tab_extend_cons(tab, 1) < 0)
 280                 goto error;
 281         if (isl_tab_add_ineq(tab, eq) < 0)
 282                 goto error;
 283
 284         sample = isl_tab_sample(tab);
 285
 286         isl_int_add_ui(eq[0], eq[0], 1);
 287         if (!up)
 288                 isl_seq_neg(eq, eq, 1 + dim);
 289
 290         if (isl_tab_rollback(tab, snap) < 0)
 291                 goto error;
 292
 293         return sample;
 294 error:
 295         isl_vec_free(sample);
 296         return NULL;
 297 }
 298
 299 struct isl_basic_set *isl_basic_set_recession_cone(struct isl_basic_set *bset)
 300 {
 301         int i;
 302
 303         bset = isl_basic_set_cow(bset);
 304         if (!bset)
 305                 return NULL;
 306         isl_assert(bset->ctx, bset->n_div == 0, goto error);
 307
 308         for (i = 0; i < bset->n_eq; ++i)
 309                 isl_int_set_si(bset->eq[i][0], 0);
 310
 311         for (i = 0; i < bset->n_ineq; ++i)
 312                 isl_int_set_si(bset->ineq[i][0], 0);
 313
 314         ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
 315         return isl_basic_set_implicit_equalities(bset);
 316 error:
 317         isl_basic_set_free(bset);
 318         return NULL;
 319 }
 320
 321 /* Extend an initial (under-)approximation of the affine hull of basic
 322  * set represented by the tableau "tab"
 323  * by looking for points that do not satisfy one of the equalities
 324  * in the current approximation and adding them to that approximation
 325  * until no such points can be found any more.
 326  *
 327  * The caller of this function ensures that "tab" is bounded or
 328  * that tab->basis and tab->n_unbounded have been set appropriately.
 329  */
 330 static struct isl_basic_set *extend_affine_hull(struct isl_tab *tab,
 331         struct isl_basic_set *hull)
 332 {
 333         int i, j, k;
 334         unsigned dim;
 335
 336         if (!tab || !hull)
 337                 goto error;
 338
 339         dim = tab->n_var;
 340
 341         if (isl_tab_extend_cons(tab, 2 * dim + 1) < 0)
 342                 goto error;
 343
 344         for (i = 0; i < dim; ++i) {
 345                 struct isl_vec *sample;
 346                 struct isl_basic_set *point;
 347                 for (j = 0; j < hull->n_eq; ++j) {
 348                         sample = outside_point(tab, hull->eq[j], 1);
 349                         if (!sample)
 350                                 goto error;
 351                         if (sample->size > 0)
 352                                 break;
 353                         isl_vec_free(sample);
 354                         sample = outside_point(tab, hull->eq[j], 0);
 355                         if (!sample)
 356                                 goto error;
 357                         if (sample->size > 0)
 358                                 break;
 359                         isl_vec_free(sample);
 360
 361                         tab = isl_tab_add_eq(tab, hull->eq[j]);
 362                         if (!tab)
 363                                 goto error;
 364                 }
 365                 if (j == hull->n_eq)
 366                         break;
 367                 if (tab->samples)
 368                         tab = isl_tab_add_sample(tab, isl_vec_copy(sample));
 369                 if (!tab)
 370                         goto error;
 371                 point = isl_basic_set_from_vec(sample);
 372                 hull = affine_hull(hull, point);
 373         }
 374
 375         return hull;
 376 error:
 377         isl_basic_set_free(hull);
 378         return NULL;
 379 }
 380
 381 /* Drop all constraints in bset that involve any of the dimensions
 382  * first to first+n-1.
 383  */
 384 static struct isl_basic_set *drop_constraints_involving
 385         (struct isl_basic_set *bset, unsigned first, unsigned n)
 386 {
 387         int i;
 388
 389         if (!bset)
 390                 return NULL;
 391
 392         bset = isl_basic_set_cow(bset);
 393
 394         for (i = bset->n_eq - 1; i >= 0; --i) {
 395                 if (isl_seq_first_non_zero(bset->eq[i] + 1 + first, n) == -1)
 396                         continue;
 397                 isl_basic_set_drop_equality(bset, i);
 398         }
 399
 400         for (i = bset->n_ineq - 1; i >= 0; --i) {
 401                 if (isl_seq_first_non_zero(bset->ineq[i] + 1 + first, n) == -1)
 402                         continue;
 403                 isl_basic_set_drop_inequality(bset, i);
 404         }
 405
 406         return bset;
 407 }
 408
 409 /* Look for all equalities satisfied by the integer points in bset,
 410  * which is assumed to be bounded.
 411  *
 412  * The equalities are obtained by successively looking for
 413  * a point that is affinely independent of the points found so far.
 414  * In particular, for each equality satisfied by the points so far,
 415  * we check if there is any point on a hyperplane parallel to the
 416  * corresponding hyperplane shifted by at least one (in either direction).
 417  */
 418 static struct isl_basic_set *uset_affine_hull_bounded(struct isl_basic_set *bset)
 419 {
 420         struct isl_vec *sample = NULL;
 421         struct isl_basic_set *hull;
 422         struct isl_tab *tab = NULL;
 423         unsigned dim;
 424
 425         if (isl_basic_set_fast_is_empty(bset))
 426                 return bset;
 427
 428         dim = isl_basic_set_n_dim(bset);
 429
 430         if (bset->sample && bset->sample->size == 1 + dim) {
 431                 int contains = isl_basic_set_contains(bset, bset->sample);
 432                 if (contains < 0)
 433                         goto error;
 434                 if (contains) {
 435                         if (dim == 0)
 436                                 return bset;
 437                         sample = isl_vec_copy(bset->sample);
 438                 } else {
 439                         isl_vec_free(bset->sample);
 440                         bset->sample = NULL;
 441                 }
 442         }
 443
 444         tab = isl_tab_from_basic_set(bset);
 445         if (!tab)
 446                 goto error;
 447         if (tab->empty) {
 448                 isl_tab_free(tab);
 449                 isl_vec_free(sample);
 450                 return isl_basic_set_set_to_empty(bset);
 451         }
 452         if (isl_tab_track_bset(tab, isl_basic_set_copy(bset)) < 0)
 453                 goto error;
 454
 455         if (!sample) {
 456                 struct isl_tab_undo *snap;
 457                 snap = isl_tab_snap(tab);
 458                 sample = isl_tab_sample(tab);
 459                 if (isl_tab_rollback(tab, snap) < 0)
 460                         goto error;
 461                 isl_vec_free(tab->bmap->sample);
 462                 tab->bmap->sample = isl_vec_copy(sample);
 463         }
 464
 465         if (!sample)
 466                 goto error;
 467         if (sample->size == 0) {
 468                 isl_tab_free(tab);
 469                 isl_vec_free(sample);
 470                 return isl_basic_set_set_to_empty(bset);
 471         }
 472
 473         hull = isl_basic_set_from_vec(sample);
 474
 475         isl_basic_set_free(bset);
 476         hull = extend_affine_hull(tab, hull);
 477         isl_tab_free(tab);
 478
 479         return hull;
 480 error:
 481         isl_vec_free(sample);
 482         isl_tab_free(tab);
 483         isl_basic_set_free(bset);
 484         return NULL;
 485 }
 486
 487 /* Given an unbounded tableau and an integer point satisfying the tableau,
 488  * construct an intial affine hull containing the recession cone
 489  * shifted to the given point.
 490  *
 491  * The unbounded directions are taken from the last rows of the basis,
 492  * which is assumed to have been initialized appropriately.
 493  */
 494 static __isl_give isl_basic_set *initial_hull(struct isl_tab *tab,
 495         __isl_take isl_vec *vec)
 496 {
 497         int i;
 498         int k;
 499         struct isl_basic_set *bset = NULL;
 500         struct isl_ctx *ctx;
 501         unsigned dim;
 502
 503         if (!vec || !tab)
 504                 return NULL;
 505         ctx = vec->ctx;
 506         isl_assert(ctx, vec->size != 0, goto error);
 507
 508         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
 509         if (!bset)
 510                 goto error;
 511         dim = isl_basic_set_n_dim(bset) - tab->n_unbounded;
 512         for (i = 0; i < dim; ++i) {
 513                 k = isl_basic_set_alloc_equality(bset);
 514                 if (k < 0)
 515                         goto error;
 516                 isl_seq_cpy(bset->eq[k] + 1, tab->basis->row[1 + i] + 1,
 517                             vec->size - 1);
 518                 isl_seq_inner_product(bset->eq[k] + 1, vec->el +1,
 519                                       vec->size - 1, &bset->eq[k][0]);
 520                 isl_int_neg(bset->eq[k][0], bset->eq[k][0]);
 521         }
 522         bset->sample = vec;
 523         bset = isl_basic_set_gauss(bset, NULL);
 524
 525         return bset;
 526 error:
 527         isl_basic_set_free(bset);
 528         isl_vec_free(vec);
 529         return NULL;
 530 }
 531
 532 /* Given a tableau of a set and a tableau of the corresponding
 533  * recession cone, detect and add all equalities to the tableau.
 534  * If the tableau is bounded, then we can simply keep the
 535  * tableau in its state after the return from extend_affine_hull.
 536  * However, if the tableau is unbounded, then
 537  * isl_tab_set_initial_basis_with_cone will add some additional
 538  * constraints to the tableau that have to be removed again.
 539  * In this case, we therefore rollback to the state before
 540  * any constraints were added and then add the eqaulities back in.
 541  */
 542 struct isl_tab *isl_tab_detect_equalities(struct isl_tab *tab,
 543         struct isl_tab *tab_cone)
 544 {
 545         int j;
 546         struct isl_vec *sample;
 547         struct isl_basic_set *hull;
 548         struct isl_tab_undo *snap;
 549
 550         if (!tab || !tab_cone)
 551                 goto error;
 552
 553         snap = isl_tab_snap(tab);
 554
 555         isl_mat_free(tab->basis);
 556         tab->basis = NULL;
 557
 558         isl_assert(tab->mat->ctx, tab->bmap, goto error);
 559         isl_assert(tab->mat->ctx, tab->samples, goto error);
 560         isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error);
 561         isl_assert(tab->mat->ctx, tab->n_sample > tab->n_outside, goto error);
 562
 563         if (isl_tab_set_initial_basis_with_cone(tab, tab_cone) < 0)
 564                 goto error;
 565
 566         sample = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var);
 567         if (!sample)
 568                 goto error;
 569
 570         isl_seq_cpy(sample->el, tab->samples->row[tab->n_outside], sample->size);
 571
 572         isl_vec_free(tab->bmap->sample);
 573         tab->bmap->sample = isl_vec_copy(sample);
 574
 575         if (tab->n_unbounded == 0)
 576                 hull = isl_basic_set_from_vec(isl_vec_copy(sample));
 577         else
 578                 hull = initial_hull(tab, isl_vec_copy(sample));
 579
 580         for (j = tab->n_outside + 1; j < tab->n_sample; ++j) {
 581                 isl_seq_cpy(sample->el, tab->samples->row[j], sample->size);
 582                 hull = affine_hull(hull,
 583                                 isl_basic_set_from_vec(isl_vec_copy(sample)));
 584         }
 585
 586         isl_vec_free(sample);
 587
 588         hull = extend_affine_hull(tab, hull);
 589         if (!hull)
 590                 goto error;
 591
 592         if (tab->n_unbounded == 0) {
 593                 isl_basic_set_free(hull);
 594                 return tab;
 595         }
 596
 597         if (isl_tab_rollback(tab, snap) < 0)
 598                 goto error;
 599
 600         if (hull->n_eq > tab->n_zero) {
 601                 for (j = 0; j < hull->n_eq; ++j) {
 602                         isl_seq_normalize(tab->mat->ctx, hull->eq[j], 1 + tab->n_var);
 603                         tab = isl_tab_add_eq(tab, hull->eq[j]);
 604                 }
 605         }
 606
 607         isl_basic_set_free(hull);
 608
 609         return tab;
 610 error:
 611         isl_tab_free(tab);
 612         return NULL;
 613 }
 614
 615 /* Compute the affine hull of "bset", where "cone" is the recession cone
 616  * of "bset".
 617  *
 618  * We first compute a unimodular transformation that puts the unbounded
 619  * directions in the last dimensions.  In particular, we take a transformation
 620  * that maps all equalities to equalities (in HNF) on the first dimensions.
 621  * Let x be the original dimensions and y the transformed, with y_1 bounded
 622  * and y_2 unbounded.
 623  *
 624  *             [ y_1 ]                  [ y_1 ]   [ Q_1 ]
 625  *      x = U  [ y_2 ]                  [ y_2 ] = [ Q_2 ] x
 626  *
 627  * Let's call the input basic set S.  We compute S' = preimage(S, U)
 628  * and drop the final dimensions including any constraints involving them.
 629  * This results in set S''.
 630  * Then we compute the affine hull A'' of S''.
 631  * Let F y_1 >= g be the constraint system of A''.  In the transformed
 632  * space the y_2 are unbounded, so we can add them back without any constraints,
 633  * resulting in
 634  *
 635  *                      [ y_1 ]
 636  *              [ F 0 ] [ y_2 ] >= g
 637  * or
 638  *                      [ Q_1 ]
 639  *              [ F 0 ] [ Q_2 ] x >= g
 640  * or
 641  *              F Q_1 x >= g
 642  *
 643  * The affine hull in the original space is then obtained as
 644  * A = preimage(A'', Q_1).
 645  */
 646 static struct isl_basic_set *affine_hull_with_cone(struct isl_basic_set *bset,
 647         struct isl_basic_set *cone)
 648 {
 649         unsigned total;
 650         unsigned cone_dim;
 651         struct isl_basic_set *hull;
 652         struct isl_mat *M, *U, *Q;
 653
 654         if (!bset || !cone)
 655                 goto error;
 656
 657         total = isl_basic_set_total_dim(cone);
 658         cone_dim = total - cone->n_eq;
 659
 660         M = isl_mat_sub_alloc(bset->ctx, cone->eq, 0, cone->n_eq, 1, total);
 661         M = isl_mat_left_hermite(M, 0, &U, &Q);
 662         if (!M)
 663                 goto error;
 664         isl_mat_free(M);
 665
 666         U = isl_mat_lin_to_aff(U);
 667         bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
 668
 669         bset = drop_constraints_involving(bset, total - cone_dim, cone_dim);
 670         bset = isl_basic_set_drop_dims(bset, total - cone_dim, cone_dim);
 671
 672         Q = isl_mat_lin_to_aff(Q);
 673         Q = isl_mat_drop_rows(Q, 1 + total - cone_dim, cone_dim);
 674
 675         if (bset && bset->sample && bset->sample->size == 1 + total)
 676                 bset->sample = isl_mat_vec_product(isl_mat_copy(Q), bset->sample);
 677
 678         hull = uset_affine_hull_bounded(bset);
 679
 680         if (!hull)
 681                 isl_mat_free(U);
 682         else {
 683                 struct isl_vec *sample = isl_vec_copy(hull->sample);
 684                 U = isl_mat_drop_cols(U, 1 + total - cone_dim, cone_dim);
 685                 if (sample && sample->size > 0)
 686                         sample = isl_mat_vec_product(U, sample);
 687                 else
 688                         isl_mat_free(U);
 689                 hull = isl_basic_set_preimage(hull, Q);
 690                 isl_vec_free(hull->sample);
 691                 hull->sample = sample;
 692         }
 693
 694         isl_basic_set_free(cone);
 695
 696         return hull;
 697 error:
 698         isl_basic_set_free(bset);
 699         isl_basic_set_free(cone);
 700         return NULL;
 701 }
 702
 703 /* Look for all equalities satisfied by the integer points in bset,
 704  * which is assumed not to have any explicit equalities.
 705  *
 706  * The equalities are obtained by successively looking for
 707  * a point that is affinely independent of the points found so far.
 708  * In particular, for each equality satisfied by the points so far,
 709  * we check if there is any point on a hyperplane parallel to the
 710  * corresponding hyperplane shifted by at least one (in either direction).
 711  *
 712  * Before looking for any outside points, we first compute the recession
 713  * cone.  The directions of this recession cone will always be part
 714  * of the affine hull, so there is no need for looking for any points
 715  * in these directions.
 716  * In particular, if the recession cone is full-dimensional, then
 717  * the affine hull is simply the whole universe.
 718  */
 719 static struct isl_basic_set *uset_affine_hull(struct isl_basic_set *bset)
 720 {
 721         struct isl_basic_set *cone;
 722
 723         if (isl_basic_set_fast_is_empty(bset))
 724                 return bset;
 725
 726         cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
 727         if (!cone)
 728                 goto error;
 729         if (cone->n_eq == 0) {
 730                 struct isl_basic_set *hull;
 731                 isl_basic_set_free(cone);
 732                 hull = isl_basic_set_universe_like(bset);
 733                 isl_basic_set_free(bset);
 734                 return hull;
 735         }
 736
 737         if (cone->n_eq < isl_basic_set_total_dim(cone))
 738                 return affine_hull_with_cone(bset, cone);
 739
 740         isl_basic_set_free(cone);
 741         return uset_affine_hull_bounded(bset);
 742 error:
 743         isl_basic_set_free(bset);
 744         return NULL;
 745 }
 746
 747 /* Look for all equalities satisfied by the integer points in bmap
 748  * that are independent of the equalities already explicitly available
 749  * in bmap.
 750  *
 751  * We first remove all equalities already explicitly available,
 752  * then look for additional equalities in the reduced space
 753  * and then transform the result to the original space.
 754  * The original equalities are _not_ added to this set.  This is
 755  * the responsibility of the calling function.
 756  * The resulting basic set has all meaning about the dimensions removed.
 757  * In particular, dimensions that correspond to existential variables
 758  * in bmap and that are found to be fixed are not removed.
 759  */
 760 static struct isl_basic_set *equalities_in_underlying_set(
 761                                                 struct isl_basic_map *bmap)
 762 {
 763         struct isl_mat *T1 = NULL;
 764         struct isl_mat *T2 = NULL;
 765         struct isl_basic_set *bset = NULL;
 766         struct isl_basic_set *hull = NULL;
 767
 768         bset = isl_basic_map_underlying_set(bmap);
 769         if (!bset)
 770                 return NULL;
 771         if (bset->n_eq)
 772                 bset = isl_basic_set_remove_equalities(bset, &T1, &T2);
 773         if (!bset)
 774                 goto error;
 775
 776         hull = uset_affine_hull(bset);
 777         if (!T2)
 778                 return hull;
 779
 780         if (!hull)
 781                 isl_mat_free(T1);
 782         else {
 783                 struct isl_vec *sample = isl_vec_copy(hull->sample);
 784                 if (sample && sample->size > 0)
 785                         sample = isl_mat_vec_product(T1, sample);
 786                 else
 787                         isl_mat_free(T1);
 788                 hull = isl_basic_set_preimage(hull, T2);
 789                 isl_vec_free(hull->sample);
 790                 hull->sample = sample;
 791         }
 792
 793         return hull;
 794 error:
 795         isl_mat_free(T2);
 796         isl_basic_set_free(bset);
 797         isl_basic_set_free(hull);
 798         return NULL;
 799 }
 800
 801 /* Detect and make explicit all equalities satisfied by the (integer)
 802  * points in bmap.
 803  */
 804 struct isl_basic_map *isl_basic_map_detect_equalities(
 805                                                 struct isl_basic_map *bmap)
 806 {
 807         int i, j;
 808         struct isl_basic_set *hull = NULL;
 809
 810         if (!bmap)
 811                 return NULL;
 812         if (bmap->n_ineq == 0)
 813                 return bmap;
 814         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
 815                 return bmap;
 816         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_ALL_EQUALITIES))
 817                 return bmap;
 818         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
 819                 return isl_basic_map_implicit_equalities(bmap);
 820
 821         hull = equalities_in_underlying_set(isl_basic_map_copy(bmap));
 822         if (!hull)
 823                 goto error;
 824         if (ISL_F_ISSET(hull, ISL_BASIC_SET_EMPTY)) {
 825                 isl_basic_set_free(hull);
 826                 return isl_basic_map_set_to_empty(bmap);
 827         }
 828         bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim), 0,
 829                                         hull->n_eq, 0);
 830         for (i = 0; i < hull->n_eq; ++i) {
 831                 j = isl_basic_map_alloc_equality(bmap);
 832                 if (j < 0)
 833                         goto error;
 834                 isl_seq_cpy(bmap->eq[j], hull->eq[i],
 835                                 1 + isl_basic_set_total_dim(hull));
 836         }
 837         isl_vec_free(bmap->sample);
 838         bmap->sample = isl_vec_copy(hull->sample);
 839         isl_basic_set_free(hull);
 840         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT | ISL_BASIC_MAP_ALL_EQUALITIES);
 841         bmap = isl_basic_map_simplify(bmap);
 842         return isl_basic_map_finalize(bmap);
 843 error:
 844         isl_basic_set_free(hull);
 845         isl_basic_map_free(bmap);
 846         return NULL;
 847 }
 848
 849 __isl_give isl_basic_set *isl_basic_set_detect_equalities(
 850                                                 __isl_take isl_basic_set *bset)
 851 {
 852         return (isl_basic_set *)
 853                 isl_basic_map_detect_equalities((isl_basic_map *)bset);
 854 }
 855
 856 struct isl_map *isl_map_detect_equalities(struct isl_map *map)
 857 {
 858         struct isl_basic_map *bmap;
 859         int i;
 860
 861         if (!map)
 862                 return NULL;
 863
 864         for (i = 0; i < map->n; ++i) {
 865                 bmap = isl_basic_map_copy(map->p[i]);
 866                 bmap = isl_basic_map_detect_equalities(bmap);
 867                 if (!bmap)
 868                         goto error;
 869                 isl_basic_map_free(map->p[i]);
 870                 map->p[i] = bmap;
 871         }
 872
 873         return map;
 874 error:
 875         isl_map_free(map);
 876         return NULL;
 877 }
 878
 879 __isl_give isl_set *isl_set_detect_equalities(__isl_take isl_set *set)
 880 {
 881         return (isl_set *)isl_map_detect_equalities((isl_map *)set);
 882 }
 883
 884 /* After computing the rational affine hull (by detecting the implicit
 885  * equalities), we compute the additional equalities satisfied by
 886  * the integer points (if any) and add the original equalities back in.
 887  */
 888 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_basic_map *bmap)
 889 {
 890         bmap = isl_basic_map_detect_equalities(bmap);
 891         bmap = isl_basic_map_cow(bmap);
 892         isl_basic_map_free_inequality(bmap, bmap->n_ineq);
 893         return bmap;
 894 }
 895
 896 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_basic_set *bset)
 897 {
 898         return (struct isl_basic_set *)
 899                 isl_basic_map_affine_hull((struct isl_basic_map *)bset);
 900 }
 901
 902 struct isl_basic_map *isl_map_affine_hull(struct isl_map *map)
 903 {
 904         int i;
 905         struct isl_basic_map *model = NULL;
 906         struct isl_basic_map *hull = NULL;
 907         struct isl_set *set;
 908
 909         if (!map)
 910                 return NULL;
 911
 912         if (map->n == 0) {
 913                 hull = isl_basic_map_empty_like_map(map);
 914                 isl_map_free(map);
 915                 return hull;
 916         }
 917
 918         map = isl_map_detect_equalities(map);
 919         map = isl_map_align_divs(map);
 920         if (!map)
 921                 return NULL;
 922         model = isl_basic_map_copy(map->p[0]);
 923         set = isl_map_underlying_set(map);
 924         set = isl_set_cow(set);
 925         if (!set)
 926                 goto error;
 927
 928         for (i = 0; i < set->n; ++i) {
 929                 set->p[i] = isl_basic_set_cow(set->p[i]);
 930                 set->p[i] = isl_basic_set_affine_hull(set->p[i]);
 931                 set->p[i] = isl_basic_set_gauss(set->p[i], NULL);
 932                 if (!set->p[i])
 933                         goto error;
 934         }
 935         set = isl_set_remove_empty_parts(set);
 936         if (set->n == 0) {
 937                 hull = isl_basic_map_empty_like(model);
 938                 isl_basic_map_free(model);
 939         } else {
 940                 struct isl_basic_set *bset;
 941                 while (set->n > 1) {
 942                         set->p[0] = affine_hull(set->p[0], set->p[--set->n]);
 943                         if (!set->p[0])
 944                                 goto error;
 945                 }
 946                 bset = isl_basic_set_copy(set->p[0]);
 947                 hull = isl_basic_map_overlying_set(bset, model);
 948         }
 949         isl_set_free(set);
 950         hull = isl_basic_map_simplify(hull);
 951         return isl_basic_map_finalize(hull);
 952 error:
 953         isl_basic_map_free(model);
 954         isl_set_free(set);
 955         return NULL;
 956 }
 957
 958 struct isl_basic_set *isl_set_affine_hull(struct isl_set *set)
 959 {
 960         return (struct isl_basic_set *)
 961                 isl_map_affine_hull((struct isl_map *)set);
 962 }