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 (Z^{n+1} -> Z^{n+1})
  76  * that maps an element x to any element that can be reached
  77  * by taking a non-negative number of steps along any of
  78  * the extended offsets v'_i = [v_i 1].
  79  * That is, construct
  80  *
  81  * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
  82  *
  83  * For any element in this relation, the number of steps taken
  84  * is equal to the difference in the final coordinates.
  85  */
  86 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
  87         __isl_keep isl_mat *steps)
  88 {
  89         int i, j, k;
  90         struct isl_basic_map *path = NULL;
  91         unsigned d;
  92         unsigned n;
  93         unsigned nparam;
  94
  95         if (!dim || !steps)
  96                 goto error;
  97
  98         d = isl_dim_size(dim, isl_dim_in);
  99         n = steps->n_row;
 100         nparam = isl_dim_size(dim, isl_dim_param);
 101
 102         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d, n);
 103
 104         for (i = 0; i < n; ++i) {
 105                 k = isl_basic_map_alloc_div(path);
 106                 if (k < 0)
 107                         goto error;
 108                 isl_assert(steps->ctx, i == k, goto error);
 109                 isl_int_set_si(path->div[k][0], 0);
 110         }
 111
 112         for (i = 0; i < d; ++i) {
 113                 k = isl_basic_map_alloc_equality(path);
 114                 if (k < 0)
 115                         goto error;
 116                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 117                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 118                 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
 119                 if (i == d - 1)
 120                         for (j = 0; j < n; ++j)
 121                                 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
 122                 else
 123                         for (j = 0; j < n; ++j)
 124                                 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
 125                                             steps->row[j][i]);
 126         }
 127
 128         for (i = 0; i < n; ++i) {
 129                 k = isl_basic_map_alloc_inequality(path);
 130                 if (k < 0)
 131                         goto error;
 132                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 133                 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
 134         }
 135
 136         isl_dim_free(dim);
 137
 138         path = isl_basic_map_simplify(path);
 139         path = isl_basic_map_finalize(path);
 140         return isl_map_from_basic_map(path);
 141 error:
 142         isl_dim_free(dim);
 143         isl_basic_map_free(path);
 144         return NULL;
 145 }
 146
 147 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
 148  * construct a map that equates the parameter to the difference
 149  * in the final coordinates and imposes that this difference is positive.
 150  * That is, construct
 151  *
 152  *      { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
 153  */
 154 static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_dim *dim,
 155         unsigned param)
 156 {
 157         struct isl_basic_map *bmap;
 158         unsigned d;
 159         unsigned nparam;
 160         int k;
 161
 162         d = isl_dim_size(dim, isl_dim_in);
 163         nparam = isl_dim_size(dim, isl_dim_param);
 164         bmap = isl_basic_map_alloc_dim(dim, 0, 1, 1);
 165         k = isl_basic_map_alloc_equality(bmap);
 166         if (k < 0)
 167                 goto error;
 168         isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
 169         isl_int_set_si(bmap->eq[k][1 + param], -1);
 170         isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
 171         isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
 172
 173         k = isl_basic_map_alloc_inequality(bmap);
 174         if (k < 0)
 175                 goto error;
 176         isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
 177         isl_int_set_si(bmap->ineq[k][1 + param], 1);
 178         isl_int_set_si(bmap->ineq[k][0], -1);
 179
 180         bmap = isl_basic_map_finalize(bmap);
 181         return isl_map_from_basic_map(bmap);
 182 error:
 183         isl_basic_map_free(bmap);
 184         return NULL;
 185 }
 186
 187 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
 188  * construct a map that is an overapproximation of the map
 189  * that takes an element from the space D to another
 190  * element from the same space, such that the difference between
 191  * them is a strictly positive sum of differences between images
 192  * and pre-images in one of the R_i.
 193  * The number of differences in the sum is equated to parameter "param".
 194  * That is, let
 195  *
 196  *      \Delta_i = { y - x | (x, y) in R_i }
 197  *
 198  * then the constructed map is an overapproximation of
 199  *
 200  *      { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
 201  *                              d = \sum_i k_i and k = \sum_i k_i > 0 }
 202  *
 203  * We first construct an extended mapping with an extra coordinate
 204  * that indicates the number of steps taken.  In particular,
 205  * the difference in the last coordinate is equal to the number
 206  * of steps taken to move from a domain element to the corresponding
 207  * image element(s).
 208  * In the final step, this difference is equated to the parameter "param"
 209  * and made positive.  The extra coordinates are subsequently projected out.
 210  *
 211  * If any of the \Delta_i contains more than one element, then
 212  * we currently simply return a universal map { (x) -> (y) }.
 213  */
 214 static __isl_give isl_map *construct_path(__isl_keep isl_map *map,
 215         unsigned param)
 216 {
 217         struct isl_mat *steps = NULL;
 218         struct isl_map *path = NULL;
 219         struct isl_map *diff;
 220         struct isl_dim *dim;
 221         unsigned d;
 222         int ok;
 223
 224         if (!map)
 225                 return NULL;
 226
 227         steps = extract_steps(map, &ok);
 228         if (!steps)
 229                 return NULL;
 230         if (!ok) {
 231                 isl_mat_free(steps);
 232                 return isl_map_universe(isl_map_get_dim(map));
 233         }
 234
 235         dim = isl_map_get_dim(map);
 236
 237         d = isl_dim_size(dim, isl_dim_in);
 238         dim = isl_dim_add(dim, isl_dim_in, 1);
 239         dim = isl_dim_add(dim, isl_dim_out, 1);
 240         path = path_along_steps(isl_dim_copy(dim), steps);
 241         diff = equate_parameter_to_length(dim, param);
 242
 243         path = isl_map_intersect(path, diff);
 244         path = isl_map_project_out(path, isl_dim_in, d, 1);
 245         path = isl_map_project_out(path, isl_dim_out, d, 1);
 246
 247         isl_mat_free(steps);
 248         return path;
 249 }
 250
 251 /* Check whether "path" is acyclic.
 252  * That is, check whether
 253  *
 254  *      { d | d = y - x and (x,y) in path }
 255  *
 256  * does not contain the origin.
 257  */
 258 static int is_acyclic(__isl_take isl_map *path)
 259 {
 260         int i;
 261         int acyclic;
 262         unsigned dim;
 263         struct isl_set *delta;
 264
 265         delta = isl_map_deltas(path);
 266         dim = isl_set_dim(delta, isl_dim_set);
 267         for (i = 0; i < dim; ++i)
 268                 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
 269
 270         acyclic = isl_set_is_empty(delta);
 271         isl_set_free(delta);
 272
 273         return acyclic;
 274 }
 275
 276 /* Shift variable at position "pos" up by one.
 277  * That is, replace the corresponding variable v by v - 1.
 278  */
 279 static __isl_give isl_basic_map *basic_map_shift_pos(
 280         __isl_take isl_basic_map *bmap, unsigned pos)
 281 {
 282         int i;
 283
 284         bmap = isl_basic_map_cow(bmap);
 285         if (!bmap)
 286                 return NULL;
 287
 288         for (i = 0; i < bmap->n_eq; ++i)
 289                 isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
 290
 291         for (i = 0; i < bmap->n_ineq; ++i)
 292                 isl_int_sub(bmap->ineq[i][0],
 293                                 bmap->ineq[i][0], bmap->ineq[i][pos]);
 294
 295         for (i = 0; i < bmap->n_div; ++i) {
 296                 if (isl_int_is_zero(bmap->div[i][0]))
 297                         continue;
 298                 isl_int_sub(bmap->div[i][1],
 299                                 bmap->div[i][1], bmap->div[i][1 + pos]);
 300         }
 301
 302         return bmap;
 303 }
 304
 305 /* Shift variable at position "pos" up by one.
 306  * That is, replace the corresponding variable v by v - 1.
 307  */
 308 static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
 309 {
 310         int i;
 311
 312         map = isl_map_cow(map);
 313         if (!map)
 314                 return NULL;
 315
 316         for (i = 0; i < map->n; ++i) {
 317                 map->p[i] = basic_map_shift_pos(map->p[i], pos);
 318                 if (!map->p[i])
 319                         goto error;
 320         }
 321         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
 322         return map;
 323 error:
 324         isl_map_free(map);
 325         return NULL;
 326 }
 327
 328 /* Check whether the overapproximation of the power of "map" is exactly
 329  * the power of "map".  Let R be "map" and A_k the overapproximation.
 330  * The approximation is exact if
 331  *
 332  *      A_1 = R
 333  *      A_k = A_{k-1} \circ R                   k >= 2
 334  *
 335  * Since A_k is known to be an overapproximation, we only need to check
 336  *
 337  *      A_1 \subset R
 338  *      A_k \subset A_{k-1} \circ R             k >= 2
 339  *
 340  */
 341 static int check_power_exactness(__isl_take isl_map *map,
 342         __isl_take isl_map *app, unsigned param)
 343 {
 344         int exact;
 345         isl_map *app_1;
 346         isl_map *app_2;
 347
 348         app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
 349
 350         exact = isl_map_is_subset(app_1, map);
 351         isl_map_free(app_1);
 352
 353         if (!exact || exact < 0) {
 354                 isl_map_free(app);
 355                 isl_map_free(map);
 356                 return exact;
 357         }
 358
 359         app_2 = isl_map_lower_bound_si(isl_map_copy(app),
 360                                         isl_dim_param, param, 2);
 361         app_1 = map_shift_pos(app, 1 + param);
 362         app_1 = isl_map_apply_range(map, app_1);
 363
 364         exact = isl_map_is_subset(app_2, app_1);
 365
 366         isl_map_free(app_1);
 367         isl_map_free(app_2);
 368
 369         return exact;
 370 }
 371
 372 /* Check whether the overapproximation of the power of "map" is exactly
 373  * the power of "map", possibly after projecting out the power (if "project"
 374  * is set).
 375  *
 376  * If "project" is set and if "steps" can only result in acyclic paths,
 377  * then we check
 378  *
 379  *      A = R \cup (A \circ R)
 380  *
 381  * where A is the overapproximation with the power projected out, i.e.,
 382  * an overapproximation of the transitive closure.
 383  * More specifically, since A is known to be an overapproximation, we check
 384  *
 385  *      A \subset R \cup (A \circ R)
 386  *
 387  * Otherwise, we check if the power is exact.
 388  */
 389 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
 390         __isl_take isl_map *path, unsigned param, int project)
 391 {
 392         isl_map *test;
 393         int exact;
 394
 395         if (project) {
 396                 project = is_acyclic(path);
 397                 if (project < 0)
 398                         goto error;
 399         } else
 400                 isl_map_free(path);
 401
 402         if (!project)
 403                 return check_power_exactness(map, app, param);
 404
 405         map = isl_map_project_out(map, isl_dim_param, param, 1);
 406         app = isl_map_project_out(app, isl_dim_param, param, 1);
 407
 408         test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
 409         test = isl_map_union(test, isl_map_copy(map));
 410
 411         exact = isl_map_is_subset(app, test);
 412
 413         isl_map_free(app);
 414         isl_map_free(test);
 415
 416         isl_map_free(map);
 417
 418         return exact;
 419 error:
 420         isl_map_free(app);
 421         isl_map_free(map);
 422         return -1;
 423 }
 424
 425 /* Compute the positive powers of "map", or an overapproximation.
 426  * The power is given by parameter "param".  If the result is exact,
 427  * then *exact is set to 1.
 428  * If project is set, then we are actually interested in the transitive
 429  * closure, so we can use a more relaxed exactness check.
 430  */
 431 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 432         int *exact, int project)
 433 {
 434         struct isl_set *domain = NULL;
 435         struct isl_set *range = NULL;
 436         struct isl_map *app = NULL;
 437         struct isl_map *path = NULL;
 438
 439         if (exact)
 440                 *exact = 1;
 441
 442         map = isl_map_remove_empty_parts(map);
 443         if (!map)
 444                 return NULL;
 445
 446         if (isl_map_fast_is_empty(map))
 447                 return map;
 448
 449         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
 450         isl_assert(map->ctx,
 451                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 452                 goto error);
 453
 454         domain = isl_map_domain(isl_map_copy(map));
 455         domain = isl_set_coalesce(domain);
 456         range = isl_map_range(isl_map_copy(map));
 457         range = isl_set_coalesce(range);
 458         app = isl_map_from_domain_and_range(isl_set_copy(domain),
 459                                             isl_set_copy(range));
 460
 461         path = construct_path(map, param);
 462         app = isl_map_intersect(app, isl_map_copy(path));
 463
 464         if (exact &&
 465             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 466                                   isl_map_copy(path), param, project)) < 0)
 467                 goto error;
 468
 469         isl_set_free(domain);
 470         isl_set_free(range);
 471         isl_map_free(path);
 472         isl_map_free(map);
 473         return app;
 474 error:
 475         isl_set_free(domain);
 476         isl_set_free(range);
 477         isl_map_free(path);
 478         isl_map_free(map);
 479         isl_map_free(app);
 480         return NULL;
 481 }
 482
 483 /* Compute the positive powers of "map", or an overapproximation.
 484  * The power is given by parameter "param".  If the result is exact,
 485  * then *exact is set to 1.
 486  */
 487 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
 488         int *exact)
 489 {
 490         return map_power(map, param, exact, 0);
 491 }
 492
 493 /* Compute the transitive closure  of "map", or an overapproximation.
 494  * If the result is exact, then *exact is set to 1.
 495  * Simply compute the powers of map and then project out the parameter
 496  * describing the power.
 497  */
 498 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
 499         int *exact)
 500 {
 501         unsigned param;
 502
 503         if (!map)
 504                 goto error;
 505
 506         param = isl_map_dim(map, isl_dim_param);
 507         map = isl_map_add(map, isl_dim_param, 1);
 508         map = map_power(map, param, exact, 1);
 509         map = isl_map_project_out(map, isl_dim_param, param, 1);
 510
 511         return map;
 512 error:
 513         isl_map_free(map);
 514         return NULL;
 515 }