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 GNU LGPLv2.1 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_dim_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, FILE *out)
  62 {
  63         int i;
  64
  65         if (!f)
  66                 return;
  67
  68         isl_morph_dump(f->morph, out);
  69         fprintf(out, "[");
  70         for (i = 0; i < f->n_group; ++i) {
  71                 if (i)
  72                         fprintf(out, ", ");
  73                 fprintf(out, "%d", f->len[i]);
  74         }
  75         fprintf(out, "]\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_dim *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_dim_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_dim(bset);
 106         dom = isl_basic_set_universe(isl_dim_copy(dim));
 107         dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
 108         dim = isl_dim_add(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_dim_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_dim_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 }