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_map *bmap1, struct isl_basic_map *bmap2,
  79         unsigned row, unsigned col)
  80 {
  81         isl_int m, c;
  82
  83         if (isl_int_eq(bmap1->eq[row][col], bmap2->eq[row][col]))
  84                 return;
  85
  86         isl_int_init(c);
  87         isl_int_init(m);
  88         isl_int_lcm(m, bmap1->eq[row][col], bmap2->eq[row][col]);
  89         isl_int_divexact(c, m, bmap1->eq[row][col]);
  90         isl_seq_scale(bmap1->eq[row], bmap1->eq[row], c, col+1);
  91         isl_int_divexact(c, m, bmap2->eq[row][col]);
  92         isl_seq_scale(bmap2->eq[row], bmap2->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_map *bmap, unsigned row)
 100 {
 101         isl_int *t;
 102         int r;
 103
 104         t = bmap->eq[row];
 105         bmap->n_eq--;
 106         for (r = row; r < bmap->n_eq; ++r)
 107                 bmap->eq[r] = bmap->eq[r+1];
 108         bmap->eq[bmap->n_eq] = t;
 109 }
 110
 111 /* Make first row entries in column col of bmap1 identical to
 112  * those of bmap2, using the fact that entry bmap1->eq[row][col]=a
 113  * is non-zero.  Initially, these elements of bmap1 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_map *bmap1, struct isl_basic_map *bmap2,
 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 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
 132         for (r = 0; r < row; ++r) {
 133                 if (isl_int_is_zero(bmap2->eq[r][col]))
 134                         continue;
 135                 isl_int_gcd(b, bmap2->eq[r][col], bmap1->eq[row][col]);
 136                 isl_int_divexact(a, bmap1->eq[row][col], b);
 137                 isl_int_divexact(b, bmap2->eq[r][col], b);
 138                 isl_seq_combine(bmap1->eq[r], a, bmap1->eq[r],
 139                                               b, bmap1->eq[row], total);
 140                 isl_seq_scale(bmap2->eq[r], bmap2->eq[r], a, total);
 141         }
 142         isl_int_clear(a);
 143         isl_int_clear(b);
 144         delete_row(bmap1, row);
 145 }
 146
 147 /* Make first row entries in column col of bmap1 identical to
 148  * those of bmap2, 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_map *bmap1, struct isl_basic_map *bmap2,
 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(bmap1->eq[t][col], bmap2->eq[t][col]))
 166                         break;
 167         if (t < 0)
 168                 return 0;
 169
 170         total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
 171         isl_int_init(a);
 172         isl_int_init(b);
 173         isl_int_init(g);
 174         isl_int_sub(b, bmap1->eq[t][col], bmap2->eq[t][col]);
 175         for (i = 0; i < t; ++i) {
 176                 isl_int_sub(a, bmap2->eq[i][col], bmap1->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(bmap1->eq[i], g, bmap1->eq[i], a, bmap1->eq[t],
 181                                 total);
 182                 isl_seq_combine(bmap2->eq[i], g, bmap2->eq[i], a, bmap2->eq[t],
 183                                 total);
 184         }
 185         isl_int_clear(a);
 186         isl_int_clear(b);
 187         isl_int_clear(g);
 188         delete_row(bmap1, t);
 189         delete_row(bmap2, 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_map *affine_hull(struct isl_ctx *ctx,
 199         struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
 200 {
 201         unsigned total;
 202         int col;
 203         int row;
 204
 205         total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
 206
 207         row = 0;
 208         for (col = total-1; col >= 0; --col) {
 209                 int is_zero1 = row >= bmap1->n_eq ||
 210                         isl_int_is_zero(bmap1->eq[row][col]);
 211                 int is_zero2 = row >= bmap2->n_eq ||
 212                         isl_int_is_zero(bmap2->eq[row][col]);
 213                 if (!is_zero1 && !is_zero2) {
 214                         set_common_multiple(bmap1, bmap2, row, col);
 215                         ++row;
 216                 } else if (!is_zero1 && is_zero2) {
 217                         construct_column(bmap1, bmap2, row, col);
 218                 } else if (is_zero1 && !is_zero2) {
 219                         construct_column(bmap2, bmap1, row, col);
 220                 } else {
 221                         if (transform_column(bmap1, bmap2, row, col))
 222                                 --row;
 223                 }
 224         }
 225         isl_basic_map_free(ctx, bmap2);
 226         return bmap1;
 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 *bmap;
 234
 235         if (!map)
 236                 return NULL;
 237
 238         if (map->n == 0) {
 239                 bmap = isl_basic_map_empty(ctx,
 240                                             map->nparam, map->n_in, map->n_out);
 241                 isl_map_free(ctx, map);
 242                 return bmap;
 243         }
 244
 245         map = isl_map_align_divs(ctx, map);
 246         map = isl_map_cow(ctx, map);
 247
 248         for (i = 0; i < map->n; ++i) {
 249                 map->p[i] = isl_basic_map_cow(ctx, map->p[i]);
 250                 map->p[i] = isl_basic_map_affine_hull(ctx, map->p[i]);
 251                 map->p[i] = isl_basic_map_gauss(ctx, map->p[i], NULL);
 252                 if (!map->p[i])
 253                         goto error;
 254         }
 255         while (map->n > 1) {
 256                 map->p[0] = affine_hull(ctx, map->p[0], map->p[--map->n]);
 257                 if (!map->p[0])
 258                         goto error;
 259         }
 260         bmap = isl_basic_map_copy(ctx, map->p[0]);
 261         isl_map_free(ctx, map);
 262         bmap = isl_basic_map_finalize(ctx, bmap);
 263         return isl_basic_map_simplify(ctx, bmap);
 264 error:
 265         isl_map_free(ctx, map);
 266         return NULL;
 267 }
 268
 269 struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
 270                                                 struct isl_set *set)
 271 {
 272         return (struct isl_basic_set *)
 273                 isl_map_affine_hull(ctx, (struct isl_map *)set);
 274 }