Home | History | Annotate | Line # | Download | only in dist
      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  */
     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 
     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 
     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 
     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 
     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 
     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, off;
    101 	isl_space *space;
    102 	isl_basic_set *dom;
    103 	isl_basic_set *ran;
    104 	isl_morph *morph;
    105 	isl_factorizer *f;
    106 	isl_mat *id;
    107 
    108 	nvar = isl_basic_set_dim(bset, isl_dim_set);
    109 	off = isl_basic_set_var_offset(bset, isl_dim_set);
    110 	if (nvar < 0 || off < 0 || !Q || !U)
    111 		goto error;
    112 
    113 	id = isl_mat_identity(bset->ctx, 1 + off);
    114 	Q = isl_mat_diagonal(isl_mat_copy(id), Q);
    115 	U = isl_mat_diagonal(id, U);
    116 
    117 	space = isl_basic_set_get_space(bset);
    118 	dom = isl_basic_set_universe(isl_space_copy(space));
    119 	space = isl_space_drop_dims(space, isl_dim_set, 0, nvar);
    120 	space = isl_space_add_dims(space, isl_dim_set, nvar);
    121 	ran = isl_basic_set_universe(space);
    122 	morph = isl_morph_alloc(dom, ran, Q, U);
    123 	f = isl_factorizer_alloc(bset, morph, n);
    124 	if (!f)
    125 		return NULL;
    126 	for (i = 0; i < n; ++i)
    127 		f->len[i] = len[i];
    128 	return f;
    129 error:
    130 	isl_mat_free(Q);
    131 	isl_mat_free(U);
    132 	return NULL;
    133 }
    134 
    135 struct isl_factor_groups {
    136 	int *pos;		/* for each column: row position of pivot */
    137 	int *group;		/* group to which a column belongs */
    138 	int *cnt;		/* number of columns in the group */
    139 	int *rowgroup;		/* group to which a constraint belongs */
    140 };
    141 
    142 /* Initialize isl_factor_groups structure: find pivot row positions,
    143  * each column initially belongs to its own group and the groups
    144  * of the constraints are still unknown.
    145  */
    146 static int init_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
    147 {
    148 	int i, j;
    149 
    150 	if (!H)
    151 		return -1;
    152 
    153 	g->pos = isl_alloc_array(H->ctx, int, H->n_col);
    154 	g->group = isl_alloc_array(H->ctx, int, H->n_col);
    155 	g->cnt = isl_alloc_array(H->ctx, int, H->n_col);
    156 	g->rowgroup = isl_alloc_array(H->ctx, int, H->n_row);
    157 
    158 	if (!g->pos || !g->group || !g->cnt || !g->rowgroup)
    159 		return -1;
    160 
    161 	for (i = 0; i < H->n_row; ++i)
    162 		g->rowgroup[i] = -1;
    163 	for (i = 0, j = 0; i < H->n_col; ++i) {
    164 		for ( ; j < H->n_row; ++j)
    165 			if (!isl_int_is_zero(H->row[j][i]))
    166 				break;
    167 		g->pos[i] = j;
    168 	}
    169 	for (i = 0; i < H->n_col; ++i) {
    170 		g->group[i] = i;
    171 		g->cnt[i] = 1;
    172 	}
    173 
    174 	return 0;
    175 }
    176 
    177 /* Update group[k] to the group column k belongs to.
    178  * When merging two groups, only the group of the current
    179  * group leader is changed.  Here we change the group of
    180  * the other members to also point to the group that the
    181  * old group leader now points to.
    182  */
    183 static void update_group(struct isl_factor_groups *g, int k)
    184 {
    185 	int p = g->group[k];
    186 	while (g->cnt[p] == 0)
    187 		p = g->group[p];
    188 	g->group[k] = p;
    189 }
    190 
    191 /* Merge group i with all groups of the subsequent columns
    192  * with non-zero coefficients in row j of H.
    193  * (The previous columns are all zero; otherwise we would have handled
    194  * the row before.)
    195  */
    196 static int update_group_i_with_row_j(struct isl_factor_groups *g, int i, int j,
    197 	__isl_keep isl_mat *H)
    198 {
    199 	int k;
    200 
    201 	g->rowgroup[j] = g->group[i];
    202 	for (k = i + 1; k < H->n_col && j >= g->pos[k]; ++k) {
    203 		update_group(g, k);
    204 		update_group(g, i);
    205 		if (g->group[k] != g->group[i] &&
    206 		    !isl_int_is_zero(H->row[j][k])) {
    207 			isl_assert(H->ctx, g->cnt[g->group[k]] != 0, return -1);
    208 			isl_assert(H->ctx, g->cnt[g->group[i]] != 0, return -1);
    209 			if (g->group[i] < g->group[k]) {
    210 				g->cnt[g->group[i]] += g->cnt[g->group[k]];
    211 				g->cnt[g->group[k]] = 0;
    212 				g->group[g->group[k]] = g->group[i];
    213 			} else {
    214 				g->cnt[g->group[k]] += g->cnt[g->group[i]];
    215 				g->cnt[g->group[i]] = 0;
    216 				g->group[g->group[i]] = g->group[k];
    217 			}
    218 		}
    219 	}
    220 
    221 	return 0;
    222 }
    223 
    224 /* Update the group information based on the constraint matrix.
    225  */
    226 static int update_groups(struct isl_factor_groups *g, __isl_keep isl_mat *H)
    227 {
    228 	int i, j;
    229 
    230 	for (i = 0; i < H->n_col && g->cnt[0] < H->n_col; ++i) {
    231 		if (g->pos[i] == H->n_row)
    232 			continue; /* A line direction */
    233 		if (g->rowgroup[g->pos[i]] == -1)
    234 			g->rowgroup[g->pos[i]] = i;
    235 		for (j = g->pos[i] + 1; j < H->n_row; ++j) {
    236 			if (isl_int_is_zero(H->row[j][i]))
    237 				continue;
    238 			if (g->rowgroup[j] != -1)
    239 				continue;
    240 			if (update_group_i_with_row_j(g, i, j, H) < 0)
    241 				return -1;
    242 		}
    243 	}
    244 	for (i = 1; i < H->n_col; ++i)
    245 		update_group(g, i);
    246 
    247 	return 0;
    248 }
    249 
    250 static void clear_groups(struct isl_factor_groups *g)
    251 {
    252 	if (!g)
    253 		return;
    254 	free(g->pos);
    255 	free(g->group);
    256 	free(g->cnt);
    257 	free(g->rowgroup);
    258 }
    259 
    260 /* Determine if the set variables of the basic set can be factorized and
    261  * return the results in an isl_factorizer.
    262  *
    263  * The algorithm works by first computing the Hermite normal form
    264  * and then grouping columns linked by one or more constraints together,
    265  * where a constraints "links" two or more columns if the constraint
    266  * has nonzero coefficients in the columns.
    267  */
    268 __isl_give isl_factorizer *isl_basic_set_factorizer(
    269 	__isl_keep isl_basic_set *bset)
    270 {
    271 	int i, j, n, done;
    272 	isl_mat *H, *U, *Q;
    273 	isl_size nvar, first;
    274 	struct isl_factor_groups g = { 0 };
    275 	isl_factorizer *f;
    276 
    277 	nvar = isl_basic_set_dim(bset, isl_dim_set);
    278 	first = isl_basic_set_var_offset(bset, isl_dim_set);
    279 	if (nvar < 0 || first < 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 + first, nvar);
    290 	isl_mat_sub_copy(bset->ctx, H->row + bset->n_eq, bset->ineq, bset->n_ineq,
    291 		0, 1 + first, 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  */
    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