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_affine_hull(struct isl_ctx *ctx,
   9                                                 struct isl_basic_map *bmap)
  10 {
  11         int i;
  12         isl_int opt;
  13
  14         bmap = isl_basic_map_cow(ctx, bmap);
  15         if (!bmap)
  16                 return NULL;
  17
  18         isl_int_init(opt);
  19         for (i = 0; i < bmap->n_ineq; ++i) {
  20                 enum isl_lp_result res;
  21                 res = isl_solve_lp(bmap, 1, bmap->ineq[i]+1, ctx->one, &opt);
  22                 if (res == isl_lp_unbounded)
  23                         continue;
  24                 if (res == isl_lp_error)
  25                         goto error;
  26                 if (res == isl_lp_empty) {
  27                         bmap = isl_basic_map_set_to_empty(ctx, bmap);
  28                         break;
  29                 }
  30                 isl_int_add(opt, opt, bmap->ineq[i][0]);
  31                 if (isl_int_is_zero(opt)) {
  32                         isl_basic_map_inequality_to_equality(ctx, bmap, i);
  33                         --i;
  34                 }
  35         }
  36         isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
  37         isl_int_clear(opt);
  38         return isl_basic_map_finalize(ctx, bmap);
  39 error:
  40         isl_int_clear(opt);
  41         isl_basic_map_free(ctx, bmap);
  42         return NULL;
  43 }
  44
  45 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
  46                                                 struct isl_basic_set *bset)
  47 {
  48         return (struct isl_basic_set *)
  49                 isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
  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_map *bmap1, struct isl_basic_map *bmap2,
  59         unsigned row, unsigned col)
  60 {
  61         isl_int m, c;
  62
  63         if (isl_int_eq(bmap1->eq[row][col], bmap2->eq[row][col]))
  64                 return;
  65
  66         isl_int_init(c);
  67         isl_int_init(m);
  68         isl_int_lcm(m, bmap1->eq[row][col], bmap2->eq[row][col]);
  69         isl_int_divexact(c, m, bmap1->eq[row][col]);
  70         isl_seq_scale(bmap1->eq[row], bmap1->eq[row], c, col+1);
  71         isl_int_divexact(c, m, bmap2->eq[row][col]);
  72         isl_seq_scale(bmap2->eq[row], bmap2->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_map *bmap, unsigned row)
  80 {
  81         isl_int *t;
  82         int r;
  83
  84         t = bmap->eq[row];
  85         bmap->n_eq--;
  86         for (r = row; r < bmap->n_eq; ++r)
  87                 bmap->eq[r] = bmap->eq[r+1];
  88         bmap->eq[bmap->n_eq] = t;
  89 }
  90
  91 /* Make first row entries in column col of bmap1 identical to
  92  * those of bmap2, using the fact that entry bmap1->eq[row][col]=a
  93  * is non-zero.  Initially, these elements of bmap1 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_map *bmap1, struct isl_basic_map *bmap2,
 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 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
 112         for (r = 0; r < row; ++r) {
 113                 if (isl_int_is_zero(bmap2->eq[r][col]))
 114                         continue;
 115                 isl_int_gcd(b, bmap2->eq[r][col], bmap1->eq[row][col]);
 116                 isl_int_divexact(a, bmap1->eq[row][col], b);
 117                 isl_int_divexact(b, bmap2->eq[r][col], b);
 118                 isl_seq_combine(bmap1->eq[r], a, bmap1->eq[r],
 119                                               b, bmap1->eq[row], total);
 120                 isl_seq_scale(bmap2->eq[r], bmap2->eq[r], a, total);
 121         }
 122         isl_int_clear(a);
 123         isl_int_clear(b);
 124         delete_row(bmap1, row);
 125 }
 126
 127 /* Make first row entries in column col of bmap1 identical to
 128  * those of bmap2, 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_map *bmap1, struct isl_basic_map *bmap2,
 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(bmap1->eq[t][col], bmap2->eq[t][col]))
 146                         break;
 147         if (t < 0)
 148                 return 0;
 149
 150         total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
 151         isl_int_init(a);
 152         isl_int_init(b);
 153         isl_int_init(g);
 154         isl_int_sub(b, bmap1->eq[t][col], bmap2->eq[t][col]);
 155         for (i = 0; i < t; ++i) {
 156                 isl_int_sub(a, bmap2->eq[i][col], bmap1->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(bmap1->eq[i], g, bmap1->eq[i], a, bmap1->eq[t],
 161                                 total);
 162                 isl_seq_combine(bmap2->eq[i], g, bmap2->eq[i], a, bmap2->eq[t],
 163                                 total);
 164         }
 165         isl_int_clear(a);
 166         isl_int_clear(b);
 167         isl_int_clear(g);
 168         delete_row(bmap1, t);
 169         delete_row(bmap2, 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_map *affine_hull(struct isl_ctx *ctx,
 179         struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
 180 {
 181         unsigned total;
 182         int col;
 183         int row;
 184
 185         total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
 186
 187         row = 0;
 188         for (col = total-1; col >= 0; --col) {
 189                 int is_zero1 = row >= bmap1->n_eq ||
 190                         isl_int_is_zero(bmap1->eq[row][col]);
 191                 int is_zero2 = row >= bmap2->n_eq ||
 192                         isl_int_is_zero(bmap2->eq[row][col]);
 193                 if (!is_zero1 && !is_zero2) {
 194                         set_common_multiple(bmap1, bmap2, row, col);
 195                         ++row;
 196                 } else if (!is_zero1 && is_zero2) {
 197                         construct_column(bmap1, bmap2, row, col);
 198                 } else if (is_zero1 && !is_zero2) {
 199                         construct_column(bmap2, bmap1, row, col);
 200                 } else {
 201                         if (transform_column(bmap1, bmap2, row, col))
 202                                 --row;
 203                 }
 204         }
 205         isl_basic_map_free(ctx, bmap2);
 206         return bmap1;
 207 }
 208
 209 struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
 210                                                 struct isl_map *map)
 211 {
 212         int i;
 213         struct isl_basic_map *bmap;
 214
 215         map = isl_map_compute_divs(ctx, map);
 216         map = isl_map_cow(ctx, map);
 217         if (!map)
 218                 return NULL;
 219
 220         if (map->n == 0) {
 221                 bmap = isl_basic_map_empty(ctx,
 222                                             map->nparam, map->n_in, map->n_out);
 223                 isl_map_free(ctx, map);
 224                 return bmap;
 225         }
 226
 227         for (i = 1; i < map->n; ++i)
 228                 map->p[0] = isl_basic_map_align_divs(ctx, map->p[0], map->p[i]);
 229         for (i = 1; i < map->n; ++i)
 230                 map->p[i] = isl_basic_map_align_divs(ctx, map->p[i], map->p[0]);
 231
 232         for (i = 0; i < map->n; ++i) {
 233                 map->p[i] = isl_basic_map_cow(ctx, map->p[i]);
 234                 map->p[i] = isl_basic_map_affine_hull(ctx, map->p[i]);
 235                 map->p[i] = isl_basic_map_gauss(ctx, map->p[i], NULL);
 236                 if (!map->p[i])
 237                         goto error;
 238         }
 239         while (map->n > 1) {
 240                 map->p[0] = affine_hull(ctx, map->p[0], map->p[--map->n]);
 241                 if (!map->p[0])
 242                         goto error;
 243         }
 244         bmap = isl_basic_map_copy(ctx, map->p[0]);
 245         isl_map_free(ctx, map);
 246         bmap = isl_basic_map_finalize(ctx, bmap);
 247         return isl_basic_map_simplify(ctx, bmap);
 248 error:
 249         isl_map_free(ctx, map);
 250         return NULL;
 251 }
 252
 253 struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
 254                                                 struct isl_set *set)
 255 {
 256         return (struct isl_basic_set *)
 257                 isl_map_affine_hull(ctx, (struct isl_map *)set);
 258 }