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
   8 struct isl_basic_map *isl_basic_map_implicit_equalities(struct isl_ctx *ctx,
   9                                                 struct isl_basic_map *bmap)
  10 {
  11         int i;
  12         isl_int opt;
  13
  14         if (!bmap)
  15                 return bmap;
  16
  17         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
  18                 return bmap;
  19         if (F_ISSET(bmap, ISL_BASIC_MAP_NO_IMPLICIT))
  20                 return bmap;
  21
  22         isl_int_init(opt);
  23         for (i = 0; i < bmap->n_ineq; ++i) {
  24                 enum isl_lp_result res;
  25                 res = isl_solve_lp(bmap, 1, bmap->ineq[i]+1, ctx->one, &opt);
  26                 if (res == isl_lp_unbounded)
  27                         continue;
  28                 if (res == isl_lp_error)
  29                         goto error;
  30                 if (res == isl_lp_empty) {
  31                         bmap = isl_basic_map_set_to_empty(ctx, bmap);
  32                         break;
  33                 }
  34                 isl_int_add(opt, opt, bmap->ineq[i][0]);
  35                 if (isl_int_is_zero(opt)) {
  36                         isl_basic_map_inequality_to_equality(ctx, bmap, i);
  37                         --i;
  38                 }
  39         }
  40         isl_int_clear(opt);
  41
  42         F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
  43         return bmap;
  44 error:
  45         isl_int_clear(opt);
  46         isl_basic_map_free(ctx, bmap);
  47         return NULL;
  48 }
  49
  50 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_ctx *ctx,
  51                                                 struct isl_basic_map *bmap)
  52 {
  53         bmap = isl_basic_map_implicit_equalities(ctx, bmap);
  54         if (!bmap)
  55                 return NULL;
  56         if (bmap->n_ineq == 0)
  57                 return bmap;
  58         bmap = isl_basic_map_cow(ctx, bmap);
  59         if (!bmap)
  60                 return NULL;
  61         isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
  62         return isl_basic_map_finalize(ctx, bmap);
  63 }
  64
  65 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
  66                                                 struct isl_basic_set *bset)
  67 {
  68         return (struct isl_basic_set *)
  69                 isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
  70 }
  71
  72 /* Make eq[row][col] of both bmaps equal so we can add the row
  73  * add the column to the common matrix.
  74  * Note that because of the echelon form, the columns of row row
  75  * after column col are zero.
  76  */
  77 static void set_common_multiple(
  78         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
  79         unsigned row, unsigned col)
  80 {
  81         isl_int m, c;
  82
  83         if (isl_int_eq(bset1->eq[row][col], bset2->eq[row][col]))
  84                 return;
  85
  86         isl_int_init(c);
  87         isl_int_init(m);
  88         isl_int_lcm(m, bset1->eq[row][col], bset2->eq[row][col]);
  89         isl_int_divexact(c, m, bset1->eq[row][col]);
  90         isl_seq_scale(bset1->eq[row], bset1->eq[row], c, col+1);
  91         isl_int_divexact(c, m, bset2->eq[row][col]);
  92         isl_seq_scale(bset2->eq[row], bset2->eq[row], c, col+1);
  93         isl_int_clear(c);
  94         isl_int_clear(m);
  95 }
  96
  97 /* Delete a given equality, moving all the following equalities one up.
  98  */
  99 static void delete_row(struct isl_basic_set *bset, unsigned row)
 100 {
 101         isl_int *t;
 102         int r;
 103
 104         t = bset->eq[row];
 105         bset->n_eq--;
 106         for (r = row; r < bset->n_eq; ++r)
 107                 bset->eq[r] = bset->eq[r+1];
 108         bset->eq[bset->n_eq] = t;
 109 }
 110
 111 /* Make first row entries in column col of bset1 identical to
 112  * those of bset2, using the fact that entry bset1->eq[row][col]=a
 113  * is non-zero.  Initially, these elements of bset1 are all zero.
 114  * For each row i < row, we set
 115  *              A[i] = a * A[i] + B[i][col] * A[row]
 116  *              B[i] = a * B[i]
 117  * so that
 118  *              A[i][col] = B[i][col] = a * old(B[i][col])
 119  */
 120 static void construct_column(
 121         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 122         unsigned row, unsigned col)
 123 {
 124         int r;
 125         isl_int a;
 126         isl_int b;
 127         unsigned total;
 128
 129         isl_int_init(a);
 130         isl_int_init(b);
 131         total = 1 + bset1->dim;
 132         for (r = 0; r < row; ++r) {
 133                 if (isl_int_is_zero(bset2->eq[r][col]))
 134                         continue;
 135                 isl_int_gcd(b, bset2->eq[r][col], bset1->eq[row][col]);
 136                 isl_int_divexact(a, bset1->eq[row][col], b);
 137                 isl_int_divexact(b, bset2->eq[r][col], b);
 138                 isl_seq_combine(bset1->eq[r], a, bset1->eq[r],
 139                                               b, bset1->eq[row], total);
 140                 isl_seq_scale(bset2->eq[r], bset2->eq[r], a, total);
 141         }
 142         isl_int_clear(a);
 143         isl_int_clear(b);
 144         delete_row(bset1, row);
 145 }
 146
 147 /* Make first row entries in column col of bset1 identical to
 148  * those of bset2, using only these entries of the two matrices.
 149  * Let t be the last row with different entries.
 150  * For each row i < t, we set
 151  *      A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
 152  *      B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
 153  * so that
 154  *      A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
 155  */
 156 static int transform_column(
 157         struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 158         unsigned row, unsigned col)
 159 {
 160         int i, t;
 161         isl_int a, b, g;
 162         unsigned total;
 163
 164         for (t = row-1; t >= 0; --t)
 165                 if (isl_int_ne(bset1->eq[t][col], bset2->eq[t][col]))
 166                         break;
 167         if (t < 0)
 168                 return 0;
 169
 170         total = 1 + bset1->dim;
 171         isl_int_init(a);
 172         isl_int_init(b);
 173         isl_int_init(g);
 174         isl_int_sub(b, bset1->eq[t][col], bset2->eq[t][col]);
 175         for (i = 0; i < t; ++i) {
 176                 isl_int_sub(a, bset2->eq[i][col], bset1->eq[i][col]);
 177                 isl_int_gcd(g, a, b);
 178                 isl_int_divexact(a, a, g);
 179                 isl_int_divexact(g, b, g);
 180                 isl_seq_combine(bset1->eq[i], g, bset1->eq[i], a, bset1->eq[t],
 181                                 total);
 182                 isl_seq_combine(bset2->eq[i], g, bset2->eq[i], a, bset2->eq[t],
 183                                 total);
 184         }
 185         isl_int_clear(a);
 186         isl_int_clear(b);
 187         isl_int_clear(g);
 188         delete_row(bset1, t);
 189         delete_row(bset2, t);
 190         return 1;
 191 }
 192
 193 /* The implementation is based on Section 5.2 of Michael Karr,
 194  * "Affine Relationships Among Variables of a Program",
 195  * except that the echelon form we use starts from the last column
 196  * and that we are dealing with integer coefficients.
 197  */
 198 static struct isl_basic_set *affine_hull(struct isl_ctx *ctx,
 199         struct isl_basic_set *bset1, struct isl_basic_set *bset2)
 200 {
 201         unsigned total;
 202         int col;
 203         int row;
 204
 205         total = 1 + bset1->dim;
 206
 207         row = 0;
 208         for (col = total-1; col >= 0; --col) {
 209                 int is_zero1 = row >= bset1->n_eq ||
 210                         isl_int_is_zero(bset1->eq[row][col]);
 211                 int is_zero2 = row >= bset2->n_eq ||
 212                         isl_int_is_zero(bset2->eq[row][col]);
 213                 if (!is_zero1 && !is_zero2) {
 214                         set_common_multiple(bset1, bset2, row, col);
 215                         ++row;
 216                 } else if (!is_zero1 && is_zero2) {
 217                         construct_column(bset1, bset2, row, col);
 218                 } else if (is_zero1 && !is_zero2) {
 219                         construct_column(bset2, bset1, row, col);
 220                 } else {
 221                         if (transform_column(bset1, bset2, row, col))
 222                                 --row;
 223                 }
 224         }
 225         isl_basic_set_free(ctx, bset2);
 226         return bset1;
 227 }
 228
 229 struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
 230                                                 struct isl_map *map)
 231 {
 232         int i;
 233         struct isl_basic_map *model = NULL;
 234         struct isl_basic_map *hull = NULL;
 235         struct isl_set *set;
 236
 237         if (!map)
 238                 return NULL;
 239
 240         if (map->n == 0) {
 241                 hull = isl_basic_map_empty(ctx,
 242                                             map->nparam, map->n_in, map->n_out);
 243                 isl_map_free(ctx, map);
 244                 return hull;
 245         }
 246
 247         map = isl_map_align_divs(ctx, map);
 248         model = isl_basic_map_copy(ctx, map->p[0]);
 249         set = isl_map_underlying_set(ctx, map);
 250         set = isl_set_cow(ctx, set);
 251         if (!set)
 252                 goto error;
 253
 254         for (i = 0; i < set->n; ++i) {
 255                 set->p[i] = isl_basic_set_cow(ctx, set->p[i]);
 256                 set->p[i] = isl_basic_set_affine_hull(ctx, set->p[i]);
 257                 set->p[i] = isl_basic_set_gauss(ctx, set->p[i], NULL);
 258                 if (!set->p[i])
 259                         goto error;
 260         }
 261         set = isl_set_remove_empty_parts(ctx, set);
 262         if (set->n == 0) {
 263                 hull = isl_basic_map_empty(ctx,
 264                                     model->nparam, model->n_in, model->n_out);
 265                 isl_basic_map_free(ctx, model);
 266         } else {
 267                 struct isl_basic_set *bset;
 268                 while (set->n > 1) {
 269                         set->p[0] = affine_hull(ctx, set->p[0],
 270                                                         set->p[--set->n]);
 271                         if (!set->p[0])
 272                                 goto error;
 273                 }
 274                 bset = isl_basic_set_copy(ctx, set->p[0]);
 275                 hull = isl_basic_map_overlying_set(ctx, bset, model);
 276         }
 277         isl_set_free(ctx, set);
 278         hull = isl_basic_map_simplify(ctx, hull);
 279         return isl_basic_map_finalize(ctx, hull);
 280 error:
 281         isl_basic_map_free(ctx, model);
 282         isl_set_free(ctx, set);
 283         return NULL;
 284 }
 285
 286 struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
 287                                                 struct isl_set *set)
 288 {
 289         return (struct isl_basic_set *)
 290                 isl_map_affine_hull(ctx, (struct isl_map *)set);
 291 }