390 lines · c
1/*2 * Copyright 2005-2007 Universiteit Leiden3 * Copyright 2008-2009 Katholieke Universiteit Leuven4 * Copyright 2010 INRIA Saclay5 *6 * Use of this software is governed by the MIT license7 *8 * Written by Sven Verdoolaege, Leiden Institute of Advanced Computer Science,9 * Universiteit Leiden, Niels Bohrweg 1, 2333 CA Leiden, The Netherlands10 * and K.U.Leuven, Departement Computerwetenschappen, Celestijnenlaan 200A,11 * B-3001 Leuven, Belgium12 * 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 */23isl_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 30static __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;56error: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 74void 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;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;130error:131 isl_mat_free(Q);132 isl_mat_free(U);133 return NULL;134}135 136struct 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 groups145 * of the constraints are still unknown.146 */147static 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 current180 * group leader is changed. Here we change the group of181 * the other members to also point to the group that the182 * old group leader now points to.183 */184static 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 columns193 * with non-zero coefficients in row j of H.194 * (The previous columns are all zero; otherwise we would have handled195 * the row before.)196 */197static 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 */227static 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 251static 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 and262 * return the results in an isl_factorizer.263 *264 * The algorithm works by first computing the Hermite normal form265 * and then grouping columns linked by one or more constraints together,266 * where a constraints "links" two or more columns if the constraint267 * has nonzero coefficients in the columns.268 */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;335error: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