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(struct isl_ctx *ctx,
  11                                                 struct isl_basic_map *bmap)
  12 {
  13         int i;
  14         isl_int opt;
  15
  16         if (!bmap)
  17                 return bmap;
  18
  19         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
  20                 return bmap;
  21         if (F_ISSET(bmap, ISL_BASIC_MAP_NO_IMPLICIT))
  22                 return bmap;
  23
  24         isl_int_init(opt);
  25         for (i = 0; i < bmap->n_ineq; ++i) {
  26                 enum isl_lp_result res;
  27                 res = isl_solve_lp(bmap, 1, bmap->ineq[i]+1, ctx->one, &opt, NULL);
  28                 if (res == isl_lp_unbounded)
  29                         continue;
  30                 if (res == isl_lp_error)
  31                         goto error;
  32                 if (res == isl_lp_empty) {
  33                         bmap = isl_basic_map_set_to_empty(ctx, bmap);
  34                         break;
  35                 }
  36                 isl_int_add(opt, opt, bmap->ineq[i][0]);
  37                 if (isl_int_is_zero(opt)) {
  38                         isl_basic_map_inequality_to_equality(ctx, bmap, i);
  39                         --i;
  40                 }
  41         }
  42         isl_int_clear(opt);
  43
  44         F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
  45         return bmap;
  46 error:
  47         isl_int_clear(opt);
  48         isl_basic_map_free(ctx, bmap);
  49         return NULL;
  50 }
  51
  52 /* Make eq[row][col] of both bmaps equal so we can add the row
  53  * add the column to the common matrix.
  54  * Note that because of the echelon form, the columns of row row
  55  * after column col are zero.
  56  */
  57 static void set_common_multiple(
  58         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
  59         unsigned row, unsigned col)
  60 {
  61         isl_int m, c;
  62
  63         if (isl_int_eq(bset1->eq[row][col], bset2->eq[row][col]))
  64                 return;
  65
  66         isl_int_init(c);
  67         isl_int_init(m);
  68         isl_int_lcm(m, bset1->eq[row][col], bset2->eq[row][col]);
  69         isl_int_divexact(c, m, bset1->eq[row][col]);
  70         isl_seq_scale(bset1->eq[row], bset1->eq[row], c, col+1);
  71         isl_int_divexact(c, m, bset2->eq[row][col]);
  72         isl_seq_scale(bset2->eq[row], bset2->eq[row], c, col+1);
  73         isl_int_clear(c);
  74         isl_int_clear(m);
  75 }
  76
  77 /* Delete a given equality, moving all the following equalities one up.
  78  */
  79 static void delete_row(struct isl_basic_set *bset, unsigned row)
  80 {
  81         isl_int *t;
  82         int r;
  83
  84         t = bset->eq[row];
  85         bset->n_eq--;
  86         for (r = row; r < bset->n_eq; ++r)
  87                 bset->eq[r] = bset->eq[r+1];
  88         bset->eq[bset->n_eq] = t;
  89 }
  90
  91 /* Make first row entries in column col of bset1 identical to
  92  * those of bset2, using the fact that entry bset1->eq[row][col]=a
  93  * is non-zero.  Initially, these elements of bset1 are all zero.
  94  * For each row i < row, we set
  95  *              A[i] = a * A[i] + B[i][col] * A[row]
  96  *              B[i] = a * B[i]
  97  * so that
  98  *              A[i][col] = B[i][col] = a * old(B[i][col])
  99  */
 100 static void construct_column(
 101         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 102         unsigned row, unsigned col)
 103 {
 104         int r;
 105         isl_int a;
 106         isl_int b;
 107         unsigned total;
 108
 109         isl_int_init(a);
 110         isl_int_init(b);
 111         total = 1 + bset1->dim;
 112         for (r = 0; r < row; ++r) {
 113                 if (isl_int_is_zero(bset2->eq[r][col]))
 114                         continue;
 115                 isl_int_gcd(b, bset2->eq[r][col], bset1->eq[row][col]);
 116                 isl_int_divexact(a, bset1->eq[row][col], b);
 117                 isl_int_divexact(b, bset2->eq[r][col], b);
 118                 isl_seq_combine(bset1->eq[r], a, bset1->eq[r],
 119                                               b, bset1->eq[row], total);
 120                 isl_seq_scale(bset2->eq[r], bset2->eq[r], a, total);
 121         }
 122         isl_int_clear(a);
 123         isl_int_clear(b);
 124         delete_row(bset1, row);
 125 }
 126
 127 /* Make first row entries in column col of bset1 identical to
 128  * those of bset2, using only these entries of the two matrices.
 129  * Let t be the last row with different entries.
 130  * For each row i < t, we set
 131  *      A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
 132  *      B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
 133  * so that
 134  *      A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
 135  */
 136 static int transform_column(
 137         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 138         unsigned row, unsigned col)
 139 {
 140         int i, t;
 141         isl_int a, b, g;
 142         unsigned total;
 143
 144         for (t = row-1; t >= 0; --t)
 145                 if (isl_int_ne(bset1->eq[t][col], bset2->eq[t][col]))
 146                         break;
 147         if (t < 0)
 148                 return 0;
 149
 150         total = 1 + bset1->dim;
 151         isl_int_init(a);
 152         isl_int_init(b);
 153         isl_int_init(g);
 154         isl_int_sub(b, bset1->eq[t][col], bset2->eq[t][col]);
 155         for (i = 0; i < t; ++i) {
 156                 isl_int_sub(a, bset2->eq[i][col], bset1->eq[i][col]);
 157                 isl_int_gcd(g, a, b);
 158                 isl_int_divexact(a, a, g);
 159                 isl_int_divexact(g, b, g);
 160                 isl_seq_combine(bset1->eq[i], g, bset1->eq[i], a, bset1->eq[t],
 161                                 total);
 162                 isl_seq_combine(bset2->eq[i], g, bset2->eq[i], a, bset2->eq[t],
 163                                 total);
 164         }
 165         isl_int_clear(a);
 166         isl_int_clear(b);
 167         isl_int_clear(g);
 168         delete_row(bset1, t);
 169         delete_row(bset2, t);
 170         return 1;
 171 }
 172
 173 /* The implementation is based on Section 5.2 of Michael Karr,
 174  * "Affine Relationships Among Variables of a Program",
 175  * except that the echelon form we use starts from the last column
 176  * and that we are dealing with integer coefficients.
 177  */
 178 static struct isl_basic_set *affine_hull(struct isl_ctx *ctx,
 179         struct isl_basic_set *bset1, struct isl_basic_set *bset2)
 180 {
 181         unsigned total;
 182         int col;
 183         int row;
 184
 185         total = 1 + bset1->dim;
 186
 187         row = 0;
 188         for (col = total-1; col >= 0; --col) {
 189                 int is_zero1 = row >= bset1->n_eq ||
 190                         isl_int_is_zero(bset1->eq[row][col]);
 191                 int is_zero2 = row >= bset2->n_eq ||
 192                         isl_int_is_zero(bset2->eq[row][col]);
 193                 if (!is_zero1 && !is_zero2) {
 194                         set_common_multiple(bset1, bset2, row, col);
 195                         ++row;
 196                 } else if (!is_zero1 && is_zero2) {
 197                         construct_column(bset1, bset2, row, col);
 198                 } else if (is_zero1 && !is_zero2) {
 199                         construct_column(bset2, bset1, row, col);
 200                 } else {
 201                         if (transform_column(bset1, bset2, row, col))
 202                                 --row;
 203                 }
 204         }
 205         isl_basic_set_free(ctx, bset2);
 206         isl_assert(ctx, row == bset1->n_eq, goto error);
 207         return bset1;
 208 error:
 209         isl_basic_set_free(ctx, bset1);
 210         return NULL;
 211 }
 212
 213 static struct isl_basic_set *isl_basic_set_from_vec(struct isl_ctx *ctx,
 214         struct isl_vec *vec)
 215 {
 216         int i;
 217         int k;
 218         struct isl_basic_set *bset = NULL;
 219
 220         if (!vec)
 221                 return NULL;
 222         isl_assert(ctx, vec->size != 0, goto error);
 223
 224         bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
 225         for (i = bset->dim - 1; i >= 0; --i) {
 226                 k = isl_basic_set_alloc_equality(ctx, bset);
 227                 if (k < 0)
 228                         goto error;
 229                 isl_seq_clr(bset->eq[k], 1 + bset->dim);
 230                 isl_int_neg(bset->eq[k][0], vec->block.data[1 + i]);
 231                 isl_int_set(bset->eq[k][1 + i], vec->block.data[0]);
 232         }
 233         isl_vec_free(ctx, vec);
 234
 235         return bset;
 236 error:
 237         isl_basic_set_free(ctx, bset);
 238         isl_vec_free(ctx, vec);
 239         return NULL;
 240 }
 241
 242 static struct isl_basic_set *outside_point(struct isl_ctx *ctx,
 243         struct isl_basic_set *bset, isl_int *eq, int up)
 244 {
 245         struct isl_basic_set *slice = NULL;
 246         struct isl_vec *sample;
 247         struct isl_basic_set *point;
 248         int k;
 249
 250         slice = isl_basic_set_copy(ctx, bset);
 251         if (!slice)
 252                 goto error;
 253         slice = isl_basic_set_extend(ctx, slice, 0, slice->dim, 0, 1, 0);
 254         k = isl_basic_set_alloc_equality(ctx, slice);
 255         if (k < 0)
 256                 goto error;
 257         isl_seq_cpy(slice->eq[k], eq, 1 + slice->dim);
 258         if (up)
 259                 isl_int_add_ui(slice->eq[k][0], slice->eq[k][0], 1);
 260         else
 261                 isl_int_sub_ui(slice->eq[k][0], slice->eq[k][0], 1);
 262
 263         sample = isl_basic_set_sample(ctx, slice);
 264         if (!sample)
 265                 goto error;
 266         if (sample->size == 0) {
 267                 isl_vec_free(ctx, sample);
 268                 point = isl_basic_set_empty(ctx, 0, bset->dim);
 269         } else
 270                 point = isl_basic_set_from_vec(ctx, sample);
 271
 272         return point;
 273 error:
 274         isl_basic_set_free(ctx, slice);
 275         return NULL;
 276 }
 277
 278 /* After computing the rational affine hull (by detecting the implicit
 279  * equalities), we remove all equalities found so far, compute
 280  * the integer affine hull of what is left, and then add the original
 281  * equalities back in.
 282  *
 283  * The integer affine hull is constructed by successively looking
 284  * a point that is affinely independent of the points found so far.
 285  * In particular, for each equality satisfied by the points so far,
 286  * we check if there is any point on the corresponding hyperplane
 287  * shifted by one (in either direction).
 288  */
 289 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_ctx *ctx,
 290                                                 struct isl_basic_map *bmap)
 291 {
 292         int i, j;
 293         struct isl_mat *T2 = NULL;
 294         struct isl_basic_set *bset = NULL;
 295         struct isl_basic_set *hull = NULL;
 296         struct isl_vec *sample;
 297
 298         bmap = isl_basic_map_implicit_equalities(ctx, bmap);
 299         if (!bmap)
 300                 return NULL;
 301         if (bmap->n_ineq == 0)
 302                 return bmap;
 303
 304         bset = isl_basic_map_underlying_set(ctx, isl_basic_map_copy(ctx, bmap));
 305         bset = isl_basic_set_remove_equalities(ctx, bset, NULL, &T2);
 306
 307         sample = isl_basic_set_sample(ctx, isl_basic_set_copy(ctx, bset));
 308         hull = isl_basic_set_from_vec(ctx, sample);
 309
 310         for (i = 0; i < bset->dim; ++i) {
 311                 struct isl_basic_set *point;
 312                 for (j = 0; j < hull->n_eq; ++j) {
 313                         point = outside_point(ctx, bset, hull->eq[j], 1);
 314                         if (!point)
 315                                 goto error;
 316                         if (!F_ISSET(point, ISL_BASIC_SET_EMPTY))
 317                                 break;
 318                         isl_basic_set_free(ctx, point);
 319                         point = outside_point(ctx, bset, hull->eq[j], 0);
 320                         if (!point)
 321                                 goto error;
 322                         if (!F_ISSET(point, ISL_BASIC_SET_EMPTY))
 323                                 break;
 324                         isl_basic_set_free(ctx, point);
 325                 }
 326                 if (j == hull->n_eq)
 327                         break;
 328                 hull = affine_hull(ctx, hull, point);
 329         }
 330
 331         isl_basic_set_free(ctx, bset);
 332         bset = NULL;
 333         bmap = isl_basic_map_cow(ctx, bmap);
 334         if (!bmap)
 335                 goto error;
 336         isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
 337         if (T2)
 338                 hull = isl_basic_set_preimage(ctx, hull, T2);
 339         bmap = isl_basic_map_intersect(ctx, bmap,
 340                                         isl_basic_map_overlying_set(ctx, hull,
 341                                                 isl_basic_map_copy(ctx, bmap)));
 342
 343         return isl_basic_map_finalize(ctx, bmap);
 344 error:
 345         isl_mat_free(ctx, T2);
 346         isl_basic_set_free(ctx, bset);
 347         isl_basic_set_free(ctx, hull);
 348         isl_basic_map_free(ctx, bmap);
 349         return NULL;
 350 }
 351
 352 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
 353                                                 struct isl_basic_set *bset)
 354 {
 355         return (struct isl_basic_set *)
 356                 isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
 357 }
 358
 359 struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
 360                                                 struct isl_map *map)
 361 {
 362         int i;
 363         struct isl_basic_map *model = NULL;
 364         struct isl_basic_map *hull = NULL;
 365         struct isl_set *set;
 366
 367         if (!map)
 368                 return NULL;
 369
 370         if (map->n == 0) {
 371                 hull = isl_basic_map_empty(ctx,
 372                                             map->nparam, map->n_in, map->n_out);
 373                 isl_map_free(ctx, map);
 374                 return hull;
 375         }
 376
 377         map = isl_map_align_divs(ctx, map);
 378         model = isl_basic_map_copy(ctx, map->p[0]);
 379         set = isl_map_underlying_set(ctx, map);
 380         set = isl_set_cow(ctx, set);
 381         if (!set)
 382                 goto error;
 383
 384         for (i = 0; i < set->n; ++i) {
 385                 set->p[i] = isl_basic_set_cow(ctx, set->p[i]);
 386                 set->p[i] = isl_basic_set_affine_hull(ctx, set->p[i]);
 387                 set->p[i] = isl_basic_set_gauss(ctx, set->p[i], NULL);
 388                 if (!set->p[i])
 389                         goto error;
 390         }
 391         set = isl_set_remove_empty_parts(ctx, set);
 392         if (set->n == 0) {
 393                 hull = isl_basic_map_empty(ctx,
 394                                     model->nparam, model->n_in, model->n_out);
 395                 isl_basic_map_free(ctx, model);
 396         } else {
 397                 struct isl_basic_set *bset;
 398                 while (set->n > 1) {
 399                         set->p[0] = affine_hull(ctx, set->p[0],
 400                                                         set->p[--set->n]);
 401                         if (!set->p[0])
 402                                 goto error;
 403                 }
 404                 bset = isl_basic_set_copy(ctx, set->p[0]);
 405                 hull = isl_basic_map_overlying_set(ctx, bset, model);
 406         }
 407         isl_set_free(ctx, set);
 408         hull = isl_basic_map_simplify(ctx, hull);
 409         return isl_basic_map_finalize(ctx, hull);
 410 error:
 411         isl_basic_map_free(ctx, model);
 412         isl_set_free(ctx, set);
 413         return NULL;
 414 }
 415
 416 struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
 417                                                 struct isl_set *set)
 418 {
 419         return (struct isl_basic_set *)
 420                 isl_map_affine_hull(ctx, (struct isl_map *)set);
 421 }