isl_factorization.c

   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 }