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
21static __isl_give isl_factorizer *isl_factorizer_alloc(
22	__isl_take isl_morph *morph, int n_group)
23{
24	isl_factorizer *f = NULL;
25	int *len = NULL;
26
27	if (!morph)
28		return NULL;
29
30	if (n_group > 0) {
31		len = isl_alloc_array(morph->dom->ctx, int, n_group);
32		if (!len)
33			goto error;
34	}
35
36	f = isl_alloc_type(morph->dom->ctx, struct isl_factorizer);
37	if (!f)
38		goto error;
39
40	f->morph = morph;
41	f->n_group = n_group;
42	f->len = len;
43
44	return f;
45error:
46	free(len);
47	isl_morph_free(morph);
48	return NULL;
49}
50
51void isl_factorizer_free(__isl_take isl_factorizer *f)
52{
53	if (!f)
54		return;
55
56	isl_morph_free(f->morph);
57	free(f->len);
58	free(f);
59}
60
61void isl_factorizer_dump(__isl_take isl_factorizer *f)
62{
63	int i;
64
65	if (!f)
66		return;
67
68	isl_morph_print_internal(f->morph, stderr);
69	fprintf(stderr, "[");
70	for (i = 0; i < f->n_group; ++i) {
71		if (i)
72			fprintf(stderr, ", ");
73		fprintf(stderr, "%d", f->len[i]);
74	}
75	fprintf(stderr, "]\n");
76}
77
78__isl_give isl_factorizer *isl_factorizer_identity(__isl_keep isl_basic_set *bset)
79{
80	return isl_factorizer_alloc(isl_morph_identity(bset), 0);
81}
82
83__isl_give isl_factorizer *isl_factorizer_groups(__isl_keep isl_basic_set *bset,
84	__isl_take isl_mat *Q, __isl_take isl_mat *U, int n, int *len)
85{
86	int i;
87	unsigned nvar;
88	unsigned ovar;
89	isl_space *dim;
90	isl_basic_set *dom;
91	isl_basic_set *ran;
92	isl_morph *morph;
93	isl_factorizer *f;
94	isl_mat *id;
95
96	if (!bset || !Q || !U)
97		goto error;
98
99	ovar = 1 + isl_space_offset(bset->dim, isl_dim_set);
100	id = isl_mat_identity(bset->ctx, ovar);
101	Q = isl_mat_diagonal(isl_mat_copy(id), Q);
102	U = isl_mat_diagonal(id, U);
103
104	nvar = isl_basic_set_dim(bset, isl_dim_set);
105	dim = isl_basic_set_get_space(bset);
106	dom = isl_basic_set_universe(isl_space_copy(dim));
107	dim = isl_space_drop_dims(dim, isl_dim_set, 0, nvar);
108	dim = isl_space_add_dims(dim, isl_dim_set, nvar);
109	ran = isl_basic_set_universe(dim);
110	morph = isl_morph_alloc(dom, ran, Q, U);
111	f = isl_factorizer_alloc(morph, n);
112	if (!f)
113		return NULL;
114	for (i = 0; i < n; ++i)
115		f->len[i] = len[i];
116	return f;
117error:
118	isl_mat_free(Q);
119	isl_mat_free(U);
120	return NULL;
121}
122
123struct isl_factor_groups {
124	int *pos;		/* for each column: row position of pivot */
125	int *group;		/* group to which a column belongs */
126	int *cnt;		/* number of columns in the group */
127	int *rowgroup;		/* group to which a constraint belongs */
128};
129
130/* Initialize isl_factor_groups structure: find pivot row positions,
131 * each column initially belongs to its own group and the groups
132 * of the constraints are still unknown.
133 */
134static int init_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
135{
136	int i, j;
137
138	if (!H)
139		return -1;
140
141	g->pos = isl_alloc_array(H->ctx, int, H->n_col);
142	g->group = isl_alloc_array(H->ctx, int, H->n_col);
143	g->cnt = isl_alloc_array(H->ctx, int, H->n_col);
144	g->rowgroup = isl_alloc_array(H->ctx, int, H->n_row);
145
146	if (!g->pos || !g->group || !g->cnt || !g->rowgroup)
147		return -1;
148
149	for (i = 0; i < H->n_row; ++i)
150		g->rowgroup[i] = -1;
151	for (i = 0, j = 0; i < H->n_col; ++i) {
152		for ( ; j < H->n_row; ++j)
153			if (!isl_int_is_zero(H->row[j][i]))
154				break;
155		g->pos[i] = j;
156	}
157	for (i = 0; i < H->n_col; ++i) {
158		g->group[i] = i;
159		g->cnt[i] = 1;
160	}
161
162	return 0;
163}
164
165/* Update group[k] to the group column k belongs to.
166 * When merging two groups, only the group of the current
167 * group leader is changed.  Here we change the group of
168 * the other members to also point to the group that the
169 * old group leader now points to.
170 */
171static void update_group(struct isl_factor_groups *g, int k)
172{
173	int p = g->group[k];
174	while (g->cnt[p] == 0)
175		p = g->group[p];
176	g->group[k] = p;
177}
178
179/* Merge group i with all groups of the subsequent columns
180 * with non-zero coefficients in row j of H.
181 * (The previous columns are all zero; otherwise we would have handled
182 * the row before.)
183 */
184static int update_group_i_with_row_j(struct isl_factor_groups *g, int i, int j,
185	__isl_keep isl_mat *H)
186{
187	int k;
188
189	g->rowgroup[j] = g->group[i];
190	for (k = i + 1; k < H->n_col && j >= g->pos[k]; ++k) {
191		update_group(g, k);
192		update_group(g, i);
193		if (g->group[k] != g->group[i] &&
194		    !isl_int_is_zero(H->row[j][k])) {
195			isl_assert(H->ctx, g->cnt[g->group[k]] != 0, return -1);
196			isl_assert(H->ctx, g->cnt[g->group[i]] != 0, return -1);
197			if (g->group[i] < g->group[k]) {
198				g->cnt[g->group[i]] += g->cnt[g->group[k]];
199				g->cnt[g->group[k]] = 0;
200				g->group[g->group[k]] = g->group[i];
201			} else {
202				g->cnt[g->group[k]] += g->cnt[g->group[i]];
203				g->cnt[g->group[i]] = 0;
204				g->group[g->group[i]] = g->group[k];
205			}
206		}
207	}
208
209	return 0;
210}
211
212/* Update the group information based on the constraint matrix.
213 */
214static int update_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
215{
216	int i, j;
217
218	for (i = 0; i < H->n_col && g->cnt[0] < H->n_col; ++i) {
219		if (g->pos[i] == H->n_row)
220			continue; /* A line direction */
221		if (g->rowgroup[g->pos[i]] == -1)
222			g->rowgroup[g->pos[i]] = i;
223		for (j = g->pos[i] + 1; j < H->n_row; ++j) {
224			if (isl_int_is_zero(H->row[j][i]))
225				continue;
226			if (g->rowgroup[j] != -1)
227				continue;
228			if (update_group_i_with_row_j(g, i, j, H) < 0)
229				return -1;
230		}
231	}
232	for (i = 1; i < H->n_col; ++i)
233		update_group(g, i);
234
235	return 0;
236}
237
238static void clear_groups(struct isl_factor_groups *g)
239{
240	if (!g)
241		return;
242	free(g->pos);
243	free(g->group);
244	free(g->cnt);
245	free(g->rowgroup);
246}
247
248/* Determine if the set variables of the basic set can be factorized and
249 * return the results in an isl_factorizer.
250 *
251 * The algorithm works by first computing the Hermite normal form
252 * and then grouping columns linked by one or more constraints together,
253 * where a constraints "links" two or more columns if the constraint
254 * has nonzero coefficients in the columns.
255 */
256__isl_give isl_factorizer *isl_basic_set_factorizer(
257	__isl_keep isl_basic_set *bset)
258{
259	int i, j, n, done;
260	isl_mat *H, *U, *Q;
261	unsigned nvar;
262	struct isl_factor_groups g = { 0 };
263	isl_factorizer *f;
264
265	if (!bset)
266		return NULL;
267
268	isl_assert(bset->ctx, isl_basic_set_dim(bset, isl_dim_div) == 0,
269		return NULL);
270
271	nvar = isl_basic_set_dim(bset, isl_dim_set);
272	if (nvar <= 1)
273		return isl_factorizer_identity(bset);
274
275	H = isl_mat_alloc(bset->ctx, bset->n_eq + bset->n_ineq, nvar);
276	if (!H)
277		return NULL;
278	isl_mat_sub_copy(bset->ctx, H->row, bset->eq, bset->n_eq,
279		0, 1 + isl_space_offset(bset->dim, isl_dim_set), nvar);
280	isl_mat_sub_copy(bset->ctx, H->row + bset->n_eq, bset->ineq, bset->n_ineq,
281		0, 1 + isl_space_offset(bset->dim, isl_dim_set), nvar);
282	H = isl_mat_left_hermite(H, 0, &U, &Q);
283
284	if (init_groups(&g, H) < 0)
285		goto error;
286	if (update_groups(&g, H) < 0)
287		goto error;
288
289	if (g.cnt[0] == nvar) {
290		isl_mat_free(H);
291		isl_mat_free(U);
292		isl_mat_free(Q);
293		clear_groups(&g);
294
295		return isl_factorizer_identity(bset);
296	}
297
298	done = 0;
299	n = 0;
300	while (done != nvar) {
301		int group = g.group[done];
302		for (i = 1; i < g.cnt[group]; ++i) {
303			if (g.group[done + i] == group)
304				continue;
305			for (j = done + g.cnt[group]; j < nvar; ++j)
306				if (g.group[j] == group)
307					break;
308			if (j == nvar)
309				isl_die(bset->ctx, isl_error_internal,
310					"internal error", goto error);
311			g.group[j] = g.group[done + i];
312			Q = isl_mat_swap_rows(Q, done + i, j);
313			U = isl_mat_swap_cols(U, done + i, j);
314		}
315		done += g.cnt[group];
316		g.pos[n++] = g.cnt[group];
317	}
318
319	f = isl_factorizer_groups(bset, Q, U, n, g.pos);
320
321	isl_mat_free(H);
322	clear_groups(&g);
323
324	return f;
325error:
326	isl_mat_free(H);
327	isl_mat_free(U);
328	isl_mat_free(Q);
329	clear_groups(&g);
330	return NULL;
331}
332