isl_affine_hull.c

   1 #include "isl_ctx.h"
   2 #include "isl_seq.h"
   3 #include "isl_set.h"
   4 #include "isl_lp.h"
   5 #include "isl_map.h"
   6 #include "isl_map_private.h"
   7 #include "isl_equalities.h"
   8 #include "isl_sample.h"
   9
  10 struct isl_basic_map *isl_basic_map_implicit_equalities(
  11                                                 struct isl_basic_map *bmap)
  12 {
  13         int i;
  14         int rational;
  15         isl_int opt;
  16         isl_int opt_denom;
  17         struct isl_ctx *ctx;
  18
  19         if (!bmap)
  20                 return bmap;
  21
  22         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
  23                 return bmap;
  24         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_IMPLICIT))
  25                 return bmap;
  26
  27         ctx = bmap->ctx;
  28         rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
  29         isl_int_init(opt);
  30         isl_int_init(opt_denom);
  31         if (!rational)
  32                 isl_int_set_si(opt_denom, 1);
  33         for (i = 0; i < bmap->n_ineq; ++i) {
  34                 enum isl_lp_result res;
  35                 res = isl_solve_lp(bmap, 1, bmap->ineq[i]+1, ctx->one,
  36                                         &opt, rational ? &opt_denom : NULL);
  37                 if (res == isl_lp_unbounded)
  38                         continue;
  39                 if (res == isl_lp_error)
  40                         goto error;
  41                 if (res == isl_lp_empty) {
  42                         bmap = isl_basic_map_set_to_empty(bmap);
  43                         break;
  44                 }
  45                 if (!isl_int_is_one(opt_denom))
  46                         continue;
  47                 isl_int_add(opt, opt, bmap->ineq[i][0]);
  48                 if (isl_int_is_zero(opt)) {
  49                         isl_basic_map_inequality_to_equality(bmap, i);
  50                         --i;
  51                 }
  52         }
  53         isl_int_clear(opt_denom);
  54         isl_int_clear(opt);
  55
  56         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
  57         return bmap;
  58 error:
  59         isl_int_clear(opt);
  60         isl_basic_map_free(bmap);
  61         return NULL;
  62 }
  63
  64 /* Make eq[row][col] of both bmaps equal so we can add the row
  65  * add the column to the common matrix.
  66  * Note that because of the echelon form, the columns of row row
  67  * after column col are zero.
  68  */
  69 static void set_common_multiple(
  70         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
  71         unsigned row, unsigned col)
  72 {
  73         isl_int m, c;
  74
  75         if (isl_int_eq(bset1->eq[row][col], bset2->eq[row][col]))
  76                 return;
  77
  78         isl_int_init(c);
  79         isl_int_init(m);
  80         isl_int_lcm(m, bset1->eq[row][col], bset2->eq[row][col]);
  81         isl_int_divexact(c, m, bset1->eq[row][col]);
  82         isl_seq_scale(bset1->eq[row], bset1->eq[row], c, col+1);
  83         isl_int_divexact(c, m, bset2->eq[row][col]);
  84         isl_seq_scale(bset2->eq[row], bset2->eq[row], c, col+1);
  85         isl_int_clear(c);
  86         isl_int_clear(m);
  87 }
  88
  89 /* Delete a given equality, moving all the following equalities one up.
  90  */
  91 static void delete_row(struct isl_basic_set *bset, unsigned row)
  92 {
  93         isl_int *t;
  94         int r;
  95
  96         t = bset->eq[row];
  97         bset->n_eq--;
  98         for (r = row; r < bset->n_eq; ++r)
  99                 bset->eq[r] = bset->eq[r+1];
 100         bset->eq[bset->n_eq] = t;
 101 }
 102
 103 /* Make first row entries in column col of bset1 identical to
 104  * those of bset2, using the fact that entry bset1->eq[row][col]=a
 105  * is non-zero.  Initially, these elements of bset1 are all zero.
 106  * For each row i < row, we set
 107  *              A[i] = a * A[i] + B[i][col] * A[row]
 108  *              B[i] = a * B[i]
 109  * so that
 110  *              A[i][col] = B[i][col] = a * old(B[i][col])
 111  */
 112 static void construct_column(
 113         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 114         unsigned row, unsigned col)
 115 {
 116         int r;
 117         isl_int a;
 118         isl_int b;
 119         unsigned total;
 120
 121         isl_int_init(a);
 122         isl_int_init(b);
 123         total = 1 + isl_basic_set_n_dim(bset1);
 124         for (r = 0; r < row; ++r) {
 125                 if (isl_int_is_zero(bset2->eq[r][col]))
 126                         continue;
 127                 isl_int_gcd(b, bset2->eq[r][col], bset1->eq[row][col]);
 128                 isl_int_divexact(a, bset1->eq[row][col], b);
 129                 isl_int_divexact(b, bset2->eq[r][col], b);
 130                 isl_seq_combine(bset1->eq[r], a, bset1->eq[r],
 131                                               b, bset1->eq[row], total);
 132                 isl_seq_scale(bset2->eq[r], bset2->eq[r], a, total);
 133         }
 134         isl_int_clear(a);
 135         isl_int_clear(b);
 136         delete_row(bset1, row);
 137 }
 138
 139 /* Make first row entries in column col of bset1 identical to
 140  * those of bset2, using only these entries of the two matrices.
 141  * Let t be the last row with different entries.
 142  * For each row i < t, we set
 143  *      A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
 144  *      B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
 145  * so that
 146  *      A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
 147  */
 148 static int transform_column(
 149         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 150         unsigned row, unsigned col)
 151 {
 152         int i, t;
 153         isl_int a, b, g;
 154         unsigned total;
 155
 156         for (t = row-1; t >= 0; --t)
 157                 if (isl_int_ne(bset1->eq[t][col], bset2->eq[t][col]))
 158                         break;
 159         if (t < 0)
 160                 return 0;
 161
 162         total = 1 + isl_basic_set_n_dim(bset1);
 163         isl_int_init(a);
 164         isl_int_init(b);
 165         isl_int_init(g);
 166         isl_int_sub(b, bset1->eq[t][col], bset2->eq[t][col]);
 167         for (i = 0; i < t; ++i) {
 168                 isl_int_sub(a, bset2->eq[i][col], bset1->eq[i][col]);
 169                 isl_int_gcd(g, a, b);
 170                 isl_int_divexact(a, a, g);
 171                 isl_int_divexact(g, b, g);
 172                 isl_seq_combine(bset1->eq[i], g, bset1->eq[i], a, bset1->eq[t],
 173                                 total);
 174                 isl_seq_combine(bset2->eq[i], g, bset2->eq[i], a, bset2->eq[t],
 175                                 total);
 176         }
 177         isl_int_clear(a);
 178         isl_int_clear(b);
 179         isl_int_clear(g);
 180         delete_row(bset1, t);
 181         delete_row(bset2, t);
 182         return 1;
 183 }
 184
 185 /* The implementation is based on Section 5.2 of Michael Karr,
 186  * "Affine Relationships Among Variables of a Program",
 187  * except that the echelon form we use starts from the last column
 188  * and that we are dealing with integer coefficients.
 189  */
 190 static struct isl_basic_set *affine_hull(
 191         struct isl_basic_set *bset1, struct isl_basic_set *bset2)
 192 {
 193         unsigned total;
 194         int col;
 195         int row;
 196
 197         total = 1 + isl_basic_set_n_dim(bset1);
 198
 199         row = 0;
 200         for (col = total-1; col >= 0; --col) {
 201                 int is_zero1 = row >= bset1->n_eq ||
 202                         isl_int_is_zero(bset1->eq[row][col]);
 203                 int is_zero2 = row >= bset2->n_eq ||
 204                         isl_int_is_zero(bset2->eq[row][col]);
 205                 if (!is_zero1 && !is_zero2) {
 206                         set_common_multiple(bset1, bset2, row, col);
 207                         ++row;
 208                 } else if (!is_zero1 && is_zero2) {
 209                         construct_column(bset1, bset2, row, col);
 210                 } else if (is_zero1 && !is_zero2) {
 211                         construct_column(bset2, bset1, row, col);
 212                 } else {
 213                         if (transform_column(bset1, bset2, row, col))
 214                                 --row;
 215                 }
 216         }
 217         isl_basic_set_free(bset2);
 218         isl_assert(ctx, row == bset1->n_eq, goto error);
 219         return bset1;
 220 error:
 221         isl_basic_set_free(bset1);
 222         return NULL;
 223 }
 224
 225 static struct isl_basic_set *isl_basic_set_from_vec(struct isl_ctx *ctx,
 226         struct isl_vec *vec)
 227 {
 228         int i;
 229         int k;
 230         struct isl_basic_set *bset = NULL;
 231         unsigned dim;
 232
 233         if (!vec)
 234                 return NULL;
 235         isl_assert(ctx, vec->size != 0, goto error);
 236
 237         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
 238         if (!bset)
 239                 goto error;
 240         dim = isl_basic_set_n_dim(bset);
 241         for (i = dim - 1; i >= 0; --i) {
 242                 k = isl_basic_set_alloc_equality(bset);
 243                 if (k < 0)
 244                         goto error;
 245                 isl_seq_clr(bset->eq[k], 1 + dim);
 246                 isl_int_neg(bset->eq[k][0], vec->block.data[1 + i]);
 247                 isl_int_set(bset->eq[k][1 + i], vec->block.data[0]);
 248         }
 249         isl_vec_free(ctx, vec);
 250
 251         return bset;
 252 error:
 253         isl_basic_set_free(bset);
 254         isl_vec_free(ctx, vec);
 255         return NULL;
 256 }
 257
 258 /* Find an integer point in "bset" that lies outside of the equality
 259  * "eq" e(x) = 0.
 260  * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
 261  * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
 262  * The point, if found, is returned as a singleton set.
 263  * If no point can be found, the empty set is returned.
 264  */
 265 static struct isl_basic_set *outside_point(struct isl_ctx *ctx,
 266         struct isl_basic_set *bset, isl_int *eq, int up)
 267 {
 268         struct isl_basic_set *slice = NULL;
 269         struct isl_vec *sample;
 270         struct isl_basic_set *point;
 271         unsigned dim;
 272         int k;
 273
 274         slice = isl_basic_set_copy(bset);
 275         if (!slice)
 276                 goto error;
 277         dim = isl_basic_set_n_dim(slice);
 278         slice = isl_basic_set_extend(slice, 0, dim, 0, 0, 1);
 279         k = isl_basic_set_alloc_inequality(slice);
 280         if (k < 0)
 281                 goto error;
 282         if (up)
 283                 isl_seq_cpy(slice->ineq[k], eq, 1 + dim);
 284         else
 285                 isl_seq_neg(slice->ineq[k], eq, 1 + dim);
 286         isl_int_sub_ui(slice->ineq[k][0], slice->ineq[k][0], 1);
 287
 288         sample = isl_basic_set_sample(slice);
 289         if (!sample)
 290                 goto error;
 291         if (sample->size == 0) {
 292                 isl_vec_free(ctx, sample);
 293                 point = isl_basic_set_empty_like(bset);
 294         } else
 295                 point = isl_basic_set_from_vec(ctx, sample);
 296
 297         return point;
 298 error:
 299         isl_basic_set_free(slice);
 300         return NULL;
 301 }
 302
 303 /* Look for all equalities satisfied by the integer points in bmap
 304  * that are independent of the equalities already explicitly available
 305  * in bmap.
 306  *
 307  * We first remove all equalities already explicitly available,
 308  * then look for additional equalities in the reduced space
 309  * and then transform the result to the original space.
 310  * The original equalities are _not_ added to this set.  This is
 311  * the responsibility of the calling function.
 312  * The resulting basic set has all meaning about the dimensions removed.
 313  * In particular, dimensions that correspond to existential variables
 314  * in bmap and that are found to be fixed are not removed.
 315  *
 316  * The additional equalities are obtained by successively looking for
 317  * a point that is affinely independent of the points found so far.
 318  * In particular, for each equality satisfied by the points so far,
 319  * we check if there is any point on a hyperplane parallel to the
 320  * corresponding hyperplane shifted by at least one (in either direction).
 321  */
 322 static struct isl_basic_set *equalities_in_underlying_set(
 323                                                 struct isl_basic_map *bmap)
 324 {
 325         int i, j;
 326         struct isl_mat *T2 = NULL;
 327         struct isl_basic_set *bset = NULL;
 328         struct isl_basic_set *hull = NULL;
 329         struct isl_vec *sample;
 330         struct isl_ctx *ctx;
 331         unsigned dim;
 332
 333         bset = isl_basic_map_underlying_set(bmap);
 334         bset = isl_basic_set_remove_equalities(bset, NULL, &T2);
 335         if (!bset)
 336                 goto error;
 337
 338         ctx = bset->ctx;
 339         sample = isl_basic_set_sample(isl_basic_set_copy(bset));
 340         if (!sample)
 341                 goto error;
 342         if (sample->size == 0) {
 343                 isl_vec_free(ctx, sample);
 344                 hull = isl_basic_set_empty_like(bset);
 345         } else
 346                 hull = isl_basic_set_from_vec(ctx, sample);
 347
 348         dim = isl_basic_set_n_dim(bset);
 349         for (i = 0; i < dim; ++i) {
 350                 struct isl_basic_set *point;
 351                 if (ISL_F_ISSET(hull, ISL_BASIC_SET_EMPTY))
 352                         break;
 353                 for (j = 0; j < hull->n_eq; ++j) {
 354                         point = outside_point(ctx, bset, hull->eq[j], 1);
 355                         if (!point)
 356                                 goto error;
 357                         if (!ISL_F_ISSET(point, ISL_BASIC_SET_EMPTY))
 358                                 break;
 359                         isl_basic_set_free(point);
 360                         point = outside_point(ctx, bset, hull->eq[j], 0);
 361                         if (!point)
 362                                 goto error;
 363                         if (!ISL_F_ISSET(point, ISL_BASIC_SET_EMPTY))
 364                                 break;
 365                         isl_basic_set_free(point);
 366                 }
 367                 if (j == hull->n_eq)
 368                         break;
 369                 hull = affine_hull(hull, point);
 370         }
 371         isl_basic_set_free(bset);
 372         if (T2)
 373                 hull = isl_basic_set_preimage(ctx, hull, T2);
 374
 375         return hull;
 376 error:
 377         isl_mat_free(ctx, T2);
 378         isl_basic_set_free(bset);
 379         isl_basic_set_free(hull);
 380         return NULL;
 381 }
 382
 383 /* Detect and make explicit all equalities satisfied by the (integer)
 384  * points in bmap.
 385  */
 386 struct isl_basic_map *isl_basic_map_detect_equalities(
 387                                                 struct isl_basic_map *bmap)
 388 {
 389         int i, j;
 390         struct isl_basic_set *hull = NULL;
 391
 392         if (!bmap)
 393                 return NULL;
 394         if (bmap->n_ineq == 0)
 395                 return bmap;
 396         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
 397                 return bmap;
 398         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_ALL_EQUALITIES))
 399                 return bmap;
 400         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
 401                 return isl_basic_map_implicit_equalities(bmap);
 402
 403         hull = equalities_in_underlying_set(isl_basic_map_copy(bmap));
 404         if (!hull)
 405                 goto error;
 406         bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim), 0,
 407                                         hull->n_eq, 0);
 408         for (i = 0; i < hull->n_eq; ++i) {
 409                 j = isl_basic_map_alloc_equality(bmap);
 410                 if (j < 0)
 411                         goto error;
 412                 isl_seq_cpy(bmap->eq[j], hull->eq[i],
 413                                 1 + isl_basic_set_total_dim(hull));
 414         }
 415         isl_basic_set_free(hull);
 416         ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT | ISL_BASIC_MAP_ALL_EQUALITIES);
 417         bmap = isl_basic_map_simplify(bmap);
 418         return isl_basic_map_finalize(bmap);
 419 error:
 420         isl_basic_set_free(hull);
 421         isl_basic_map_free(bmap);
 422         return NULL;
 423 }
 424
 425 /* After computing the rational affine hull (by detecting the implicit
 426  * equalities), we compute the additional equalities satisfied by
 427  * the integer points (if any) and add the original equalities back in.
 428  */
 429 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_basic_map *bmap)
 430 {
 431         struct isl_basic_set *hull = NULL;
 432
 433         bmap = isl_basic_map_implicit_equalities(bmap);
 434         if (!bmap)
 435                 return NULL;
 436         if (bmap->n_ineq == 0)
 437                 return bmap;
 438
 439         if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
 440                 bmap = isl_basic_map_cow(bmap);
 441                 isl_basic_map_free_inequality(bmap, bmap->n_ineq);
 442                 return bmap;
 443         }
 444
 445         hull = equalities_in_underlying_set(isl_basic_map_copy(bmap));
 446         if (!hull)
 447                 goto error;
 448
 449         bmap = isl_basic_map_cow(bmap);
 450         if (!bmap)
 451                 goto error;
 452         isl_basic_map_free_inequality(bmap, bmap->n_ineq);
 453         bmap = isl_basic_map_intersect(bmap,
 454                                         isl_basic_map_overlying_set(hull,
 455                                                 isl_basic_map_copy(bmap)));
 456
 457         return isl_basic_map_finalize(bmap);
 458 error:
 459         isl_basic_set_free(hull);
 460         isl_basic_map_free(bmap);
 461         return NULL;
 462 }
 463
 464 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_basic_set *bset)
 465 {
 466         return (struct isl_basic_set *)
 467                 isl_basic_map_affine_hull((struct isl_basic_map *)bset);
 468 }
 469
 470 struct isl_basic_map *isl_map_affine_hull(struct isl_map *map)
 471 {
 472         int i;
 473         struct isl_basic_map *model = NULL;
 474         struct isl_basic_map *hull = NULL;
 475         struct isl_set *set;
 476
 477         if (!map)
 478                 return NULL;
 479
 480         if (map->n == 0) {
 481                 hull = isl_basic_map_empty_like_map(map);
 482                 isl_map_free(map);
 483                 return hull;
 484         }
 485
 486         map = isl_map_align_divs(map);
 487         model = isl_basic_map_copy(map->p[0]);
 488         set = isl_map_underlying_set(map);
 489         set = isl_set_cow(set);
 490         if (!set)
 491                 goto error;
 492
 493         for (i = 0; i < set->n; ++i) {
 494                 set->p[i] = isl_basic_set_cow(set->p[i]);
 495                 set->p[i] = isl_basic_set_affine_hull(set->p[i]);
 496                 set->p[i] = isl_basic_set_gauss(set->p[i], NULL);
 497                 if (!set->p[i])
 498                         goto error;
 499         }
 500         set = isl_set_remove_empty_parts(set);
 501         if (set->n == 0) {
 502                 hull = isl_basic_map_empty_like(model);
 503                 isl_basic_map_free(model);
 504         } else {
 505                 struct isl_basic_set *bset;
 506                 while (set->n > 1) {
 507                         set->p[0] = affine_hull(set->p[0], set->p[--set->n]);
 508                         if (!set->p[0])
 509                                 goto error;
 510                 }
 511                 bset = isl_basic_set_copy(set->p[0]);
 512                 hull = isl_basic_map_overlying_set(bset, model);
 513         }
 514         isl_set_free(set);
 515         hull = isl_basic_map_simplify(hull);
 516         return isl_basic_map_finalize(hull);
 517 error:
 518         isl_basic_map_free(model);
 519         isl_set_free(set);
 520         return NULL;
 521 }
 522
 523 struct isl_basic_set *isl_set_affine_hull(struct isl_set *set)
 524 {
 525         return (struct isl_basic_set *)
 526                 isl_map_affine_hull((struct isl_map *)set);
 527 }