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  */
  83 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
  84         __isl_keep isl_mat *steps, unsigned param)
  85 {
  86         int i, j, k;
  87         struct isl_basic_map *path = NULL;
  88         unsigned d;
  89         unsigned n;
  90         unsigned nparam;
  91
  92         if (!dim || !steps)
  93                 goto error;
  94
  95         d = isl_dim_size(dim, isl_dim_in);
  96         n = steps->n_row;
  97         nparam = isl_dim_size(dim, isl_dim_param);
  98
  99         path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d + 1, n + 1);
 100
 101         for (i = 0; i < n; ++i) {
 102                 k = isl_basic_map_alloc_div(path);
 103                 if (k < 0)
 104                         goto error;
 105                 isl_assert(steps->ctx, i == k, goto error);
 106                 isl_int_set_si(path->div[k][0], 0);
 107         }
 108
 109         for (i = 0; i < d; ++i) {
 110                 k = isl_basic_map_alloc_equality(path);
 111                 if (k < 0)
 112                         goto error;
 113                 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 114                 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
 115                 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
 116                 for (j = 0; j < n; ++j)
 117                         isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
 118                                     steps->row[j][i]);
 119         }
 120
 121         k = isl_basic_map_alloc_equality(path);
 122         if (k < 0)
 123                 goto error;
 124         isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
 125         isl_int_set_si(path->eq[k][1 + param], -1);
 126         for (j = 0; j < n; ++j)
 127                 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
 128
 129         for (i = 0; i < n; ++i) {
 130                 k = isl_basic_map_alloc_inequality(path);
 131                 if (k < 0)
 132                         goto error;
 133                 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 134                 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
 135         }
 136
 137         k = isl_basic_map_alloc_inequality(path);
 138         if (k < 0)
 139                 goto error;
 140         isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
 141         isl_int_set_si(path->ineq[k][1 + param], 1);
 142         isl_int_set_si(path->ineq[k][0], -1);
 143
 144         isl_dim_free(dim);
 145
 146         path = isl_basic_map_simplify(path);
 147         path = isl_basic_map_finalize(path);
 148         return isl_map_from_basic_map(path);
 149 error:
 150         isl_dim_free(dim);
 151         isl_basic_map_free(path);
 152         return NULL;
 153 }
 154
 155 /* Check whether any path of length at least one along "steps" is acyclic.
 156  * That is, check whether
 157  *
 158  *      \sum_i k_i \delta_i = 0
 159  *      \sum_i k_i >= 1
 160  *      k_i >= 0
 161  *
 162  * with \delta_i the rows of "steps", is infeasible.
 163  */
 164 static int is_acyclic(__isl_keep isl_mat *steps)
 165 {
 166         int acyclic;
 167         int i, j, k;
 168         struct isl_basic_set *test;
 169
 170         if (!steps)
 171                 return -1;
 172
 173         test = isl_basic_set_alloc(steps->ctx, 0, steps->n_row, 0,
 174                                         steps->n_col, steps->n_row + 1);
 175
 176         for (i = 0; i < steps->n_col; ++i) {
 177                 k = isl_basic_set_alloc_equality(test);
 178                 if (k < 0)
 179                         goto error;
 180                 isl_int_set_si(test->eq[k][0], 0);
 181                 for (j = 0; j < steps->n_row; ++j)
 182                         isl_int_set(test->eq[k][1 + j], steps->row[j][i]);
 183         }
 184         for (j = 0; j < steps->n_row; ++j) {
 185                 k = isl_basic_set_alloc_inequality(test);
 186                 if (k < 0)
 187                         goto error;
 188                 isl_seq_clr(test->ineq[k], 1 + steps->n_row);
 189                 isl_int_set_si(test->ineq[k][1 + j], 1);
 190         }
 191
 192         k = isl_basic_set_alloc_inequality(test);
 193         if (k < 0)
 194                 goto error;
 195         isl_int_set_si(test->ineq[k][0], -1);
 196         for (j = 0; j < steps->n_row; ++j)
 197                 isl_int_set_si(test->ineq[k][1 + j], 1);
 198
 199         acyclic = isl_basic_set_is_empty(test);
 200
 201         isl_basic_set_free(test);
 202         return acyclic;
 203 error:
 204         isl_basic_set_free(test);
 205         return -1;
 206 }
 207
 208 /* Shift variable at position "pos" up by one.
 209  * That is, replace the corresponding variable v by v - 1.
 210  */
 211 static __isl_give isl_basic_map *basic_map_shift_pos(
 212         __isl_take isl_basic_map *bmap, unsigned pos)
 213 {
 214         int i;
 215
 216         bmap = isl_basic_map_cow(bmap);
 217         if (!bmap)
 218                 return NULL;
 219
 220         for (i = 0; i < bmap->n_eq; ++i)
 221                 isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
 222
 223         for (i = 0; i < bmap->n_ineq; ++i)
 224                 isl_int_sub(bmap->ineq[i][0],
 225                                 bmap->ineq[i][0], bmap->ineq[i][pos]);
 226
 227         for (i = 0; i < bmap->n_div; ++i) {
 228                 if (isl_int_is_zero(bmap->div[i][0]))
 229                         continue;
 230                 isl_int_sub(bmap->div[i][1],
 231                                 bmap->div[i][1], bmap->div[i][1 + pos]);
 232         }
 233
 234         return bmap;
 235 }
 236
 237 /* Shift variable at position "pos" up by one.
 238  * That is, replace the corresponding variable v by v - 1.
 239  */
 240 static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
 241 {
 242         int i;
 243
 244         map = isl_map_cow(map);
 245         if (!map)
 246                 return NULL;
 247
 248         for (i = 0; i < map->n; ++i) {
 249                 map->p[i] = basic_map_shift_pos(map->p[i], pos);
 250                 if (!map->p[i])
 251                         goto error;
 252         }
 253         ISL_F_CLR(map, ISL_MAP_NORMALIZED);
 254         return map;
 255 error:
 256         isl_map_free(map);
 257         return NULL;
 258 }
 259
 260 /* Check whether the overapproximation of the power of "map" is exactly
 261  * the power of "map".  Let R be "map" and A_k the overapproximation.
 262  * The approximation is exact if
 263  *
 264  *      A_1 = R
 265  *      A_k = A_{k-1} \circ R                   k >= 2
 266  *
 267  * Since A_k is known to be an overapproximation, we only need to check
 268  *
 269  *      A_1 \subset R
 270  *      A_k \subset A_{k-1} \circ R             k >= 2
 271  *
 272  */
 273 static int check_power_exactness(__isl_take isl_map *map,
 274         __isl_take isl_map *app, unsigned param)
 275 {
 276         int exact;
 277         isl_map *app_1;
 278         isl_map *app_2;
 279
 280         app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
 281
 282         exact = isl_map_is_subset(app_1, map);
 283         isl_map_free(app_1);
 284
 285         if (!exact || exact < 0) {
 286                 isl_map_free(app);
 287                 isl_map_free(map);
 288                 return exact;
 289         }
 290
 291         app_2 = isl_map_lower_bound_si(isl_map_copy(app),
 292                                         isl_dim_param, param, 2);
 293         app_1 = map_shift_pos(app, 1 + param);
 294         app_1 = isl_map_apply_range(map, app_1);
 295
 296         exact = isl_map_is_subset(app_2, app_1);
 297
 298         isl_map_free(app_1);
 299         isl_map_free(app_2);
 300
 301         return exact;
 302 }
 303
 304 /* Check whether the overapproximation of the power of "map" is exactly
 305  * the power of "map", possibly after projecting out the power (if "project"
 306  * is set).
 307  *
 308  * If "project" is set and if "steps" can only result in acyclic paths,
 309  * then we check
 310  *
 311  *      A = R \cup (A \circ R)
 312  *
 313  * where A is the overapproximation with the power projected out, i.e.,
 314  * an overapproximation of the transitive closure.
 315  * More specifically, since A is known to be an overapproximation, we check
 316  *
 317  *      A \subset R \cup (A \circ R)
 318  *
 319  * Otherwise, we check if the power is exact.
 320  */
 321 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
 322         __isl_keep isl_mat *steps, unsigned param, int project)
 323 {
 324         isl_map *test;
 325         int exact;
 326
 327         if (project) {
 328                 project = is_acyclic(steps);
 329                 if (project < 0)
 330                         goto error;
 331         }
 332
 333         if (!project)
 334                 return check_power_exactness(map, app, param);
 335
 336         map = isl_map_project_out(map, isl_dim_param, param, 1);
 337         app = isl_map_project_out(app, isl_dim_param, param, 1);
 338
 339         test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
 340         test = isl_map_union(test, isl_map_copy(map));
 341
 342         exact = isl_map_is_subset(app, test);
 343
 344         isl_map_free(app);
 345         isl_map_free(test);
 346
 347         isl_map_free(map);
 348
 349         return exact;
 350 error:
 351         isl_map_free(app);
 352         isl_map_free(map);
 353         return -1;
 354 }
 355
 356 /* Compute the positive powers of "map", or an overapproximation.
 357  * The power is given by parameter "param".  If the result is exact,
 358  * then *exact is set to 1.
 359  * If project is set, then we are actually interested in the transitive
 360  * closure, so we can use a more relaxed exactness check.
 361  */
 362 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
 363         int *exact, int project)
 364 {
 365         struct isl_mat *steps = NULL;
 366         struct isl_set *domain = NULL;
 367         struct isl_set *range = NULL;
 368         struct isl_map *app = NULL;
 369         struct isl_map *path = NULL;
 370         int ok;
 371
 372         if (exact)
 373                 *exact = 1;
 374
 375         map = isl_map_remove_empty_parts(map);
 376         if (!map)
 377                 return NULL;
 378
 379         if (isl_map_fast_is_empty(map))
 380                 return map;
 381
 382         isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
 383         isl_assert(map->ctx,
 384                 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
 385                 goto error);
 386
 387         domain = isl_map_domain(isl_map_copy(map));
 388         domain = isl_set_coalesce(domain);
 389         range = isl_map_range(isl_map_copy(map));
 390         range = isl_set_coalesce(range);
 391         app = isl_map_from_domain_and_range(isl_set_copy(domain),
 392                                             isl_set_copy(range));
 393
 394         steps = extract_steps(map, &ok);
 395         if (!ok)
 396                 goto not_exact;
 397
 398         path = path_along_steps(isl_map_get_dim(map), steps, param);
 399         app = isl_map_intersect(app, isl_map_copy(path));
 400
 401         if (exact &&
 402             (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
 403                                   steps, param, project)) < 0)
 404                 goto error;
 405
 406         if (0) {
 407 not_exact:
 408                 if (exact)
 409                         *exact = 0;
 410         }
 411         isl_set_free(domain);
 412         isl_set_free(range);
 413         isl_mat_free(steps);
 414         isl_map_free(path);
 415         isl_map_free(map);
 416         return app;
 417 error:
 418         isl_set_free(domain);
 419         isl_set_free(range);
 420         isl_mat_free(steps);
 421         isl_map_free(path);
 422         isl_map_free(map);
 423         isl_map_free(app);
 424         return NULL;
 425 }
 426
 427 /* Compute the positive powers of "map", or an overapproximation.
 428  * The power is given by parameter "param".  If the result is exact,
 429  * then *exact is set to 1.
 430  */
 431 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
 432         int *exact)
 433 {
 434         return map_power(map, param, exact, 0);
 435 }
 436
 437 /* Compute the transitive closure  of "map", or an overapproximation.
 438  * If the result is exact, then *exact is set to 1.
 439  * Simply compute the powers of map and then project out the parameter
 440  * describing the power.
 441  */
 442 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
 443         int *exact)
 444 {
 445         unsigned param;
 446
 447         if (!map)
 448                 goto error;
 449
 450         param = isl_map_dim(map, isl_dim_param);
 451         map = isl_map_add(map, isl_dim_param, 1);
 452         map = map_power(map, param, exact, 1);
 453         map = isl_map_project_out(map, isl_dim_param, param, 1);
 454
 455         return map;
 456 error:
 457         isl_map_free(map);
 458         return NULL;
 459 }