2 * Copyright 2005-2007 Universiteit Leiden
3 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Copyright 2010 INRIA Saclay
6 * Use of this software is governed by the GNU LGPLv2.1 license
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
16 #include <isl_factorization.h>
17 #include <isl_dim_private.h>
18 #include <isl_mat_private.h>
20 static __isl_give isl_factorizer *isl_factorizer_alloc(
21 __isl_take isl_morph *morph, int n_group)
23 isl_factorizer *f = NULL;
30 len = isl_alloc_array(morph->dom->ctx, int, n_group);
35 f = isl_alloc_type(morph->dom->ctx, struct isl_factorizer);
46 isl_morph_free(morph);
50 void isl_factorizer_free(__isl_take isl_factorizer *f)
55 isl_morph_free(f->morph);
60 void isl_factorizer_dump(__isl_take isl_factorizer *f, FILE *out)
67 isl_morph_dump(f->morph, out);
69 for (i = 0; i < f->n_group; ++i) {
72 fprintf(out, "%d", f->len[i]);
77 __isl_give isl_factorizer *isl_factorizer_identity(__isl_keep isl_basic_set *bset)
79 return isl_factorizer_alloc(isl_morph_identity(bset), 0);
82 __isl_give isl_factorizer *isl_factorizer_groups(__isl_keep isl_basic_set *bset,
83 __isl_take isl_mat *Q, __isl_take isl_mat *U, int n, int *len)
95 if (!bset || !Q || !U)
98 ovar = 1 + isl_dim_offset(bset->dim, isl_dim_set);
99 id = isl_mat_identity(bset->ctx, ovar);
100 Q = isl_mat_diagonal(isl_mat_copy(id), Q);
101 U = isl_mat_diagonal(id, U);
103 nvar = isl_basic_set_dim(bset, isl_dim_set);
104 dim = isl_basic_set_get_dim(bset);
105 dom = isl_basic_set_universe(isl_dim_copy(dim));
106 dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
107 dim = isl_dim_add(dim, isl_dim_set, nvar);
108 ran = isl_basic_set_universe(dim);
109 morph = isl_morph_alloc(dom, ran, Q, U);
110 f = isl_factorizer_alloc(morph, n);
113 for (i = 0; i < n; ++i)
122 struct isl_factor_groups {
123 int *pos; /* for each column: row position of pivot */
124 int *group; /* group to which a column belongs */
125 int *cnt; /* number of columns in the group */
126 int *rowgroup; /* group to which a constraint belongs */
129 /* Initialize isl_factor_groups structure: find pivot row positions,
130 * each column initially belongs to its own group and the groups
131 * of the constraints are still unknown.
133 static int init_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
140 g->pos = isl_alloc_array(H->ctx, int, H->n_col);
141 g->group = isl_alloc_array(H->ctx, int, H->n_col);
142 g->cnt = isl_alloc_array(H->ctx, int, H->n_col);
143 g->rowgroup = isl_alloc_array(H->ctx, int, H->n_row);
145 if (!g->pos || !g->group || !g->cnt || !g->rowgroup)
148 for (i = 0; i < H->n_row; ++i)
150 for (i = 0, j = 0; i < H->n_col; ++i) {
151 for ( ; j < H->n_row; ++j)
152 if (!isl_int_is_zero(H->row[j][i]))
156 for (i = 0; i < H->n_col; ++i) {
164 /* Update group[k] to the group column k belongs to.
165 * When merging two groups, only the group of the current
166 * group leader is changed. Here we change the group of
167 * the other members to also point to the group that the
168 * old group leader now points to.
170 static void update_group(struct isl_factor_groups *g, int k)
173 while (g->cnt[p] == 0)
178 /* Merge group i with all groups of the subsequent columns
179 * with non-zero coefficients in row j of H.
180 * (The previous columns are all zero; otherwise we would have handled
183 static int update_group_i_with_row_j(struct isl_factor_groups *g, int i, int j,
184 __isl_keep isl_mat *H)
188 g->rowgroup[j] = g->group[i];
189 for (k = i + 1; k < H->n_col && j >= g->pos[k]; ++k) {
192 if (g->group[k] != g->group[i] &&
193 !isl_int_is_zero(H->row[j][k])) {
194 isl_assert(H->ctx, g->cnt[g->group[k]] != 0, return -1);
195 isl_assert(H->ctx, g->cnt[g->group[i]] != 0, return -1);
196 if (g->group[i] < g->group[k]) {
197 g->cnt[g->group[i]] += g->cnt[g->group[k]];
198 g->cnt[g->group[k]] = 0;
199 g->group[g->group[k]] = g->group[i];
201 g->cnt[g->group[k]] += g->cnt[g->group[i]];
202 g->cnt[g->group[i]] = 0;
203 g->group[g->group[i]] = g->group[k];
211 /* Update the group information based on the constraint matrix.
213 static int update_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
217 for (i = 0; i < H->n_col && g->cnt[0] < H->n_col; ++i) {
218 if (g->pos[i] == H->n_row)
219 continue; /* A line direction */
220 if (g->rowgroup[g->pos[i]] == -1)
221 g->rowgroup[g->pos[i]] = i;
222 for (j = g->pos[i] + 1; j < H->n_row; ++j) {
223 if (isl_int_is_zero(H->row[j][i]))
225 if (g->rowgroup[j] != -1)
227 if (update_group_i_with_row_j(g, i, j, H) < 0)
235 static void clear_groups(struct isl_factor_groups *g)
245 /* Determine if the set variables of the basic set can be factorized and
246 * return the results in an isl_factorizer.
248 * The algorithm works by first computing the Hermite normal form
249 * and then grouping columns linked by one or more constraints together,
250 * where a constraints "links" two or more columns if the constraint
251 * has nonzero coefficients in the columns.
253 __isl_give isl_factorizer *isl_basic_set_factorizer(
254 __isl_keep isl_basic_set *bset)
259 struct isl_factor_groups g = { 0 };
265 isl_assert(bset->ctx, isl_basic_set_dim(bset, isl_dim_div) == 0,
268 nvar = isl_basic_set_dim(bset, isl_dim_set);
270 return isl_factorizer_identity(bset);
272 H = isl_mat_alloc(bset->ctx, bset->n_eq + bset->n_ineq, nvar);
275 isl_mat_sub_copy(bset->ctx, H->row, bset->eq, bset->n_eq,
276 0, 1 + isl_dim_offset(bset->dim, isl_dim_set), nvar);
277 isl_mat_sub_copy(bset->ctx, H->row + bset->n_eq, bset->ineq, bset->n_ineq,
278 0, 1 + isl_dim_offset(bset->dim, isl_dim_set), nvar);
279 H = isl_mat_left_hermite(H, 0, &U, &Q);
281 if (init_groups(&g, H) < 0)
283 if (update_groups(&g, H) < 0)
286 if (g.cnt[0] == nvar) {
292 return isl_factorizer_identity(bset);
297 while (done != nvar) {
298 int group = g.group[done];
299 for (i = 1; i < g.cnt[group]; ++i) {
300 if (g.group[done + i] == group)
302 for (j = done + g.cnt[group]; j < nvar; ++j)
303 if (g.group[j] == group)
305 g.group[j] = g.group[done + i];
306 Q = isl_mat_swap_rows(Q, done + i, j);
307 U = isl_mat_swap_cols(U, done + i, j);
309 done += g.cnt[group];
310 g.pos[n++] = g.cnt[group];
313 f = isl_factorizer_groups(bset, Q, U, n, g.pos);