1 /*
2 * Copyright 2005-2007 Universiteit Leiden
3 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Copyright 2010 INRIA Saclay
5 *
6 * Use of this software is governed by the MIT license
7 *
8 * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science,
9 * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands
10 * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A,
11 * B-3001 Leuven, Belgium
12 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
13 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
14 */
15
16 #include <isl_map_private.h>
17 #include <isl_factorization.h>
18 #include <isl_space_private.h>
19 #include <isl_mat_private.h>
20
21 /* Return the isl_ctx to which "f" belongs.
22 */
isl_factorizer_get_ctx(__isl_keep isl_factorizer * f)23 isl_ctx *isl_factorizer_get_ctx(__isl_keep isl_factorizer *f)
24 {
25 if (!f)
26 return NULL;
27 return isl_basic_set_get_ctx(f->bset);
28 }
29
isl_factorizer_alloc(__isl_keep isl_basic_set * bset,__isl_take isl_morph * morph,int n_group)30 static __isl_give isl_factorizer *isl_factorizer_alloc(
31 __isl_keep isl_basic_set *bset, __isl_take isl_morph *morph,
32 int n_group)
33 {
34 isl_factorizer *f = NULL;
35 int *len = NULL;
36
37 if (!morph)
38 return NULL;
39
40 if (n_group > 0) {
41 len = isl_alloc_array(morph->dom->ctx, int, n_group);
42 if (!len)
43 goto error;
44 }
45
46 f = isl_alloc_type(morph->dom->ctx, struct isl_factorizer);
47 if (!f)
48 goto error;
49
50 f->bset = isl_basic_set_copy(bset);
51 f->morph = morph;
52 f->n_group = n_group;
53 f->len = len;
54
55 return f;
56 error:
57 free(len);
58 isl_morph_free(morph);
59 return NULL;
60 }
61
isl_factorizer_free(__isl_take isl_factorizer * f)62 __isl_null isl_factorizer *isl_factorizer_free(__isl_take isl_factorizer *f)
63 {
64 if (!f)
65 return NULL;
66
67 isl_basic_set_free(f->bset);
68 isl_morph_free(f->morph);
69 free(f->len);
70 free(f);
71 return NULL;
72 }
73
isl_factorizer_dump(__isl_take isl_factorizer * f)74 void isl_factorizer_dump(__isl_take isl_factorizer *f)
75 {
76 int i;
77
78 if (!f)
79 return;
80
81 isl_morph_print_internal(f->morph, stderr);
82 fprintf(stderr, "[");
83 for (i = 0; i < f->n_group; ++i) {
84 if (i)
85 fprintf(stderr, ", ");
86 fprintf(stderr, "%d", f->len[i]);
87 }
88 fprintf(stderr, "]\n");
89 }
90
isl_factorizer_identity(__isl_keep isl_basic_set * bset)91 __isl_give isl_factorizer *isl_factorizer_identity(__isl_keep isl_basic_set *bset)
92 {
93 return isl_factorizer_alloc(bset, isl_morph_identity(bset), 0);
94 }
95
isl_factorizer_groups(__isl_keep isl_basic_set * bset,__isl_take isl_mat * Q,__isl_take isl_mat * U,int n,int * len)96 __isl_give isl_factorizer *isl_factorizer_groups(__isl_keep isl_basic_set *bset,
97 __isl_take isl_mat *Q, __isl_take isl_mat *U, int n, int *len)
98 {
99 int i;
100 isl_size nvar;
101 unsigned ovar;
102 isl_space *space;
103 isl_basic_set *dom;
104 isl_basic_set *ran;
105 isl_morph *morph;
106 isl_factorizer *f;
107 isl_mat *id;
108
109 nvar = isl_basic_set_dim(bset, isl_dim_set);
110 if (nvar < 0 || !Q || !U)
111 goto error;
112
113 ovar = 1 + isl_space_offset(bset->dim, isl_dim_set);
114 id = isl_mat_identity(bset->ctx, ovar);
115 Q = isl_mat_diagonal(isl_mat_copy(id), Q);
116 U = isl_mat_diagonal(id, U);
117
118 space = isl_basic_set_get_space(bset);
119 dom = isl_basic_set_universe(isl_space_copy(space));
120 space = isl_space_drop_dims(space, isl_dim_set, 0, nvar);
121 space = isl_space_add_dims(space, isl_dim_set, nvar);
122 ran = isl_basic_set_universe(space);
123 morph = isl_morph_alloc(dom, ran, Q, U);
124 f = isl_factorizer_alloc(bset, morph, n);
125 if (!f)
126 return NULL;
127 for (i = 0; i < n; ++i)
128 f->len[i] = len[i];
129 return f;
130 error:
131 isl_mat_free(Q);
132 isl_mat_free(U);
133 return NULL;
134 }
135
136 struct isl_factor_groups {
137 int *pos; /* for each column: row position of pivot */
138 int *group; /* group to which a column belongs */
139 int *cnt; /* number of columns in the group */
140 int *rowgroup; /* group to which a constraint belongs */
141 };
142
143 /* Initialize isl_factor_groups structure: find pivot row positions,
144 * each column initially belongs to its own group and the groups
145 * of the constraints are still unknown.
146 */
init_groups(struct isl_factor_groups * g,__isl_keep isl_mat * H)147 static int init_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
148 {
149 int i, j;
150
151 if (!H)
152 return -1;
153
154 g->pos = isl_alloc_array(H->ctx, int, H->n_col);
155 g->group = isl_alloc_array(H->ctx, int, H->n_col);
156 g->cnt = isl_alloc_array(H->ctx, int, H->n_col);
157 g->rowgroup = isl_alloc_array(H->ctx, int, H->n_row);
158
159 if (!g->pos || !g->group || !g->cnt || !g->rowgroup)
160 return -1;
161
162 for (i = 0; i < H->n_row; ++i)
163 g->rowgroup[i] = -1;
164 for (i = 0, j = 0; i < H->n_col; ++i) {
165 for ( ; j < H->n_row; ++j)
166 if (!isl_int_is_zero(H->row[j][i]))
167 break;
168 g->pos[i] = j;
169 }
170 for (i = 0; i < H->n_col; ++i) {
171 g->group[i] = i;
172 g->cnt[i] = 1;
173 }
174
175 return 0;
176 }
177
178 /* Update group[k] to the group column k belongs to.
179 * When merging two groups, only the group of the current
180 * group leader is changed. Here we change the group of
181 * the other members to also point to the group that the
182 * old group leader now points to.
183 */
update_group(struct isl_factor_groups * g,int k)184 static void update_group(struct isl_factor_groups *g, int k)
185 {
186 int p = g->group[k];
187 while (g->cnt[p] == 0)
188 p = g->group[p];
189 g->group[k] = p;
190 }
191
192 /* Merge group i with all groups of the subsequent columns
193 * with non-zero coefficients in row j of H.
194 * (The previous columns are all zero; otherwise we would have handled
195 * the row before.)
196 */
update_group_i_with_row_j(struct isl_factor_groups * g,int i,int j,__isl_keep isl_mat * H)197 static int update_group_i_with_row_j(struct isl_factor_groups *g, int i, int j,
198 __isl_keep isl_mat *H)
199 {
200 int k;
201
202 g->rowgroup[j] = g->group[i];
203 for (k = i + 1; k < H->n_col && j >= g->pos[k]; ++k) {
204 update_group(g, k);
205 update_group(g, i);
206 if (g->group[k] != g->group[i] &&
207 !isl_int_is_zero(H->row[j][k])) {
208 isl_assert(H->ctx, g->cnt[g->group[k]] != 0, return -1);
209 isl_assert(H->ctx, g->cnt[g->group[i]] != 0, return -1);
210 if (g->group[i] < g->group[k]) {
211 g->cnt[g->group[i]] += g->cnt[g->group[k]];
212 g->cnt[g->group[k]] = 0;
213 g->group[g->group[k]] = g->group[i];
214 } else {
215 g->cnt[g->group[k]] += g->cnt[g->group[i]];
216 g->cnt[g->group[i]] = 0;
217 g->group[g->group[i]] = g->group[k];
218 }
219 }
220 }
221
222 return 0;
223 }
224
225 /* Update the group information based on the constraint matrix.
226 */
update_groups(struct isl_factor_groups * g,__isl_keep isl_mat * H)227 static int update_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
228 {
229 int i, j;
230
231 for (i = 0; i < H->n_col && g->cnt[0] < H->n_col; ++i) {
232 if (g->pos[i] == H->n_row)
233 continue; /* A line direction */
234 if (g->rowgroup[g->pos[i]] == -1)
235 g->rowgroup[g->pos[i]] = i;
236 for (j = g->pos[i] + 1; j < H->n_row; ++j) {
237 if (isl_int_is_zero(H->row[j][i]))
238 continue;
239 if (g->rowgroup[j] != -1)
240 continue;
241 if (update_group_i_with_row_j(g, i, j, H) < 0)
242 return -1;
243 }
244 }
245 for (i = 1; i < H->n_col; ++i)
246 update_group(g, i);
247
248 return 0;
249 }
250
clear_groups(struct isl_factor_groups * g)251 static void clear_groups(struct isl_factor_groups *g)
252 {
253 if (!g)
254 return;
255 free(g->pos);
256 free(g->group);
257 free(g->cnt);
258 free(g->rowgroup);
259 }
260
261 /* Determine if the set variables of the basic set can be factorized and
262 * return the results in an isl_factorizer.
263 *
264 * The algorithm works by first computing the Hermite normal form
265 * and then grouping columns linked by one or more constraints together,
266 * where a constraints "links" two or more columns if the constraint
267 * has nonzero coefficients in the columns.
268 */
isl_basic_set_factorizer(__isl_keep isl_basic_set * bset)269 __isl_give isl_factorizer *isl_basic_set_factorizer(
270 __isl_keep isl_basic_set *bset)
271 {
272 int i, j, n, done;
273 isl_mat *H, *U, *Q;
274 isl_size nvar;
275 struct isl_factor_groups g = { 0 };
276 isl_factorizer *f;
277
278 nvar = isl_basic_set_dim(bset, isl_dim_set);
279 if (nvar < 0 || isl_basic_set_check_no_locals(bset) < 0)
280 return NULL;
281
282 if (nvar <= 1)
283 return isl_factorizer_identity(bset);
284
285 H = isl_mat_alloc(bset->ctx, bset->n_eq + bset->n_ineq, nvar);
286 if (!H)
287 return NULL;
288 isl_mat_sub_copy(bset->ctx, H->row, bset->eq, bset->n_eq,
289 0, 1 + isl_space_offset(bset->dim, isl_dim_set), nvar);
290 isl_mat_sub_copy(bset->ctx, H->row + bset->n_eq, bset->ineq, bset->n_ineq,
291 0, 1 + isl_space_offset(bset->dim, isl_dim_set), nvar);
292 H = isl_mat_left_hermite(H, 0, &U, &Q);
293
294 if (init_groups(&g, H) < 0)
295 goto error;
296 if (update_groups(&g, H) < 0)
297 goto error;
298
299 if (g.cnt[0] == nvar) {
300 isl_mat_free(H);
301 isl_mat_free(U);
302 isl_mat_free(Q);
303 clear_groups(&g);
304
305 return isl_factorizer_identity(bset);
306 }
307
308 done = 0;
309 n = 0;
310 while (done != nvar) {
311 int group = g.group[done];
312 for (i = 1; i < g.cnt[group]; ++i) {
313 if (g.group[done + i] == group)
314 continue;
315 for (j = done + g.cnt[group]; j < nvar; ++j)
316 if (g.group[j] == group)
317 break;
318 if (j == nvar)
319 isl_die(bset->ctx, isl_error_internal,
320 "internal error", goto error);
321 g.group[j] = g.group[done + i];
322 Q = isl_mat_swap_rows(Q, done + i, j);
323 U = isl_mat_swap_cols(U, done + i, j);
324 }
325 done += g.cnt[group];
326 g.pos[n++] = g.cnt[group];
327 }
328
329 f = isl_factorizer_groups(bset, Q, U, n, g.pos);
330
331 isl_mat_free(H);
332 clear_groups(&g);
333
334 return f;
335 error:
336 isl_mat_free(H);
337 isl_mat_free(U);
338 isl_mat_free(Q);
339 clear_groups(&g);
340 return NULL;
341 }
342
343 /* Given the factorizer "f" of a basic set,
344 * call "test" on each resulting factor as long as each call succeeds.
345 */
isl_factorizer_every_factor_basic_set(__isl_keep isl_factorizer * f,isl_bool (* test)(__isl_keep isl_basic_set * bset,void * user),void * user)346 __isl_give isl_bool isl_factorizer_every_factor_basic_set(
347 __isl_keep isl_factorizer *f,
348 isl_bool (*test)(__isl_keep isl_basic_set *bset, void *user),
349 void *user)
350 {
351 int i, n;
352 isl_bool every = isl_bool_true;
353 isl_size nparam, nvar;
354 isl_basic_set *bset;
355
356 if (!f)
357 return isl_bool_error;
358 nparam = isl_basic_set_dim(f->bset, isl_dim_param);
359 nvar = isl_basic_set_dim(f->bset, isl_dim_set);
360 if (nparam < 0 || nvar < 0)
361 return isl_bool_error;
362
363 bset = isl_basic_set_copy(f->bset);
364 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
365
366 for (i = 0, n = 0; i < f->n_group; ++i) {
367 isl_basic_set *factor;
368
369 factor = isl_basic_set_copy(bset);
370 factor = isl_basic_set_drop_constraints_involving(factor,
371 nparam + n + f->len[i], nvar - n - f->len[i]);
372 factor = isl_basic_set_drop_constraints_involving(factor,
373 nparam, n);
374 factor = isl_basic_set_drop(factor, isl_dim_set,
375 n + f->len[i], nvar - n - f->len[i]);
376 factor = isl_basic_set_drop(factor, isl_dim_set, 0, n);
377 every = test(factor, user);
378 isl_basic_set_free(factor);
379
380 if (every < 0 || !every)
381 break;
382
383 n += f->len[i];
384 }
385
386 isl_basic_set_free(bset);
387
388 return every;
389 }
390