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 static __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; 45 error: 46 free(len); 47 isl_morph_free(morph); 48 return NULL; 49 } 50 51 void 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 61 void 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; 117 error: 118 isl_mat_free(Q); 119 isl_mat_free(U); 120 return NULL; 121 } 122 123 struct 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 */ 134 static 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 */ 171 static 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 */ 184 static 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 */ 214 static 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 238 static 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; 325 error: 326 isl_mat_free(H); 327 isl_mat_free(U); 328 isl_mat_free(Q); 329 clear_groups(&g); 330 return NULL; 331 } 332