6 #include "isl_map_private.h"
8 struct isl_basic_map *isl_basic_map_affine_hull(struct isl_ctx *ctx,
9 struct isl_basic_map *bmap)
14 bmap = isl_basic_map_cow(ctx, bmap);
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)
24 if (res == isl_lp_error)
26 if (res == isl_lp_empty) {
27 bmap = isl_basic_map_set_to_empty(ctx, bmap);
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);
36 isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
38 return isl_basic_map_finalize(ctx, bmap);
41 isl_basic_map_free(ctx, bmap);
45 struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
46 struct isl_basic_set *bset)
48 return (struct isl_basic_set *)
49 isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
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.
57 static void set_common_multiple(
58 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
59 unsigned row, unsigned col)
63 if (isl_int_eq(bmap1->eq[row][col], bmap2->eq[row][col]))
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);
77 /* Delete a given equality, moving all the following equalities one up.
79 static void delete_row(struct isl_basic_map *bmap, unsigned row)
86 for (r = row; r < bmap->n_eq; ++r)
87 bmap->eq[r] = bmap->eq[r+1];
88 bmap->eq[bmap->n_eq] = t;
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]
98 * A[i][col] = B[i][col] = a * old(B[i][col])
100 static void construct_column(
101 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
102 unsigned row, unsigned col)
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]))
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);
124 delete_row(bmap1, row);
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]
134 * A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
136 static int transform_column(
137 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
138 unsigned row, unsigned col)
144 for (t = row-1; t >= 0; --t)
145 if (isl_int_ne(bmap1->eq[t][col], bmap2->eq[t][col]))
150 total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
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],
162 isl_seq_combine(bmap2->eq[i], g, bmap2->eq[i], a, bmap2->eq[t],
168 delete_row(bmap1, t);
169 delete_row(bmap2, t);
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.
178 static struct isl_basic_map *affine_hull(struct isl_ctx *ctx,
179 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
185 total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
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);
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);
201 if (transform_column(bmap1, bmap2, row, col))
205 isl_basic_map_free(ctx, bmap2);
209 struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
213 struct isl_basic_map *bmap;
215 map = isl_map_compute_divs(ctx, map);
216 map = isl_map_cow(ctx, map);
221 bmap = isl_basic_map_empty(ctx,
222 map->nparam, map->n_in, map->n_out);
223 isl_map_free(ctx, map);
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]);
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);
240 map->p[0] = affine_hull(ctx, map->p[0], map->p[--map->n]);
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);
249 isl_map_free(ctx, map);
253 struct isl_basic_set *isl_set_affine_hull(struct isl_ctx *ctx,
256 return (struct isl_basic_set *)
257 isl_map_affine_hull(ctx, (struct isl_map *)set);