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