isl_transitive_closure.c

   1 /*
   2  * Copyright 2010      INRIA Saclay
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
   7  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
   8  * 91893 Orsay, France
   9  */
  10
  11 #include "isl_map.h"
  12 #include "isl_map_private.h"
  13
  14 /*
  15  * The transitive closure implementation is based on the paper
  16  * "Computing the Transitive Closure of a Union of Affine Integer
  17  * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
  18  * Albert Cohen.
  19  */
  20
  21 /* Given a union of translations (uniform dependences), return a matrix
  22  * with as many rows as there are disjuncts in the union and as many
  23  * columns as there are input dimensions (which should be equal to
  24  * the number of output dimensions).
  25  * Each row contains the translation performed by the corresponding disjunct.
  26  * If "map" turns out not to be a union of translations, then the contents
  27  * of the returned matrix are undefined and *ok is set to 0.
  28  */
  29 static __isl_give isl_mat *extract_steps(__isl_keep isl_map *map, int *ok)
  30 {
  31         int i, j;
  32         struct isl_mat *steps;
  33         unsigned dim = isl_map_dim(map, isl_dim_in);
  34
  35         *ok = 1;
  36
  37         steps = isl_mat_alloc(map->ctx, map->n, dim);
  38         if (!steps)
  39                 return NULL;
  40
  41         for (i = 0; i < map->n; ++i) {
  42                 struct isl_basic_set *delta;
  43
  44                 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
  45
  46                 for (j = 0; j < dim; ++j) {
  47                         int fixed;
  48
  49                         fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
  50                                                             &steps->row[i][j]);
  51                         if (fixed < 0) {
  52                                 isl_basic_set_free(delta);
  53                                 goto error;
  54                         }
  55                         if (!fixed)
  56                                 break;
  57                 }
  58
  59                 isl_basic_set_free(delta);
  60
  61                 if (j < dim)
  62                         break;
  63         }
  64
  65         if (i < map->n)
  66                 *ok = 0;
  67
  68         return steps;
  69 error:
  70         isl_mat_free(steps);
  71         return NULL;
  72 }
  73
  74 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
  75  * of the given dimension specification that maps a element x to any
  76  * element that can be reached by taking a positive number of steps
  77  * along any of the offsets, where the number of steps k is equal to
  78  * parameter "param".  That is, construct
  79  *
  80  * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v_i, k = \sum_i k_i > 0 }
  81  *
  82  * If strict is non-negative, then at least one step should be taken
  83  * along the given offset v_strict, i.e., k_strict > 0.
  84  */
  85 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
  86         __isl_keep isl_mat *steps, unsigned param, int strict)
  87 {
  88         int i, j, k;
  89         struct isl_basic_map *path = NULL;
  90         unsigned d;
  91         unsigned n;
  92         unsigned nparam;
  93
  94         if (!dim || !steps)
  95                 goto error;
  96
  97         d = isl_dim_size(dim, isl_dim_in);
  98         n = steps->n_row;
  99         nparam = isl_dim_size(dim, isl_dim_param);
 100
 101         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d + 1, n + 1);
 102
 103         for (i = 0; i < n; ++i) {
 104                 k = isl_basic_map_alloc_div(path);
 105                 if (k < 0)
 106                         goto error;
 107                 isl_assert(steps->ctx, i == k, goto error);
 108                 isl_int_set_si(path->div[k][0], 0);
 109         }
 110
 111         for (i = 0; i < d; ++i) {
 112                 k = isl_basic_map_alloc_equality(path);
 113                 if (k < 0)
 114                         goto error;
 115                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 116                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 117                 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
 118                 for (j = 0; j < n; ++j)
 119                         isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
 120                                     steps->row[j][i]);
 121         }
 122
 123         k = isl_basic_map_alloc_equality(path);
 124         if (k < 0)
 125                 goto error;
 126         isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 127         isl_int_set_si(path->eq[k][1 + param], -1);
 128         for (j = 0; j < n; ++j)
 129                 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
 130
 131         for (i = 0; i < n; ++i) {
 132                 k = isl_basic_map_alloc_inequality(path);
 133                 if (k < 0)
 134                         goto error;
 135                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 136                 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
 137                 if (i == strict)
 138                         isl_int_set_si(path->ineq[k][0], -1);
 139         }
 140
 141         k = isl_basic_map_alloc_inequality(path);
 142         if (k < 0)
 143                 goto error;
 144         isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 145         isl_int_set_si(path->ineq[k][1 + param], 1);
 146         isl_int_set_si(path->ineq[k][0], -1);
 147
 148         isl_dim_free(dim);
 149
 150         path = isl_basic_map_simplify(path);
 151         path = isl_basic_map_finalize(path);
 152         return isl_map_from_basic_map(path);
 153 error:
 154         isl_dim_free(dim);
 155         isl_basic_map_free(path);
 156         return NULL;
 157 }
 158
 159 /* Check whether the overapproximation of the power of "map" is exactly
 160  * the power of "map".  In particular, for each path of a given length
 161  * that starts of in domain or range and ends up in the range,
 162  * check whether there is at least one path of the same length
 163  * with a valid first segment, i.e., one in "map".
 164  * If "project" is set, then we drop the condition that the lengths
 165  * should be the same.
 166  *
 167  * "domain" and "range" are the domain and range of "map"
 168  * "steps" represents the translations of "map"
 169  * "path" is a path along "steps"
 170  *
 171  * "domain", "range", "steps" and "path" have been precomputed by the calling
 172  * function.
 173  */
 174 static int check_path_exactness(__isl_take isl_map *map,
 175         __isl_take isl_set *domain, __isl_take isl_set *range,
 176         __isl_take isl_map *path, __isl_keep isl_mat *steps, unsigned param,
 177         int project)
 178 {
 179         isl_map *test;
 180         int ok;
 181         int i;
 182
 183         if (!map)
 184                 goto error;
 185
 186         test = isl_map_empty(isl_map_get_dim(map));
 187         for (i = 0; i < map->n; ++i) {
 188                 struct isl_map *path_i;
 189                 struct isl_set *dom_i;
 190                 path_i = path_along_steps(isl_map_get_dim(map), steps, param, i);
 191                 dom_i = isl_set_from_basic_set(
 192                         isl_basic_map_domain(isl_basic_map_copy(map->p[i])));
 193                 path_i = isl_map_intersect_domain(path_i, dom_i);
 194                 test = isl_map_union(test, path_i);
 195         }
 196         isl_map_free(map);
 197         test = isl_map_intersect_range(test, isl_set_copy(range));
 198
 199         domain = isl_set_union(domain, isl_set_copy(range));
 200         path = isl_map_intersect_domain(path, domain);
 201         path = isl_map_intersect_range(path, range);
 202
 203         if (project) {
 204                 path = isl_map_project_out(path, isl_dim_param, param, 1);
 205                 test = isl_map_project_out(test, isl_dim_param, param, 1);
 206         }
 207
 208         ok = isl_map_is_subset(path, test);
 209
 210         isl_map_free(path);
 211         isl_map_free(test);
 212
 213         return ok;
 214 error:
 215         isl_map_free(map);
 216         isl_set_free(domain);
 217         isl_set_free(range);
 218         isl_map_free(path);
 219         return -1;
 220 }
 221
 222 /* Check whether any path of length at least one along "steps" is acyclic.
 223  * That is, check whether
 224  *
 225  *      \sum_i k_i \delta_i = 0
 226  *      \sum_i k_i >= 1
 227  *      k_i >= 0
 228  *
 229  * with \delta_i the rows of "steps", is infeasible.
 230  */
 231 static int is_acyclic(__isl_keep isl_mat *steps)
 232 {
 233         int acyclic;
 234         int i, j, k;
 235         struct isl_basic_set *test;
 236
 237         if (!steps)
 238                 return -1;
 239
 240         test = isl_basic_set_alloc(steps->ctx, 0, steps->n_row, 0,
 241                                         steps->n_col, steps->n_row + 1);
 242
 243         for (i = 0; i < steps->n_col; ++i) {
 244                 k = isl_basic_set_alloc_equality(test);
 245                 if (k < 0)
 246                         goto error;
 247                 isl_int_set_si(test->eq[k][0], 0);
 248                 for (j = 0; j < steps->n_row; ++j)
 249                         isl_int_set(test->eq[k][1 + j], steps->row[j][i]);
 250         }
 251         for (j = 0; j < steps->n_row; ++j) {
 252                 k = isl_basic_set_alloc_inequality(test);
 253                 if (k < 0)
 254                         goto error;
 255                 isl_seq_clr(test->ineq[k], 1 + steps->n_row);
 256                 isl_int_set_si(test->ineq[k][1 + j], 1);
 257         }
 258
 259         k = isl_basic_set_alloc_inequality(test);
 260         if (k < 0)
 261                 goto error;
 262         isl_int_set_si(test->ineq[k][0], -1);
 263         for (j = 0; j < steps->n_row; ++j)
 264                 isl_int_set_si(test->ineq[k][1 + j], 1);
 265
 266         acyclic = isl_basic_set_is_empty(test);
 267
 268         isl_basic_set_free(test);
 269         return acyclic;
 270 error:
 271         isl_basic_set_free(test);
 272         return -1;
 273 }
 274
 275 /* Check whether the overapproximation of the power of "map" is exactly
 276  * the power of "map", possibly after projecting out the power (if "project"
 277  * is set).
 278  *
 279  * If "project" is not set, then we simply check for each power if there
 280  * is a path of the given length with valid initial segment.
 281  * If "project" is set, then we check if "steps" can only result in acyclic
 282  * paths.  If so, we only need to check that there is a path of _some_
 283  * length >= 1.  Otherwise, we perform the standard check, i.e., whether
 284  * there is a path of the given length.
 285  */
 286 static int check_exactness(__isl_take isl_map *map, __isl_take isl_set *domain,
 287         __isl_take isl_set *range, __isl_take isl_map *path,
 288         __isl_keep isl_mat *steps, unsigned param, int project)
 289 {
 290         int acyclic;
 291
 292         if (!project)
 293                 return check_path_exactness(map, domain, range, path, steps,
 294                                                 param, 0);
 295
 296         acyclic = is_acyclic(steps);
 297         if (acyclic < 0)
 298                 goto error;
 299
 300         return check_path_exactness(map, domain, range, path, steps,
 301                                         param, acyclic);
 302 error:
 303         isl_map_free(map);
 304         isl_set_free(domain);
 305         isl_set_free(range);
 306         isl_map_free(path);
 307         return -1;
 308 }
 309
 310 /* Compute the positive powers of "map", or an overapproximation.
 311  * The power is given by parameter "param".  If the result is exact,
 312  * then *exact is set to 1.
 313  * If project is set, then we are actually interested in the transitive
 314  * closure, so we can use a more relaxed exactness check.
 315  */
 316 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 317         int *exact, int project)
 318 {
 319         struct isl_mat *steps = NULL;
 320         struct isl_set *domain = NULL;
 321         struct isl_set *range = NULL;
 322         struct isl_map *app = NULL;
 323         struct isl_map *path = NULL;
 324         int ok;
 325
 326         if (exact)
 327                 *exact = 1;
 328
 329         map = isl_map_remove_empty_parts(map);
 330         if (!map)
 331                 return NULL;
 332
 333         if (isl_map_fast_is_empty(map))
 334                 return map;
 335
 336         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
 337         isl_assert(map->ctx,
 338                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 339                 goto error);
 340
 341         domain = isl_map_domain(isl_map_copy(map));
 342         range = isl_map_range(isl_map_copy(map));
 343         app = isl_map_from_domain_and_range(isl_set_copy(domain),
 344                                             isl_set_copy(range));
 345
 346         steps = extract_steps(map, &ok);
 347         if (!ok)
 348                 goto not_exact;
 349
 350         path = path_along_steps(isl_map_get_dim(map), steps, param, -1);
 351         app = isl_map_intersect(app, isl_map_copy(path));
 352
 353         if (exact &&
 354             (*exact = check_exactness(isl_map_copy(map), isl_set_copy(domain),
 355                                   isl_set_copy(range), isl_map_copy(path),
 356                                   steps, param, project)) < 0)
 357                 goto error;
 358
 359         if (0) {
 360 not_exact:
 361                 if (exact)
 362                         *exact = 0;
 363         }
 364         isl_set_free(domain);
 365         isl_set_free(range);
 366         isl_mat_free(steps);
 367         isl_map_free(path);
 368         isl_map_free(map);
 369         return app;
 370 error:
 371         isl_set_free(domain);
 372         isl_set_free(range);
 373         isl_mat_free(steps);
 374         isl_map_free(path);
 375         isl_map_free(map);
 376         isl_map_free(app);
 377         return NULL;
 378 }
 379
 380 /* Compute the positive powers of "map", or an overapproximation.
 381  * The power is given by parameter "param".  If the result is exact,
 382  * then *exact is set to 1.
 383  */
 384 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
 385         int *exact)
 386 {
 387         return map_power(map, param, exact, 0);
 388 }
 389
 390 /* Compute the transitive closure  of "map", or an overapproximation.
 391  * If the result is exact, then *exact is set to 1.
 392  * Simply compute the powers of map and then project out the parameter
 393  * describing the power.
 394  */
 395 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
 396         int *exact)
 397 {
 398         unsigned param;
 399
 400         if (!map)
 401                 goto error;
 402
 403         param = isl_map_dim(map, isl_dim_param);
 404         map = isl_map_add(map, isl_dim_param, 1);
 405         map = map_power(map, param, exact, 1);
 406         map = isl_map_project_out(map, isl_dim_param, param, 1);
 407
 408         return map;
 409 error:
 410         isl_map_free(map);
 411         return NULL;
 412 }