1/* 2 * Copyright 2010 INRIA Saclay 3 * 4 * Use of this software is governed by the MIT 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_ctx_private.h> 12#include <isl_map_private.h> 13#include <isl/map.h> 14#include <isl/seq.h> 15#include <isl_space_private.h> 16#include <isl/lp.h> 17#include <isl/union_map.h> 18#include <isl_mat_private.h> 19#include <isl_options_private.h> 20#include <isl_tarjan.h> 21 22int isl_map_is_transitively_closed(__isl_keep isl_map *map) 23{ 24 isl_map *map2; 25 int closed; 26 27 map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map)); 28 closed = isl_map_is_subset(map2, map); 29 isl_map_free(map2); 30 31 return closed; 32} 33 34int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap) 35{ 36 isl_union_map *umap2; 37 int closed; 38 39 umap2 = isl_union_map_apply_range(isl_union_map_copy(umap), 40 isl_union_map_copy(umap)); 41 closed = isl_union_map_is_subset(umap2, umap); 42 isl_union_map_free(umap2); 43 44 return closed; 45} 46 47/* Given a map that represents a path with the length of the path 48 * encoded as the difference between the last output coordindate 49 * and the last input coordinate, set this length to either 50 * exactly "length" (if "exactly" is set) or at least "length" 51 * (if "exactly" is not set). 52 */ 53static __isl_give isl_map *set_path_length(__isl_take isl_map *map, 54 int exactly, int length) 55{ 56 isl_space *dim; 57 struct isl_basic_map *bmap; 58 unsigned d; 59 unsigned nparam; 60 int k; 61 isl_int *c; 62 63 if (!map) 64 return NULL; 65 66 dim = isl_map_get_space(map); 67 d = isl_space_dim(dim, isl_dim_in); 68 nparam = isl_space_dim(dim, isl_dim_param); 69 bmap = isl_basic_map_alloc_space(dim, 0, 1, 1); 70 if (exactly) { 71 k = isl_basic_map_alloc_equality(bmap); 72 c = bmap->eq[k]; 73 } else { 74 k = isl_basic_map_alloc_inequality(bmap); 75 c = bmap->ineq[k]; 76 } 77 if (k < 0) 78 goto error; 79 isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap)); 80 isl_int_set_si(c[0], -length); 81 isl_int_set_si(c[1 + nparam + d - 1], -1); 82 isl_int_set_si(c[1 + nparam + d + d - 1], 1); 83 84 bmap = isl_basic_map_finalize(bmap); 85 map = isl_map_intersect(map, isl_map_from_basic_map(bmap)); 86 87 return map; 88error: 89 isl_basic_map_free(bmap); 90 isl_map_free(map); 91 return NULL; 92} 93 94/* Check whether the overapproximation of the power of "map" is exactly 95 * the power of "map". Let R be "map" and A_k the overapproximation. 96 * The approximation is exact if 97 * 98 * A_1 = R 99 * A_k = A_{k-1} \circ R k >= 2 100 * 101 * Since A_k is known to be an overapproximation, we only need to check 102 * 103 * A_1 \subset R 104 * A_k \subset A_{k-1} \circ R k >= 2 105 * 106 * In practice, "app" has an extra input and output coordinate 107 * to encode the length of the path. So, we first need to add 108 * this coordinate to "map" and set the length of the path to 109 * one. 110 */ 111static int check_power_exactness(__isl_take isl_map *map, 112 __isl_take isl_map *app) 113{ 114 int exact; 115 isl_map *app_1; 116 isl_map *app_2; 117 118 map = isl_map_add_dims(map, isl_dim_in, 1); 119 map = isl_map_add_dims(map, isl_dim_out, 1); 120 map = set_path_length(map, 1, 1); 121 122 app_1 = set_path_length(isl_map_copy(app), 1, 1); 123 124 exact = isl_map_is_subset(app_1, map); 125 isl_map_free(app_1); 126 127 if (!exact || exact < 0) { 128 isl_map_free(app); 129 isl_map_free(map); 130 return exact; 131 } 132 133 app_1 = set_path_length(isl_map_copy(app), 0, 1); 134 app_2 = set_path_length(app, 0, 2); 135 app_1 = isl_map_apply_range(map, app_1); 136 137 exact = isl_map_is_subset(app_2, app_1); 138 139 isl_map_free(app_1); 140 isl_map_free(app_2); 141 142 return exact; 143} 144 145/* Check whether the overapproximation of the power of "map" is exactly 146 * the power of "map", possibly after projecting out the power (if "project" 147 * is set). 148 * 149 * If "project" is set and if "steps" can only result in acyclic paths, 150 * then we check 151 * 152 * A = R \cup (A \circ R) 153 * 154 * where A is the overapproximation with the power projected out, i.e., 155 * an overapproximation of the transitive closure. 156 * More specifically, since A is known to be an overapproximation, we check 157 * 158 * A \subset R \cup (A \circ R) 159 * 160 * Otherwise, we check if the power is exact. 161 * 162 * Note that "app" has an extra input and output coordinate to encode 163 * the length of the part. If we are only interested in the transitive 164 * closure, then we can simply project out these coordinates first. 165 */ 166static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app, 167 int project) 168{ 169 isl_map *test; 170 int exact; 171 unsigned d; 172 173 if (!project) 174 return check_power_exactness(map, app); 175 176 d = isl_map_dim(map, isl_dim_in); 177 app = set_path_length(app, 0, 1); 178 app = isl_map_project_out(app, isl_dim_in, d, 1); 179 app = isl_map_project_out(app, isl_dim_out, d, 1); 180 181 app = isl_map_reset_space(app, isl_map_get_space(map)); 182 183 test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app)); 184 test = isl_map_union(test, isl_map_copy(map)); 185 186 exact = isl_map_is_subset(app, test); 187 188 isl_map_free(app); 189 isl_map_free(test); 190 191 isl_map_free(map); 192 193 return exact; 194} 195 196/* 197 * The transitive closure implementation is based on the paper 198 * "Computing the Transitive Closure of a Union of Affine Integer 199 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and 200 * Albert Cohen. 201 */ 202 203/* Given a set of n offsets v_i (the rows of "steps"), construct a relation 204 * of the given dimension specification (Z^{n+1} -> Z^{n+1}) 205 * that maps an element x to any element that can be reached 206 * by taking a non-negative number of steps along any of 207 * the extended offsets v'_i = [v_i 1]. 208 * That is, construct 209 * 210 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i } 211 * 212 * For any element in this relation, the number of steps taken 213 * is equal to the difference in the final coordinates. 214 */ 215static __isl_give isl_map *path_along_steps(__isl_take isl_space *dim, 216 __isl_keep isl_mat *steps) 217{ 218 int i, j, k; 219 struct isl_basic_map *path = NULL; 220 unsigned d; 221 unsigned n; 222 unsigned nparam; 223 224 if (!dim || !steps) 225 goto error; 226 227 d = isl_space_dim(dim, isl_dim_in); 228 n = steps->n_row; 229 nparam = isl_space_dim(dim, isl_dim_param); 230 231 path = isl_basic_map_alloc_space(isl_space_copy(dim), n, d, n); 232 233 for (i = 0; i < n; ++i) { 234 k = isl_basic_map_alloc_div(path); 235 if (k < 0) 236 goto error; 237 isl_assert(steps->ctx, i == k, goto error); 238 isl_int_set_si(path->div[k][0], 0); 239 } 240 241 for (i = 0; i < d; ++i) { 242 k = isl_basic_map_alloc_equality(path); 243 if (k < 0) 244 goto error; 245 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path)); 246 isl_int_set_si(path->eq[k][1 + nparam + i], 1); 247 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1); 248 if (i == d - 1) 249 for (j = 0; j < n; ++j) 250 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1); 251 else 252 for (j = 0; j < n; ++j) 253 isl_int_set(path->eq[k][1 + nparam + 2 * d + j], 254 steps->row[j][i]); 255 } 256 257 for (i = 0; i < n; ++i) { 258 k = isl_basic_map_alloc_inequality(path); 259 if (k < 0) 260 goto error; 261 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path)); 262 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1); 263 } 264 265 isl_space_free(dim); 266 267 path = isl_basic_map_simplify(path); 268 path = isl_basic_map_finalize(path); 269 return isl_map_from_basic_map(path); 270error: 271 isl_space_free(dim); 272 isl_basic_map_free(path); 273 return NULL; 274} 275 276#define IMPURE 0 277#define PURE_PARAM 1 278#define PURE_VAR 2 279#define MIXED 3 280 281/* Check whether the parametric constant term of constraint c is never 282 * positive in "bset". 283 */ 284static int parametric_constant_never_positive(__isl_keep isl_basic_set *bset, 285 isl_int *c, int *div_purity) 286{ 287 unsigned d; 288 unsigned n_div; 289 unsigned nparam; 290 int i; 291 int k; 292 int empty; 293 294 n_div = isl_basic_set_dim(bset, isl_dim_div); 295 d = isl_basic_set_dim(bset, isl_dim_set); 296 nparam = isl_basic_set_dim(bset, isl_dim_param); 297 298 bset = isl_basic_set_copy(bset); 299 bset = isl_basic_set_cow(bset); 300 bset = isl_basic_set_extend_constraints(bset, 0, 1); 301 k = isl_basic_set_alloc_inequality(bset); 302 if (k < 0) 303 goto error; 304 isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset)); 305 isl_seq_cpy(bset->ineq[k], c, 1 + nparam); 306 for (i = 0; i < n_div; ++i) { 307 if (div_purity[i] != PURE_PARAM) 308 continue; 309 isl_int_set(bset->ineq[k][1 + nparam + d + i], 310 c[1 + nparam + d + i]); 311 } 312 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1); 313 empty = isl_basic_set_is_empty(bset); 314 isl_basic_set_free(bset); 315 316 return empty; 317error: 318 isl_basic_set_free(bset); 319 return -1; 320} 321 322/* Return PURE_PARAM if only the coefficients of the parameters are non-zero. 323 * Return PURE_VAR if only the coefficients of the set variables are non-zero. 324 * Return MIXED if only the coefficients of the parameters and the set 325 * variables are non-zero and if moreover the parametric constant 326 * can never attain positive values. 327 * Return IMPURE otherwise. 328 */ 329static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity, 330 int eq) 331{ 332 unsigned d; 333 unsigned n_div; 334 unsigned nparam; 335 int empty; 336 int i; 337 int p = 0, v = 0; 338 339 n_div = isl_basic_set_dim(bset, isl_dim_div); 340 d = isl_basic_set_dim(bset, isl_dim_set); 341 nparam = isl_basic_set_dim(bset, isl_dim_param); 342 343 for (i = 0; i < n_div; ++i) { 344 if (isl_int_is_zero(c[1 + nparam + d + i])) 345 continue; 346 switch (div_purity[i]) { 347 case PURE_PARAM: p = 1; break; 348 case PURE_VAR: v = 1; break; 349 default: return IMPURE; 350 } 351 } 352 if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1) 353 return PURE_VAR; 354 if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1) 355 return PURE_PARAM; 356 357 empty = parametric_constant_never_positive(bset, c, div_purity); 358 if (eq && empty >= 0 && !empty) { 359 isl_seq_neg(c, c, 1 + nparam + d + n_div); 360 empty = parametric_constant_never_positive(bset, c, div_purity); 361 } 362 363 return empty < 0 ? -1 : empty ? MIXED : IMPURE; 364} 365 366/* Return an array of integers indicating the type of each div in bset. 367 * If the div is (recursively) defined in terms of only the parameters, 368 * then the type is PURE_PARAM. 369 * If the div is (recursively) defined in terms of only the set variables, 370 * then the type is PURE_VAR. 371 * Otherwise, the type is IMPURE. 372 */ 373static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset) 374{ 375 int i, j; 376 int *div_purity; 377 unsigned d; 378 unsigned n_div; 379 unsigned nparam; 380 381 if (!bset) 382 return NULL; 383 384 n_div = isl_basic_set_dim(bset, isl_dim_div); 385 d = isl_basic_set_dim(bset, isl_dim_set); 386 nparam = isl_basic_set_dim(bset, isl_dim_param); 387 388 div_purity = isl_alloc_array(bset->ctx, int, n_div); 389 if (n_div && !div_purity) 390 return NULL; 391 392 for (i = 0; i < bset->n_div; ++i) { 393 int p = 0, v = 0; 394 if (isl_int_is_zero(bset->div[i][0])) { 395 div_purity[i] = IMPURE; 396 continue; 397 } 398 if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1) 399 p = 1; 400 if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1) 401 v = 1; 402 for (j = 0; j < i; ++j) { 403 if (isl_int_is_zero(bset->div[i][2 + nparam + d + j])) 404 continue; 405 switch (div_purity[j]) { 406 case PURE_PARAM: p = 1; break; 407 case PURE_VAR: v = 1; break; 408 default: p = v = 1; break; 409 } 410 } 411 div_purity[i] = v ? p ? IMPURE : PURE_VAR : PURE_PARAM; 412 } 413 414 return div_purity; 415} 416 417/* Given a path with the as yet unconstrained length at position "pos", 418 * check if setting the length to zero results in only the identity 419 * mapping. 420 */ 421static int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos) 422{ 423 isl_basic_map *test = NULL; 424 isl_basic_map *id = NULL; 425 int k; 426 int is_id; 427 428 test = isl_basic_map_copy(path); 429 test = isl_basic_map_extend_constraints(test, 1, 0); 430 k = isl_basic_map_alloc_equality(test); 431 if (k < 0) 432 goto error; 433 isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test)); 434 isl_int_set_si(test->eq[k][pos], 1); 435 id = isl_basic_map_identity(isl_basic_map_get_space(path)); 436 is_id = isl_basic_map_is_equal(test, id); 437 isl_basic_map_free(test); 438 isl_basic_map_free(id); 439 return is_id; 440error: 441 isl_basic_map_free(test); 442 return -1; 443} 444 445/* If any of the constraints is found to be impure then this function 446 * sets *impurity to 1. 447 * 448 * If impurity is NULL then we are dealing with a non-parametric set 449 * and so the constraints are obviously PURE_VAR. 450 */ 451static __isl_give isl_basic_map *add_delta_constraints( 452 __isl_take isl_basic_map *path, 453 __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam, 454 unsigned d, int *div_purity, int eq, int *impurity) 455{ 456 int i, k; 457 int n = eq ? delta->n_eq : delta->n_ineq; 458 isl_int **delta_c = eq ? delta->eq : delta->ineq; 459 unsigned n_div; 460 461 n_div = isl_basic_set_dim(delta, isl_dim_div); 462 463 for (i = 0; i < n; ++i) { 464 isl_int *path_c; 465 int p = PURE_VAR; 466 if (impurity) 467 p = purity(delta, delta_c[i], div_purity, eq); 468 if (p < 0) 469 goto error; 470 if (p != PURE_VAR && p != PURE_PARAM && !*impurity) 471 *impurity = 1; 472 if (p == IMPURE) 473 continue; 474 if (eq && p != MIXED) { 475 k = isl_basic_map_alloc_equality(path); 476 path_c = path->eq[k]; 477 } else { 478 k = isl_basic_map_alloc_inequality(path); 479 path_c = path->ineq[k]; 480 } 481 if (k < 0) 482 goto error; 483 isl_seq_clr(path_c, 1 + isl_basic_map_total_dim(path)); 484 if (p == PURE_VAR) { 485 isl_seq_cpy(path_c + off, 486 delta_c[i] + 1 + nparam, d); 487 isl_int_set(path_c[off + d], delta_c[i][0]); 488 } else if (p == PURE_PARAM) { 489 isl_seq_cpy(path_c, delta_c[i], 1 + nparam); 490 } else { 491 isl_seq_cpy(path_c + off, 492 delta_c[i] + 1 + nparam, d); 493 isl_seq_cpy(path_c, delta_c[i], 1 + nparam); 494 } 495 isl_seq_cpy(path_c + off - n_div, 496 delta_c[i] + 1 + nparam + d, n_div); 497 } 498 499 return path; 500error: 501 isl_basic_map_free(path); 502 return NULL; 503} 504 505/* Given a set of offsets "delta", construct a relation of the 506 * given dimension specification (Z^{n+1} -> Z^{n+1}) that 507 * is an overapproximation of the relations that 508 * maps an element x to any element that can be reached 509 * by taking a non-negative number of steps along any of 510 * the elements in "delta". 511 * That is, construct an approximation of 512 * 513 * { [x] -> [y] : exists f \in \delta, k \in Z : 514 * y = x + k [f, 1] and k >= 0 } 515 * 516 * For any element in this relation, the number of steps taken 517 * is equal to the difference in the final coordinates. 518 * 519 * In particular, let delta be defined as 520 * 521 * \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and 522 * C x + C'p + c >= 0 and 523 * D x + D'p + d >= 0 } 524 * 525 * where the constraints C x + C'p + c >= 0 are such that the parametric 526 * constant term of each constraint j, "C_j x + C'_j p + c_j", 527 * can never attain positive values, then the relation is constructed as 528 * 529 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and 530 * A f + k a >= 0 and B p + b >= 0 and 531 * C f + C'p + c >= 0 and k >= 1 } 532 * union { [x] -> [x] } 533 * 534 * If the zero-length paths happen to correspond exactly to the identity 535 * mapping, then we return 536 * 537 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and 538 * A f + k a >= 0 and B p + b >= 0 and 539 * C f + C'p + c >= 0 and k >= 0 } 540 * 541 * instead. 542 * 543 * Existentially quantified variables in \delta are handled by 544 * classifying them as independent of the parameters, purely 545 * parameter dependent and others. Constraints containing 546 * any of the other existentially quantified variables are removed. 547 * This is safe, but leads to an additional overapproximation. 548 * 549 * If there are any impure constraints, then we also eliminate 550 * the parameters from \delta, resulting in a set 551 * 552 * \delta' = { [x] : E x + e >= 0 } 553 * 554 * and add the constraints 555 * 556 * E f + k e >= 0 557 * 558 * to the constructed relation. 559 */ 560static __isl_give isl_map *path_along_delta(__isl_take isl_space *dim, 561 __isl_take isl_basic_set *delta) 562{ 563 isl_basic_map *path = NULL; 564 unsigned d; 565 unsigned n_div; 566 unsigned nparam; 567 unsigned off; 568 int i, k; 569 int is_id; 570 int *div_purity = NULL; 571 int impurity = 0; 572 573 if (!delta) 574 goto error; 575 n_div = isl_basic_set_dim(delta, isl_dim_div); 576 d = isl_basic_set_dim(delta, isl_dim_set); 577 nparam = isl_basic_set_dim(delta, isl_dim_param); 578 path = isl_basic_map_alloc_space(isl_space_copy(dim), n_div + d + 1, 579 d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1); 580 off = 1 + nparam + 2 * (d + 1) + n_div; 581 582 for (i = 0; i < n_div + d + 1; ++i) { 583 k = isl_basic_map_alloc_div(path); 584 if (k < 0) 585 goto error; 586 isl_int_set_si(path->div[k][0], 0); 587 } 588 589 for (i = 0; i < d + 1; ++i) { 590 k = isl_basic_map_alloc_equality(path); 591 if (k < 0) 592 goto error; 593 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path)); 594 isl_int_set_si(path->eq[k][1 + nparam + i], 1); 595 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1); 596 isl_int_set_si(path->eq[k][off + i], 1); 597 } 598 599 div_purity = get_div_purity(delta); 600 if (n_div && !div_purity) 601 goto error; 602 603 path = add_delta_constraints(path, delta, off, nparam, d, 604 div_purity, 1, &impurity); 605 path = add_delta_constraints(path, delta, off, nparam, d, 606 div_purity, 0, &impurity); 607 if (impurity) { 608 isl_space *dim = isl_basic_set_get_space(delta); 609 delta = isl_basic_set_project_out(delta, 610 isl_dim_param, 0, nparam); 611 delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam); 612 delta = isl_basic_set_reset_space(delta, dim); 613 if (!delta) 614 goto error; 615 path = isl_basic_map_extend_constraints(path, delta->n_eq, 616 delta->n_ineq + 1); 617 path = add_delta_constraints(path, delta, off, nparam, d, 618 NULL, 1, NULL); 619 path = add_delta_constraints(path, delta, off, nparam, d, 620 NULL, 0, NULL); 621 path = isl_basic_map_gauss(path, NULL); 622 } 623 624 is_id = empty_path_is_identity(path, off + d); 625 if (is_id < 0) 626 goto error; 627 628 k = isl_basic_map_alloc_inequality(path); 629 if (k < 0) 630 goto error; 631 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path)); 632 if (!is_id) 633 isl_int_set_si(path->ineq[k][0], -1); 634 isl_int_set_si(path->ineq[k][off + d], 1); 635 636 free(div_purity); 637 isl_basic_set_free(delta); 638 path = isl_basic_map_finalize(path); 639 if (is_id) { 640 isl_space_free(dim); 641 return isl_map_from_basic_map(path); 642 } 643 return isl_basic_map_union(path, isl_basic_map_identity(dim)); 644error: 645 free(div_purity); 646 isl_space_free(dim); 647 isl_basic_set_free(delta); 648 isl_basic_map_free(path); 649 return NULL; 650} 651 652/* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param", 653 * construct a map that equates the parameter to the difference 654 * in the final coordinates and imposes that this difference is positive. 655 * That is, construct 656 * 657 * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 } 658 */ 659static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_space *dim, 660 unsigned param) 661{ 662 struct isl_basic_map *bmap; 663 unsigned d; 664 unsigned nparam; 665 int k; 666 667 d = isl_space_dim(dim, isl_dim_in); 668 nparam = isl_space_dim(dim, isl_dim_param); 669 bmap = isl_basic_map_alloc_space(dim, 0, 1, 1); 670 k = isl_basic_map_alloc_equality(bmap); 671 if (k < 0) 672 goto error; 673 isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap)); 674 isl_int_set_si(bmap->eq[k][1 + param], -1); 675 isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1); 676 isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1); 677 678 k = isl_basic_map_alloc_inequality(bmap); 679 if (k < 0) 680 goto error; 681 isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap)); 682 isl_int_set_si(bmap->ineq[k][1 + param], 1); 683 isl_int_set_si(bmap->ineq[k][0], -1); 684 685 bmap = isl_basic_map_finalize(bmap); 686 return isl_map_from_basic_map(bmap); 687error: 688 isl_basic_map_free(bmap); 689 return NULL; 690} 691 692/* Check whether "path" is acyclic, where the last coordinates of domain 693 * and range of path encode the number of steps taken. 694 * That is, check whether 695 * 696 * { d | d = y - x and (x,y) in path } 697 * 698 * does not contain any element with positive last coordinate (positive length) 699 * and zero remaining coordinates (cycle). 700 */ 701static int is_acyclic(__isl_take isl_map *path) 702{ 703 int i; 704 int acyclic; 705 unsigned dim; 706 struct isl_set *delta; 707 708 delta = isl_map_deltas(path); 709 dim = isl_set_dim(delta, isl_dim_set); 710 for (i = 0; i < dim; ++i) { 711 if (i == dim -1) 712 delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1); 713 else 714 delta = isl_set_fix_si(delta, isl_dim_set, i, 0); 715 } 716 717 acyclic = isl_set_is_empty(delta); 718 isl_set_free(delta); 719 720 return acyclic; 721} 722 723/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D 724 * and a dimension specification (Z^{n+1} -> Z^{n+1}), 725 * construct a map that is an overapproximation of the map 726 * that takes an element from the space D \times Z to another 727 * element from the same space, such that the first n coordinates of the 728 * difference between them is a sum of differences between images 729 * and pre-images in one of the R_i and such that the last coordinate 730 * is equal to the number of steps taken. 731 * That is, let 732 * 733 * \Delta_i = { y - x | (x, y) in R_i } 734 * 735 * then the constructed map is an overapproximation of 736 * 737 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 738 * d = (\sum_i k_i \delta_i, \sum_i k_i) } 739 * 740 * The elements of the singleton \Delta_i's are collected as the 741 * rows of the steps matrix. For all these \Delta_i's together, 742 * a single path is constructed. 743 * For each of the other \Delta_i's, we compute an overapproximation 744 * of the paths along elements of \Delta_i. 745 * Since each of these paths performs an addition, composition is 746 * symmetric and we can simply compose all resulting paths in any order. 747 */ 748static __isl_give isl_map *construct_extended_path(__isl_take isl_space *dim, 749 __isl_keep isl_map *map, int *project) 750{ 751 struct isl_mat *steps = NULL; 752 struct isl_map *path = NULL; 753 unsigned d; 754 int i, j, n; 755 756 d = isl_map_dim(map, isl_dim_in); 757 758 path = isl_map_identity(isl_space_copy(dim)); 759 760 steps = isl_mat_alloc(map->ctx, map->n, d); 761 if (!steps) 762 goto error; 763 764 n = 0; 765 for (i = 0; i < map->n; ++i) { 766 struct isl_basic_set *delta; 767 768 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i])); 769 770 for (j = 0; j < d; ++j) { 771 int fixed; 772 773 fixed = isl_basic_set_plain_dim_is_fixed(delta, j, 774 &steps->row[n][j]); 775 if (fixed < 0) { 776 isl_basic_set_free(delta); 777 goto error; 778 } 779 if (!fixed) 780 break; 781 } 782 783 784 if (j < d) { 785 path = isl_map_apply_range(path, 786 path_along_delta(isl_space_copy(dim), delta)); 787 path = isl_map_coalesce(path); 788 } else { 789 isl_basic_set_free(delta); 790 ++n; 791 } 792 } 793 794 if (n > 0) { 795 steps->n_row = n; 796 path = isl_map_apply_range(path, 797 path_along_steps(isl_space_copy(dim), steps)); 798 } 799 800 if (project && *project) { 801 *project = is_acyclic(isl_map_copy(path)); 802 if (*project < 0) 803 goto error; 804 } 805 806 isl_space_free(dim); 807 isl_mat_free(steps); 808 return path; 809error: 810 isl_space_free(dim); 811 isl_mat_free(steps); 812 isl_map_free(path); 813 return NULL; 814} 815 816static int isl_set_overlaps(__isl_keep isl_set *set1, __isl_keep isl_set *set2) 817{ 818 isl_set *i; 819 int no_overlap; 820 821 if (!isl_space_tuple_match(set1->dim, isl_dim_set, set2->dim, isl_dim_set)) 822 return 0; 823 824 i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2)); 825 no_overlap = isl_set_is_empty(i); 826 isl_set_free(i); 827 828 return no_overlap < 0 ? -1 : !no_overlap; 829} 830 831/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D 832 * and a dimension specification (Z^{n+1} -> Z^{n+1}), 833 * construct a map that is an overapproximation of the map 834 * that takes an element from the dom R \times Z to an 835 * element from ran R \times Z, such that the first n coordinates of the 836 * difference between them is a sum of differences between images 837 * and pre-images in one of the R_i and such that the last coordinate 838 * is equal to the number of steps taken. 839 * That is, let 840 * 841 * \Delta_i = { y - x | (x, y) in R_i } 842 * 843 * then the constructed map is an overapproximation of 844 * 845 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 846 * d = (\sum_i k_i \delta_i, \sum_i k_i) and 847 * x in dom R and x + d in ran R and 848 * \sum_i k_i >= 1 } 849 */ 850static __isl_give isl_map *construct_component(__isl_take isl_space *dim, 851 __isl_keep isl_map *map, int *exact, int project) 852{ 853 struct isl_set *domain = NULL; 854 struct isl_set *range = NULL; 855 struct isl_map *app = NULL; 856 struct isl_map *path = NULL; 857 858 domain = isl_map_domain(isl_map_copy(map)); 859 domain = isl_set_coalesce(domain); 860 range = isl_map_range(isl_map_copy(map)); 861 range = isl_set_coalesce(range); 862 if (!isl_set_overlaps(domain, range)) { 863 isl_set_free(domain); 864 isl_set_free(range); 865 isl_space_free(dim); 866 867 map = isl_map_copy(map); 868 map = isl_map_add_dims(map, isl_dim_in, 1); 869 map = isl_map_add_dims(map, isl_dim_out, 1); 870 map = set_path_length(map, 1, 1); 871 return map; 872 } 873 app = isl_map_from_domain_and_range(domain, range); 874 app = isl_map_add_dims(app, isl_dim_in, 1); 875 app = isl_map_add_dims(app, isl_dim_out, 1); 876 877 path = construct_extended_path(isl_space_copy(dim), map, 878 exact && *exact ? &project : NULL); 879 app = isl_map_intersect(app, path); 880 881 if (exact && *exact && 882 (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app), 883 project)) < 0) 884 goto error; 885 886 isl_space_free(dim); 887 app = set_path_length(app, 0, 1); 888 return app; 889error: 890 isl_space_free(dim); 891 isl_map_free(app); 892 return NULL; 893} 894 895/* Call construct_component and, if "project" is set, project out 896 * the final coordinates. 897 */ 898static __isl_give isl_map *construct_projected_component( 899 __isl_take isl_space *dim, 900 __isl_keep isl_map *map, int *exact, int project) 901{ 902 isl_map *app; 903 unsigned d; 904 905 if (!dim) 906 return NULL; 907 d = isl_space_dim(dim, isl_dim_in); 908 909 app = construct_component(dim, map, exact, project); 910 if (project) { 911 app = isl_map_project_out(app, isl_dim_in, d - 1, 1); 912 app = isl_map_project_out(app, isl_dim_out, d - 1, 1); 913 } 914 return app; 915} 916 917/* Compute an extended version, i.e., with path lengths, of 918 * an overapproximation of the transitive closure of "bmap" 919 * with path lengths greater than or equal to zero and with 920 * domain and range equal to "dom". 921 */ 922static __isl_give isl_map *q_closure(__isl_take isl_space *dim, 923 __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact) 924{ 925 int project = 1; 926 isl_map *path; 927 isl_map *map; 928 isl_map *app; 929 930 dom = isl_set_add_dims(dom, isl_dim_set, 1); 931 app = isl_map_from_domain_and_range(dom, isl_set_copy(dom)); 932 map = isl_map_from_basic_map(isl_basic_map_copy(bmap)); 933 path = construct_extended_path(dim, map, &project); 934 app = isl_map_intersect(app, path); 935 936 if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0) 937 goto error; 938 939 return app; 940error: 941 isl_map_free(app); 942 return NULL; 943} 944 945/* Check whether qc has any elements of length at least one 946 * with domain and/or range outside of dom and ran. 947 */ 948static int has_spurious_elements(__isl_keep isl_map *qc, 949 __isl_keep isl_set *dom, __isl_keep isl_set *ran) 950{ 951 isl_set *s; 952 int subset; 953 unsigned d; 954 955 if (!qc || !dom || !ran) 956 return -1; 957 958 d = isl_map_dim(qc, isl_dim_in); 959 960 qc = isl_map_copy(qc); 961 qc = set_path_length(qc, 0, 1); 962 qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1); 963 qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1); 964 965 s = isl_map_domain(isl_map_copy(qc)); 966 subset = isl_set_is_subset(s, dom); 967 isl_set_free(s); 968 if (subset < 0) 969 goto error; 970 if (!subset) { 971 isl_map_free(qc); 972 return 1; 973 } 974 975 s = isl_map_range(qc); 976 subset = isl_set_is_subset(s, ran); 977 isl_set_free(s); 978 979 return subset < 0 ? -1 : !subset; 980error: 981 isl_map_free(qc); 982 return -1; 983} 984 985#define LEFT 2 986#define RIGHT 1 987 988/* For each basic map in "map", except i, check whether it combines 989 * with the transitive closure that is reflexive on C combines 990 * to the left and to the right. 991 * 992 * In particular, if 993 * 994 * dom map_j \subseteq C 995 * 996 * then right[j] is set to 1. Otherwise, if 997 * 998 * ran map_i \cap dom map_j = \emptyset 999 * 1000 * then right[j] is set to 0. Otherwise, composing to the right 1001 * is impossible. 1002 * 1003 * Similar, for composing to the left, we have if 1004 * 1005 * ran map_j \subseteq C 1006 * 1007 * then left[j] is set to 1. Otherwise, if 1008 * 1009 * dom map_i \cap ran map_j = \emptyset 1010 * 1011 * then left[j] is set to 0. Otherwise, composing to the left 1012 * is impossible. 1013 * 1014 * The return value is or'd with LEFT if composing to the left 1015 * is possible and with RIGHT if composing to the right is possible. 1016 */ 1017static int composability(__isl_keep isl_set *C, int i, 1018 isl_set **dom, isl_set **ran, int *left, int *right, 1019 __isl_keep isl_map *map) 1020{ 1021 int j; 1022 int ok; 1023 1024 ok = LEFT | RIGHT; 1025 for (j = 0; j < map->n && ok; ++j) { 1026 int overlaps, subset; 1027 if (j == i) 1028 continue; 1029 1030 if (ok & RIGHT) { 1031 if (!dom[j]) 1032 dom[j] = isl_set_from_basic_set( 1033 isl_basic_map_domain( 1034 isl_basic_map_copy(map->p[j]))); 1035 if (!dom[j]) 1036 return -1; 1037 overlaps = isl_set_overlaps(ran[i], dom[j]); 1038 if (overlaps < 0) 1039 return -1; 1040 if (!overlaps) 1041 right[j] = 0; 1042 else { 1043 subset = isl_set_is_subset(dom[j], C); 1044 if (subset < 0) 1045 return -1; 1046 if (subset) 1047 right[j] = 1; 1048 else 1049 ok &= ~RIGHT; 1050 } 1051 } 1052 1053 if (ok & LEFT) { 1054 if (!ran[j]) 1055 ran[j] = isl_set_from_basic_set( 1056 isl_basic_map_range( 1057 isl_basic_map_copy(map->p[j]))); 1058 if (!ran[j]) 1059 return -1; 1060 overlaps = isl_set_overlaps(dom[i], ran[j]); 1061 if (overlaps < 0) 1062 return -1; 1063 if (!overlaps) 1064 left[j] = 0; 1065 else { 1066 subset = isl_set_is_subset(ran[j], C); 1067 if (subset < 0) 1068 return -1; 1069 if (subset) 1070 left[j] = 1; 1071 else 1072 ok &= ~LEFT; 1073 } 1074 } 1075 } 1076 1077 return ok; 1078} 1079 1080static __isl_give isl_map *anonymize(__isl_take isl_map *map) 1081{ 1082 map = isl_map_reset(map, isl_dim_in); 1083 map = isl_map_reset(map, isl_dim_out); 1084 return map; 1085} 1086 1087/* Return a map that is a union of the basic maps in "map", except i, 1088 * composed to left and right with qc based on the entries of "left" 1089 * and "right". 1090 */ 1091static __isl_give isl_map *compose(__isl_keep isl_map *map, int i, 1092 __isl_take isl_map *qc, int *left, int *right) 1093{ 1094 int j; 1095 isl_map *comp; 1096 1097 comp = isl_map_empty(isl_map_get_space(map)); 1098 for (j = 0; j < map->n; ++j) { 1099 isl_map *map_j; 1100 1101 if (j == i) 1102 continue; 1103 1104 map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j])); 1105 map_j = anonymize(map_j); 1106 if (left && left[j]) 1107 map_j = isl_map_apply_range(map_j, isl_map_copy(qc)); 1108 if (right && right[j]) 1109 map_j = isl_map_apply_range(isl_map_copy(qc), map_j); 1110 comp = isl_map_union(comp, map_j); 1111 } 1112 1113 comp = isl_map_compute_divs(comp); 1114 comp = isl_map_coalesce(comp); 1115 1116 isl_map_free(qc); 1117 1118 return comp; 1119} 1120 1121/* Compute the transitive closure of "map" incrementally by 1122 * computing 1123 * 1124 * map_i^+ \cup qc^+ 1125 * 1126 * or 1127 * 1128 * map_i^+ \cup ((id \cup map_i^) \circ qc^+) 1129 * 1130 * or 1131 * 1132 * map_i^+ \cup (qc^+ \circ (id \cup map_i^)) 1133 * 1134 * depending on whether left or right are NULL. 1135 */ 1136static __isl_give isl_map *compute_incremental( 1137 __isl_take isl_space *dim, __isl_keep isl_map *map, 1138 int i, __isl_take isl_map *qc, int *left, int *right, int *exact) 1139{ 1140 isl_map *map_i; 1141 isl_map *tc; 1142 isl_map *rtc = NULL; 1143 1144 if (!map) 1145 goto error; 1146 isl_assert(map->ctx, left || right, goto error); 1147 1148 map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i])); 1149 tc = construct_projected_component(isl_space_copy(dim), map_i, 1150 exact, 1); 1151 isl_map_free(map_i); 1152 1153 if (*exact) 1154 qc = isl_map_transitive_closure(qc, exact); 1155 1156 if (!*exact) { 1157 isl_space_free(dim); 1158 isl_map_free(tc); 1159 isl_map_free(qc); 1160 return isl_map_universe(isl_map_get_space(map)); 1161 } 1162 1163 if (!left || !right) 1164 rtc = isl_map_union(isl_map_copy(tc), 1165 isl_map_identity(isl_map_get_space(tc))); 1166 if (!right) 1167 qc = isl_map_apply_range(rtc, qc); 1168 if (!left) 1169 qc = isl_map_apply_range(qc, rtc); 1170 qc = isl_map_union(tc, qc); 1171 1172 isl_space_free(dim); 1173 1174 return qc; 1175error: 1176 isl_space_free(dim); 1177 isl_map_free(qc); 1178 return NULL; 1179} 1180 1181/* Given a map "map", try to find a basic map such that 1182 * map^+ can be computed as 1183 * 1184 * map^+ = map_i^+ \cup 1185 * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+ 1186 * 1187 * with C the simple hull of the domain and range of the input map. 1188 * map_i^ \cup Id_C is computed by allowing the path lengths to be zero 1189 * and by intersecting domain and range with C. 1190 * Of course, we need to check that this is actually equal to map_i^ \cup Id_C. 1191 * Also, we only use the incremental computation if all the transitive 1192 * closures are exact and if the number of basic maps in the union, 1193 * after computing the integer divisions, is smaller than the number 1194 * of basic maps in the input map. 1195 */ 1196static int incemental_on_entire_domain(__isl_keep isl_space *dim, 1197 __isl_keep isl_map *map, 1198 isl_set **dom, isl_set **ran, int *left, int *right, 1199 __isl_give isl_map **res) 1200{ 1201 int i; 1202 isl_set *C; 1203 unsigned d; 1204 1205 *res = NULL; 1206 1207 C = isl_set_union(isl_map_domain(isl_map_copy(map)), 1208 isl_map_range(isl_map_copy(map))); 1209 C = isl_set_from_basic_set(isl_set_simple_hull(C)); 1210 if (!C) 1211 return -1; 1212 if (C->n != 1) { 1213 isl_set_free(C); 1214 return 0; 1215 } 1216 1217 d = isl_map_dim(map, isl_dim_in); 1218 1219 for (i = 0; i < map->n; ++i) { 1220 isl_map *qc; 1221 int exact_i, spurious; 1222 int j; 1223 dom[i] = isl_set_from_basic_set(isl_basic_map_domain( 1224 isl_basic_map_copy(map->p[i]))); 1225 ran[i] = isl_set_from_basic_set(isl_basic_map_range( 1226 isl_basic_map_copy(map->p[i]))); 1227 qc = q_closure(isl_space_copy(dim), isl_set_copy(C), 1228 map->p[i], &exact_i); 1229 if (!qc) 1230 goto error; 1231 if (!exact_i) { 1232 isl_map_free(qc); 1233 continue; 1234 } 1235 spurious = has_spurious_elements(qc, dom[i], ran[i]); 1236 if (spurious) { 1237 isl_map_free(qc); 1238 if (spurious < 0) 1239 goto error; 1240 continue; 1241 } 1242 qc = isl_map_project_out(qc, isl_dim_in, d, 1); 1243 qc = isl_map_project_out(qc, isl_dim_out, d, 1); 1244 qc = isl_map_compute_divs(qc); 1245 for (j = 0; j < map->n; ++j) 1246 left[j] = right[j] = 1; 1247 qc = compose(map, i, qc, left, right); 1248 if (!qc) 1249 goto error; 1250 if (qc->n >= map->n) { 1251 isl_map_free(qc); 1252 continue; 1253 } 1254 *res = compute_incremental(isl_space_copy(dim), map, i, qc, 1255 left, right, &exact_i); 1256 if (!*res) 1257 goto error; 1258 if (exact_i) 1259 break; 1260 isl_map_free(*res); 1261 *res = NULL; 1262 } 1263 1264 isl_set_free(C); 1265 1266 return *res != NULL; 1267error: 1268 isl_set_free(C); 1269 return -1; 1270} 1271 1272/* Try and compute the transitive closure of "map" as 1273 * 1274 * map^+ = map_i^+ \cup 1275 * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+ 1276 * 1277 * with C either the simple hull of the domain and range of the entire 1278 * map or the simple hull of domain and range of map_i. 1279 */ 1280static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim, 1281 __isl_keep isl_map *map, int *exact, int project) 1282{ 1283 int i; 1284 isl_set **dom = NULL; 1285 isl_set **ran = NULL; 1286 int *left = NULL; 1287 int *right = NULL; 1288 isl_set *C; 1289 unsigned d; 1290 isl_map *res = NULL; 1291 1292 if (!project) 1293 return construct_projected_component(dim, map, exact, project); 1294 1295 if (!map) 1296 goto error; 1297 if (map->n <= 1) 1298 return construct_projected_component(dim, map, exact, project); 1299 1300 d = isl_map_dim(map, isl_dim_in); 1301 1302 dom = isl_calloc_array(map->ctx, isl_set *, map->n); 1303 ran = isl_calloc_array(map->ctx, isl_set *, map->n); 1304 left = isl_calloc_array(map->ctx, int, map->n); 1305 right = isl_calloc_array(map->ctx, int, map->n); 1306 if (!ran || !dom || !left || !right) 1307 goto error; 1308 1309 if (incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 0) 1310 goto error; 1311 1312 for (i = 0; !res && i < map->n; ++i) { 1313 isl_map *qc; 1314 int exact_i, spurious, comp; 1315 if (!dom[i]) 1316 dom[i] = isl_set_from_basic_set( 1317 isl_basic_map_domain( 1318 isl_basic_map_copy(map->p[i]))); 1319 if (!dom[i]) 1320 goto error; 1321 if (!ran[i]) 1322 ran[i] = isl_set_from_basic_set( 1323 isl_basic_map_range( 1324 isl_basic_map_copy(map->p[i]))); 1325 if (!ran[i]) 1326 goto error; 1327 C = isl_set_union(isl_set_copy(dom[i]), 1328 isl_set_copy(ran[i])); 1329 C = isl_set_from_basic_set(isl_set_simple_hull(C)); 1330 if (!C) 1331 goto error; 1332 if (C->n != 1) { 1333 isl_set_free(C); 1334 continue; 1335 } 1336 comp = composability(C, i, dom, ran, left, right, map); 1337 if (!comp || comp < 0) { 1338 isl_set_free(C); 1339 if (comp < 0) 1340 goto error; 1341 continue; 1342 } 1343 qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i); 1344 if (!qc) 1345 goto error; 1346 if (!exact_i) { 1347 isl_map_free(qc); 1348 continue; 1349 } 1350 spurious = has_spurious_elements(qc, dom[i], ran[i]); 1351 if (spurious) { 1352 isl_map_free(qc); 1353 if (spurious < 0) 1354 goto error; 1355 continue; 1356 } 1357 qc = isl_map_project_out(qc, isl_dim_in, d, 1); 1358 qc = isl_map_project_out(qc, isl_dim_out, d, 1); 1359 qc = isl_map_compute_divs(qc); 1360 qc = compose(map, i, qc, (comp & LEFT) ? left : NULL, 1361 (comp & RIGHT) ? right : NULL); 1362 if (!qc) 1363 goto error; 1364 if (qc->n >= map->n) { 1365 isl_map_free(qc); 1366 continue; 1367 } 1368 res = compute_incremental(isl_space_copy(dim), map, i, qc, 1369 (comp & LEFT) ? left : NULL, 1370 (comp & RIGHT) ? right : NULL, &exact_i); 1371 if (!res) 1372 goto error; 1373 if (exact_i) 1374 break; 1375 isl_map_free(res); 1376 res = NULL; 1377 } 1378 1379 for (i = 0; i < map->n; ++i) { 1380 isl_set_free(dom[i]); 1381 isl_set_free(ran[i]); 1382 } 1383 free(dom); 1384 free(ran); 1385 free(left); 1386 free(right); 1387 1388 if (res) { 1389 isl_space_free(dim); 1390 return res; 1391 } 1392 1393 return construct_projected_component(dim, map, exact, project); 1394error: 1395 if (dom) 1396 for (i = 0; i < map->n; ++i) 1397 isl_set_free(dom[i]); 1398 free(dom); 1399 if (ran) 1400 for (i = 0; i < map->n; ++i) 1401 isl_set_free(ran[i]); 1402 free(ran); 1403 free(left); 1404 free(right); 1405 isl_space_free(dim); 1406 return NULL; 1407} 1408 1409/* Given an array of sets "set", add "dom" at position "pos" 1410 * and search for elements at earlier positions that overlap with "dom". 1411 * If any can be found, then merge all of them, together with "dom", into 1412 * a single set and assign the union to the first in the array, 1413 * which becomes the new group leader for all groups involved in the merge. 1414 * During the search, we only consider group leaders, i.e., those with 1415 * group[i] = i, as the other sets have already been combined 1416 * with one of the group leaders. 1417 */ 1418static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos) 1419{ 1420 int i; 1421 1422 group[pos] = pos; 1423 set[pos] = isl_set_copy(dom); 1424 1425 for (i = pos - 1; i >= 0; --i) { 1426 int o; 1427 1428 if (group[i] != i) 1429 continue; 1430 1431 o = isl_set_overlaps(set[i], dom); 1432 if (o < 0) 1433 goto error; 1434 if (!o) 1435 continue; 1436 1437 set[i] = isl_set_union(set[i], set[group[pos]]); 1438 set[group[pos]] = NULL; 1439 if (!set[i]) 1440 goto error; 1441 group[group[pos]] = i; 1442 group[pos] = i; 1443 } 1444 1445 isl_set_free(dom); 1446 return 0; 1447error: 1448 isl_set_free(dom); 1449 return -1; 1450} 1451 1452/* Replace each entry in the n by n grid of maps by the cross product 1453 * with the relation { [i] -> [i + 1] }. 1454 */ 1455static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n) 1456{ 1457 int i, j, k; 1458 isl_space *dim; 1459 isl_basic_map *bstep; 1460 isl_map *step; 1461 unsigned nparam; 1462 1463 if (!map) 1464 return -1; 1465 1466 dim = isl_map_get_space(map); 1467 nparam = isl_space_dim(dim, isl_dim_param); 1468 dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in)); 1469 dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out)); 1470 dim = isl_space_add_dims(dim, isl_dim_in, 1); 1471 dim = isl_space_add_dims(dim, isl_dim_out, 1); 1472 bstep = isl_basic_map_alloc_space(dim, 0, 1, 0); 1473 k = isl_basic_map_alloc_equality(bstep); 1474 if (k < 0) { 1475 isl_basic_map_free(bstep); 1476 return -1; 1477 } 1478 isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep)); 1479 isl_int_set_si(bstep->eq[k][0], 1); 1480 isl_int_set_si(bstep->eq[k][1 + nparam], 1); 1481 isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1); 1482 bstep = isl_basic_map_finalize(bstep); 1483 step = isl_map_from_basic_map(bstep); 1484 1485 for (i = 0; i < n; ++i) 1486 for (j = 0; j < n; ++j) 1487 grid[i][j] = isl_map_product(grid[i][j], 1488 isl_map_copy(step)); 1489 1490 isl_map_free(step); 1491 1492 return 0; 1493} 1494 1495/* The core of the Floyd-Warshall algorithm. 1496 * Updates the given n x x matrix of relations in place. 1497 * 1498 * The algorithm iterates over all vertices. In each step, the whole 1499 * matrix is updated to include all paths that go to the current vertex, 1500 * possibly stay there a while (including passing through earlier vertices) 1501 * and then come back. At the start of each iteration, the diagonal 1502 * element corresponding to the current vertex is replaced by its 1503 * transitive closure to account for all indirect paths that stay 1504 * in the current vertex. 1505 */ 1506static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact) 1507{ 1508 int r, p, q; 1509 1510 for (r = 0; r < n; ++r) { 1511 int r_exact; 1512 grid[r][r] = isl_map_transitive_closure(grid[r][r], 1513 (exact && *exact) ? &r_exact : NULL); 1514 if (exact && *exact && !r_exact) 1515 *exact = 0; 1516 1517 for (p = 0; p < n; ++p) 1518 for (q = 0; q < n; ++q) { 1519 isl_map *loop; 1520 if (p == r && q == r) 1521 continue; 1522 loop = isl_map_apply_range( 1523 isl_map_copy(grid[p][r]), 1524 isl_map_copy(grid[r][q])); 1525 grid[p][q] = isl_map_union(grid[p][q], loop); 1526 loop = isl_map_apply_range( 1527 isl_map_copy(grid[p][r]), 1528 isl_map_apply_range( 1529 isl_map_copy(grid[r][r]), 1530 isl_map_copy(grid[r][q]))); 1531 grid[p][q] = isl_map_union(grid[p][q], loop); 1532 grid[p][q] = isl_map_coalesce(grid[p][q]); 1533 } 1534 } 1535} 1536 1537/* Given a partition of the domains and ranges of the basic maps in "map", 1538 * apply the Floyd-Warshall algorithm with the elements in the partition 1539 * as vertices. 1540 * 1541 * In particular, there are "n" elements in the partition and "group" is 1542 * an array of length 2 * map->n with entries in [0,n-1]. 1543 * 1544 * We first construct a matrix of relations based on the partition information, 1545 * apply Floyd-Warshall on this matrix of relations and then take the 1546 * union of all entries in the matrix as the final result. 1547 * 1548 * If we are actually computing the power instead of the transitive closure, 1549 * i.e., when "project" is not set, then the result should have the 1550 * path lengths encoded as the difference between an extra pair of 1551 * coordinates. We therefore apply the nested transitive closures 1552 * to relations that include these lengths. In particular, we replace 1553 * the input relation by the cross product with the unit length relation 1554 * { [i] -> [i + 1] }. 1555 */ 1556static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_space *dim, 1557 __isl_keep isl_map *map, int *exact, int project, int *group, int n) 1558{ 1559 int i, j, k; 1560 isl_map ***grid = NULL; 1561 isl_map *app; 1562 1563 if (!map) 1564 goto error; 1565 1566 if (n == 1) { 1567 free(group); 1568 return incremental_closure(dim, map, exact, project); 1569 } 1570 1571 grid = isl_calloc_array(map->ctx, isl_map **, n); 1572 if (!grid) 1573 goto error; 1574 for (i = 0; i < n; ++i) { 1575 grid[i] = isl_calloc_array(map->ctx, isl_map *, n); 1576 if (!grid[i]) 1577 goto error; 1578 for (j = 0; j < n; ++j) 1579 grid[i][j] = isl_map_empty(isl_map_get_space(map)); 1580 } 1581 1582 for (k = 0; k < map->n; ++k) { 1583 i = group[2 * k]; 1584 j = group[2 * k + 1]; 1585 grid[i][j] = isl_map_union(grid[i][j], 1586 isl_map_from_basic_map( 1587 isl_basic_map_copy(map->p[k]))); 1588 } 1589 1590 if (!project && add_length(map, grid, n) < 0) 1591 goto error; 1592 1593 floyd_warshall_iterate(grid, n, exact); 1594 1595 app = isl_map_empty(isl_map_get_space(map)); 1596 1597 for (i = 0; i < n; ++i) { 1598 for (j = 0; j < n; ++j) 1599 app = isl_map_union(app, grid[i][j]); 1600 free(grid[i]); 1601 } 1602 free(grid); 1603 1604 free(group); 1605 isl_space_free(dim); 1606 1607 return app; 1608error: 1609 if (grid) 1610 for (i = 0; i < n; ++i) { 1611 if (!grid[i]) 1612 continue; 1613 for (j = 0; j < n; ++j) 1614 isl_map_free(grid[i][j]); 1615 free(grid[i]); 1616 } 1617 free(grid); 1618 free(group); 1619 isl_space_free(dim); 1620 return NULL; 1621} 1622 1623/* Partition the domains and ranges of the n basic relations in list 1624 * into disjoint cells. 1625 * 1626 * To find the partition, we simply consider all of the domains 1627 * and ranges in turn and combine those that overlap. 1628 * "set" contains the partition elements and "group" indicates 1629 * to which partition element a given domain or range belongs. 1630 * The domain of basic map i corresponds to element 2 * i in these arrays, 1631 * while the domain corresponds to element 2 * i + 1. 1632 * During the construction group[k] is either equal to k, 1633 * in which case set[k] contains the union of all the domains and 1634 * ranges in the corresponding group, or is equal to some l < k, 1635 * with l another domain or range in the same group. 1636 */ 1637static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n, 1638 isl_set ***set, int *n_group) 1639{ 1640 int i; 1641 int *group = NULL; 1642 int g; 1643 1644 *set = isl_calloc_array(ctx, isl_set *, 2 * n); 1645 group = isl_alloc_array(ctx, int, 2 * n); 1646 1647 if (!*set || !group) 1648 goto error; 1649 1650 for (i = 0; i < n; ++i) { 1651 isl_set *dom; 1652 dom = isl_set_from_basic_set(isl_basic_map_domain( 1653 isl_basic_map_copy(list[i]))); 1654 if (merge(*set, group, dom, 2 * i) < 0) 1655 goto error; 1656 dom = isl_set_from_basic_set(isl_basic_map_range( 1657 isl_basic_map_copy(list[i]))); 1658 if (merge(*set, group, dom, 2 * i + 1) < 0) 1659 goto error; 1660 } 1661 1662 g = 0; 1663 for (i = 0; i < 2 * n; ++i) 1664 if (group[i] == i) { 1665 if (g != i) { 1666 (*set)[g] = (*set)[i]; 1667 (*set)[i] = NULL; 1668 } 1669 group[i] = g++; 1670 } else 1671 group[i] = group[group[i]]; 1672 1673 *n_group = g; 1674 1675 return group; 1676error: 1677 if (*set) { 1678 for (i = 0; i < 2 * n; ++i) 1679 isl_set_free((*set)[i]); 1680 free(*set); 1681 *set = NULL; 1682 } 1683 free(group); 1684 return NULL; 1685} 1686 1687/* Check if the domains and ranges of the basic maps in "map" can 1688 * be partitioned, and if so, apply Floyd-Warshall on the elements 1689 * of the partition. Note that we also apply this algorithm 1690 * if we want to compute the power, i.e., when "project" is not set. 1691 * However, the results are unlikely to be exact since the recursive 1692 * calls inside the Floyd-Warshall algorithm typically result in 1693 * non-linear path lengths quite quickly. 1694 */ 1695static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim, 1696 __isl_keep isl_map *map, int *exact, int project) 1697{ 1698 int i; 1699 isl_set **set = NULL; 1700 int *group = NULL; 1701 int n; 1702 1703 if (!map) 1704 goto error; 1705 if (map->n <= 1) 1706 return incremental_closure(dim, map, exact, project); 1707 1708 group = setup_groups(map->ctx, map->p, map->n, &set, &n); 1709 if (!group) 1710 goto error; 1711 1712 for (i = 0; i < 2 * map->n; ++i) 1713 isl_set_free(set[i]); 1714 1715 free(set); 1716 1717 return floyd_warshall_with_groups(dim, map, exact, project, group, n); 1718error: 1719 isl_space_free(dim); 1720 return NULL; 1721} 1722 1723/* Structure for representing the nodes of the graph of which 1724 * strongly connected components are being computed. 1725 * 1726 * list contains the actual nodes 1727 * check_closed is set if we may have used the fact that 1728 * a pair of basic maps can be interchanged 1729 */ 1730struct isl_tc_follows_data { 1731 isl_basic_map **list; 1732 int check_closed; 1733}; 1734 1735/* Check whether in the computation of the transitive closure 1736 * "list[i]" (R_1) should follow (or be part of the same component as) 1737 * "list[j]" (R_2). 1738 * 1739 * That is check whether 1740 * 1741 * R_1 \circ R_2 1742 * 1743 * is a subset of 1744 * 1745 * R_2 \circ R_1 1746 * 1747 * If so, then there is no reason for R_1 to immediately follow R_2 1748 * in any path. 1749 * 1750 * *check_closed is set if the subset relation holds while 1751 * R_1 \circ R_2 is not empty. 1752 */ 1753static int basic_map_follows(int i, int j, void *user) 1754{ 1755 struct isl_tc_follows_data *data = user; 1756 struct isl_map *map12 = NULL; 1757 struct isl_map *map21 = NULL; 1758 int subset; 1759 1760 if (!isl_space_tuple_match(data->list[i]->dim, isl_dim_in, 1761 data->list[j]->dim, isl_dim_out)) 1762 return 0; 1763 1764 map21 = isl_map_from_basic_map( 1765 isl_basic_map_apply_range( 1766 isl_basic_map_copy(data->list[j]), 1767 isl_basic_map_copy(data->list[i]))); 1768 subset = isl_map_is_empty(map21); 1769 if (subset < 0) 1770 goto error; 1771 if (subset) { 1772 isl_map_free(map21); 1773 return 0; 1774 } 1775 1776 if (!isl_space_tuple_match(data->list[i]->dim, isl_dim_in, 1777 data->list[i]->dim, isl_dim_out) || 1778 !isl_space_tuple_match(data->list[j]->dim, isl_dim_in, 1779 data->list[j]->dim, isl_dim_out)) { 1780 isl_map_free(map21); 1781 return 1; 1782 } 1783 1784 map12 = isl_map_from_basic_map( 1785 isl_basic_map_apply_range( 1786 isl_basic_map_copy(data->list[i]), 1787 isl_basic_map_copy(data->list[j]))); 1788 1789 subset = isl_map_is_subset(map21, map12); 1790 1791 isl_map_free(map12); 1792 isl_map_free(map21); 1793 1794 if (subset) 1795 data->check_closed = 1; 1796 1797 return subset < 0 ? -1 : !subset; 1798error: 1799 isl_map_free(map21); 1800 return -1; 1801} 1802 1803/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D 1804 * and a dimension specification (Z^{n+1} -> Z^{n+1}), 1805 * construct a map that is an overapproximation of the map 1806 * that takes an element from the dom R \times Z to an 1807 * element from ran R \times Z, such that the first n coordinates of the 1808 * difference between them is a sum of differences between images 1809 * and pre-images in one of the R_i and such that the last coordinate 1810 * is equal to the number of steps taken. 1811 * If "project" is set, then these final coordinates are not included, 1812 * i.e., a relation of type Z^n -> Z^n is returned. 1813 * That is, let 1814 * 1815 * \Delta_i = { y - x | (x, y) in R_i } 1816 * 1817 * then the constructed map is an overapproximation of 1818 * 1819 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1820 * d = (\sum_i k_i \delta_i, \sum_i k_i) and 1821 * x in dom R and x + d in ran R } 1822 * 1823 * or 1824 * 1825 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1826 * d = (\sum_i k_i \delta_i) and 1827 * x in dom R and x + d in ran R } 1828 * 1829 * if "project" is set. 1830 * 1831 * We first split the map into strongly connected components, perform 1832 * the above on each component and then join the results in the correct 1833 * order, at each join also taking in the union of both arguments 1834 * to allow for paths that do not go through one of the two arguments. 1835 */ 1836static __isl_give isl_map *construct_power_components(__isl_take isl_space *dim, 1837 __isl_keep isl_map *map, int *exact, int project) 1838{ 1839 int i, n, c; 1840 struct isl_map *path = NULL; 1841 struct isl_tc_follows_data data; 1842 struct isl_tarjan_graph *g = NULL; 1843 int *orig_exact; 1844 int local_exact; 1845 1846 if (!map) 1847 goto error; 1848 if (map->n <= 1) 1849 return floyd_warshall(dim, map, exact, project); 1850 1851 data.list = map->p; 1852 data.check_closed = 0; 1853 g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data); 1854 if (!g) 1855 goto error; 1856 1857 orig_exact = exact; 1858 if (data.check_closed && !exact) 1859 exact = &local_exact; 1860 1861 c = 0; 1862 i = 0; 1863 n = map->n; 1864 if (project) 1865 path = isl_map_empty(isl_map_get_space(map)); 1866 else 1867 path = isl_map_empty(isl_space_copy(dim)); 1868 path = anonymize(path); 1869 while (n) { 1870 struct isl_map *comp; 1871 isl_map *path_comp, *path_comb; 1872 comp = isl_map_alloc_space(isl_map_get_space(map), n, 0); 1873 while (g->order[i] != -1) { 1874 comp = isl_map_add_basic_map(comp, 1875 isl_basic_map_copy(map->p[g->order[i]])); 1876 --n; 1877 ++i; 1878 } 1879 path_comp = floyd_warshall(isl_space_copy(dim), 1880 comp, exact, project); 1881 path_comp = anonymize(path_comp); 1882 path_comb = isl_map_apply_range(isl_map_copy(path), 1883 isl_map_copy(path_comp)); 1884 path = isl_map_union(path, path_comp); 1885 path = isl_map_union(path, path_comb); 1886 isl_map_free(comp); 1887 ++i; 1888 ++c; 1889 } 1890 1891 if (c > 1 && data.check_closed && !*exact) { 1892 int closed; 1893 1894 closed = isl_map_is_transitively_closed(path); 1895 if (closed < 0) 1896 goto error; 1897 if (!closed) { 1898 isl_tarjan_graph_free(g); 1899 isl_map_free(path); 1900 return floyd_warshall(dim, map, orig_exact, project); 1901 } 1902 } 1903 1904 isl_tarjan_graph_free(g); 1905 isl_space_free(dim); 1906 1907 return path; 1908error: 1909 isl_tarjan_graph_free(g); 1910 isl_space_free(dim); 1911 isl_map_free(path); 1912 return NULL; 1913} 1914 1915/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D, 1916 * construct a map that is an overapproximation of the map 1917 * that takes an element from the space D to another 1918 * element from the same space, such that the difference between 1919 * them is a strictly positive sum of differences between images 1920 * and pre-images in one of the R_i. 1921 * The number of differences in the sum is equated to parameter "param". 1922 * That is, let 1923 * 1924 * \Delta_i = { y - x | (x, y) in R_i } 1925 * 1926 * then the constructed map is an overapproximation of 1927 * 1928 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1929 * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 } 1930 * or 1931 * 1932 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1933 * d = \sum_i k_i \delta_i and \sum_i k_i > 0 } 1934 * 1935 * if "project" is set. 1936 * 1937 * If "project" is not set, then 1938 * we construct an extended mapping with an extra coordinate 1939 * that indicates the number of steps taken. In particular, 1940 * the difference in the last coordinate is equal to the number 1941 * of steps taken to move from a domain element to the corresponding 1942 * image element(s). 1943 */ 1944static __isl_give isl_map *construct_power(__isl_keep isl_map *map, 1945 int *exact, int project) 1946{ 1947 struct isl_map *app = NULL; 1948 isl_space *dim = NULL; 1949 unsigned d; 1950 1951 if (!map) 1952 return NULL; 1953 1954 dim = isl_map_get_space(map); 1955 1956 d = isl_space_dim(dim, isl_dim_in); 1957 dim = isl_space_add_dims(dim, isl_dim_in, 1); 1958 dim = isl_space_add_dims(dim, isl_dim_out, 1); 1959 1960 app = construct_power_components(isl_space_copy(dim), map, 1961 exact, project); 1962 1963 isl_space_free(dim); 1964 1965 return app; 1966} 1967 1968/* Compute the positive powers of "map", or an overapproximation. 1969 * If the result is exact, then *exact is set to 1. 1970 * 1971 * If project is set, then we are actually interested in the transitive 1972 * closure, so we can use a more relaxed exactness check. 1973 * The lengths of the paths are also projected out instead of being 1974 * encoded as the difference between an extra pair of final coordinates. 1975 */ 1976static __isl_give isl_map *map_power(__isl_take isl_map *map, 1977 int *exact, int project) 1978{ 1979 struct isl_map *app = NULL; 1980 1981 if (exact) 1982 *exact = 1; 1983 1984 if (!map) 1985 return NULL; 1986 1987 isl_assert(map->ctx, 1988 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out), 1989 goto error); 1990 1991 app = construct_power(map, exact, project); 1992 1993 isl_map_free(map); 1994 return app; 1995error: 1996 isl_map_free(map); 1997 isl_map_free(app); 1998 return NULL; 1999} 2000 2001/* Compute the positive powers of "map", or an overapproximation. 2002 * The result maps the exponent to a nested copy of the corresponding power. 2003 * If the result is exact, then *exact is set to 1. 2004 * map_power constructs an extended relation with the path lengths 2005 * encoded as the difference between the final coordinates. 2006 * In the final step, this difference is equated to an extra parameter 2007 * and made positive. The extra coordinates are subsequently projected out 2008 * and the parameter is turned into the domain of the result. 2009 */ 2010__isl_give isl_map *isl_map_power(__isl_take isl_map *map, int *exact) 2011{ 2012 isl_space *target_dim; 2013 isl_space *dim; 2014 isl_map *diff; 2015 unsigned d; 2016 unsigned param; 2017 2018 if (!map) 2019 return NULL; 2020 2021 d = isl_map_dim(map, isl_dim_in); 2022 param = isl_map_dim(map, isl_dim_param); 2023 2024 map = isl_map_compute_divs(map); 2025 map = isl_map_coalesce(map); 2026 2027 if (isl_map_plain_is_empty(map)) { 2028 map = isl_map_from_range(isl_map_wrap(map)); 2029 map = isl_map_add_dims(map, isl_dim_in, 1); 2030 map = isl_map_set_dim_name(map, isl_dim_in, 0, "k"); 2031 return map; 2032 } 2033 2034 target_dim = isl_map_get_space(map); 2035 target_dim = isl_space_from_range(isl_space_wrap(target_dim)); 2036 target_dim = isl_space_add_dims(target_dim, isl_dim_in, 1); 2037 target_dim = isl_space_set_dim_name(target_dim, isl_dim_in, 0, "k"); 2038 2039 map = map_power(map, exact, 0); 2040 2041 map = isl_map_add_dims(map, isl_dim_param, 1); 2042 dim = isl_map_get_space(map); 2043 diff = equate_parameter_to_length(dim, param); 2044 map = isl_map_intersect(map, diff); 2045 map = isl_map_project_out(map, isl_dim_in, d, 1); 2046 map = isl_map_project_out(map, isl_dim_out, d, 1); 2047 map = isl_map_from_range(isl_map_wrap(map)); 2048 map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1); 2049 2050 map = isl_map_reset_space(map, target_dim); 2051 2052 return map; 2053} 2054 2055/* Compute a relation that maps each element in the range of the input 2056 * relation to the lengths of all paths composed of edges in the input 2057 * relation that end up in the given range element. 2058 * The result may be an overapproximation, in which case *exact is set to 0. 2059 * The resulting relation is very similar to the power relation. 2060 * The difference are that the domain has been projected out, the 2061 * range has become the domain and the exponent is the range instead 2062 * of a parameter. 2063 */ 2064__isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map, 2065 int *exact) 2066{ 2067 isl_space *dim; 2068 isl_map *diff; 2069 unsigned d; 2070 unsigned param; 2071 2072 if (!map) 2073 return NULL; 2074 2075 d = isl_map_dim(map, isl_dim_in); 2076 param = isl_map_dim(map, isl_dim_param); 2077 2078 map = isl_map_compute_divs(map); 2079 map = isl_map_coalesce(map); 2080 2081 if (isl_map_plain_is_empty(map)) { 2082 if (exact) 2083 *exact = 1; 2084 map = isl_map_project_out(map, isl_dim_out, 0, d); 2085 map = isl_map_add_dims(map, isl_dim_out, 1); 2086 return map; 2087 } 2088 2089 map = map_power(map, exact, 0); 2090 2091 map = isl_map_add_dims(map, isl_dim_param, 1); 2092 dim = isl_map_get_space(map); 2093 diff = equate_parameter_to_length(dim, param); 2094 map = isl_map_intersect(map, diff); 2095 map = isl_map_project_out(map, isl_dim_in, 0, d + 1); 2096 map = isl_map_project_out(map, isl_dim_out, d, 1); 2097 map = isl_map_reverse(map); 2098 map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1); 2099 2100 return map; 2101} 2102 2103/* Check whether equality i of bset is a pure stride constraint 2104 * on a single dimensions, i.e., of the form 2105 * 2106 * v = k e 2107 * 2108 * with k a constant and e an existentially quantified variable. 2109 */ 2110static int is_eq_stride(__isl_keep isl_basic_set *bset, int i) 2111{ 2112 unsigned nparam; 2113 unsigned d; 2114 unsigned n_div; 2115 int pos1; 2116 int pos2; 2117 2118 if (!bset) 2119 return -1; 2120 2121 if (!isl_int_is_zero(bset->eq[i][0])) 2122 return 0; 2123 2124 nparam = isl_basic_set_dim(bset, isl_dim_param); 2125 d = isl_basic_set_dim(bset, isl_dim_set); 2126 n_div = isl_basic_set_dim(bset, isl_dim_div); 2127 2128 if (isl_seq_first_non_zero(bset->eq[i] + 1, nparam) != -1) 2129 return 0; 2130 pos1 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam, d); 2131 if (pos1 == -1) 2132 return 0; 2133 if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + pos1 + 1, 2134 d - pos1 - 1) != -1) 2135 return 0; 2136 2137 pos2 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d, n_div); 2138 if (pos2 == -1) 2139 return 0; 2140 if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d + pos2 + 1, 2141 n_div - pos2 - 1) != -1) 2142 return 0; 2143 if (!isl_int_is_one(bset->eq[i][1 + nparam + pos1]) && 2144 !isl_int_is_negone(bset->eq[i][1 + nparam + pos1])) 2145 return 0; 2146 2147 return 1; 2148} 2149 2150/* Given a map, compute the smallest superset of this map that is of the form 2151 * 2152 * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } 2153 * 2154 * (where p ranges over the (non-parametric) dimensions), 2155 * compute the transitive closure of this map, i.e., 2156 * 2157 * { i -> j : exists k > 0: 2158 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2159 * 2160 * and intersect domain and range of this transitive closure with 2161 * the given domain and range. 2162 * 2163 * If with_id is set, then try to include as much of the identity mapping 2164 * as possible, by computing 2165 * 2166 * { i -> j : exists k >= 0: 2167 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2168 * 2169 * instead (i.e., allow k = 0). 2170 * 2171 * In practice, we compute the difference set 2172 * 2173 * delta = { j - i | i -> j in map }, 2174 * 2175 * look for stride constraint on the individual dimensions and compute 2176 * (constant) lower and upper bounds for each individual dimension, 2177 * adding a constraint for each bound not equal to infinity. 2178 */ 2179static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map, 2180 __isl_take isl_set *dom, __isl_take isl_set *ran, int with_id) 2181{ 2182 int i; 2183 int k; 2184 unsigned d; 2185 unsigned nparam; 2186 unsigned total; 2187 isl_space *dim; 2188 isl_set *delta; 2189 isl_map *app = NULL; 2190 isl_basic_set *aff = NULL; 2191 isl_basic_map *bmap = NULL; 2192 isl_vec *obj = NULL; 2193 isl_int opt; 2194 2195 isl_int_init(opt); 2196 2197 delta = isl_map_deltas(isl_map_copy(map)); 2198 2199 aff = isl_set_affine_hull(isl_set_copy(delta)); 2200 if (!aff) 2201 goto error; 2202 dim = isl_map_get_space(map); 2203 d = isl_space_dim(dim, isl_dim_in); 2204 nparam = isl_space_dim(dim, isl_dim_param); 2205 total = isl_space_dim(dim, isl_dim_all); 2206 bmap = isl_basic_map_alloc_space(dim, 2207 aff->n_div + 1, aff->n_div, 2 * d + 1); 2208 for (i = 0; i < aff->n_div + 1; ++i) { 2209 k = isl_basic_map_alloc_div(bmap); 2210 if (k < 0) 2211 goto error; 2212 isl_int_set_si(bmap->div[k][0], 0); 2213 } 2214 for (i = 0; i < aff->n_eq; ++i) { 2215 if (!is_eq_stride(aff, i)) 2216 continue; 2217 k = isl_basic_map_alloc_equality(bmap); 2218 if (k < 0) 2219 goto error; 2220 isl_seq_clr(bmap->eq[k], 1 + nparam); 2221 isl_seq_cpy(bmap->eq[k] + 1 + nparam + d, 2222 aff->eq[i] + 1 + nparam, d); 2223 isl_seq_neg(bmap->eq[k] + 1 + nparam, 2224 aff->eq[i] + 1 + nparam, d); 2225 isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d, 2226 aff->eq[i] + 1 + nparam + d, aff->n_div); 2227 isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0); 2228 } 2229 obj = isl_vec_alloc(map->ctx, 1 + nparam + d); 2230 if (!obj) 2231 goto error; 2232 isl_seq_clr(obj->el, 1 + nparam + d); 2233 for (i = 0; i < d; ++ i) { 2234 enum isl_lp_result res; 2235 2236 isl_int_set_si(obj->el[1 + nparam + i], 1); 2237 2238 res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt, 2239 NULL, NULL); 2240 if (res == isl_lp_error) 2241 goto error; 2242 if (res == isl_lp_ok) { 2243 k = isl_basic_map_alloc_inequality(bmap); 2244 if (k < 0) 2245 goto error; 2246 isl_seq_clr(bmap->ineq[k], 2247 1 + nparam + 2 * d + bmap->n_div); 2248 isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1); 2249 isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1); 2250 isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt); 2251 } 2252 2253 res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt, 2254 NULL, NULL); 2255 if (res == isl_lp_error) 2256 goto error; 2257 if (res == isl_lp_ok) { 2258 k = isl_basic_map_alloc_inequality(bmap); 2259 if (k < 0) 2260 goto error; 2261 isl_seq_clr(bmap->ineq[k], 2262 1 + nparam + 2 * d + bmap->n_div); 2263 isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1); 2264 isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1); 2265 isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt); 2266 } 2267 2268 isl_int_set_si(obj->el[1 + nparam + i], 0); 2269 } 2270 k = isl_basic_map_alloc_inequality(bmap); 2271 if (k < 0) 2272 goto error; 2273 isl_seq_clr(bmap->ineq[k], 2274 1 + nparam + 2 * d + bmap->n_div); 2275 if (!with_id) 2276 isl_int_set_si(bmap->ineq[k][0], -1); 2277 isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1); 2278 2279 app = isl_map_from_domain_and_range(dom, ran); 2280 2281 isl_vec_free(obj); 2282 isl_basic_set_free(aff); 2283 isl_map_free(map); 2284 bmap = isl_basic_map_finalize(bmap); 2285 isl_set_free(delta); 2286 isl_int_clear(opt); 2287 2288 map = isl_map_from_basic_map(bmap); 2289 map = isl_map_intersect(map, app); 2290 2291 return map; 2292error: 2293 isl_vec_free(obj); 2294 isl_basic_map_free(bmap); 2295 isl_basic_set_free(aff); 2296 isl_set_free(dom); 2297 isl_set_free(ran); 2298 isl_map_free(map); 2299 isl_set_free(delta); 2300 isl_int_clear(opt); 2301 return NULL; 2302} 2303 2304/* Given a map, compute the smallest superset of this map that is of the form 2305 * 2306 * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } 2307 * 2308 * (where p ranges over the (non-parametric) dimensions), 2309 * compute the transitive closure of this map, i.e., 2310 * 2311 * { i -> j : exists k > 0: 2312 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2313 * 2314 * and intersect domain and range of this transitive closure with 2315 * domain and range of the original map. 2316 */ 2317static __isl_give isl_map *box_closure(__isl_take isl_map *map) 2318{ 2319 isl_set *domain; 2320 isl_set *range; 2321 2322 domain = isl_map_domain(isl_map_copy(map)); 2323 domain = isl_set_coalesce(domain); 2324 range = isl_map_range(isl_map_copy(map)); 2325 range = isl_set_coalesce(range); 2326 2327 return box_closure_on_domain(map, domain, range, 0); 2328} 2329 2330/* Given a map, compute the smallest superset of this map that is of the form 2331 * 2332 * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } 2333 * 2334 * (where p ranges over the (non-parametric) dimensions), 2335 * compute the transitive and partially reflexive closure of this map, i.e., 2336 * 2337 * { i -> j : exists k >= 0: 2338 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2339 * 2340 * and intersect domain and range of this transitive closure with 2341 * the given domain. 2342 */ 2343static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map, 2344 __isl_take isl_set *dom) 2345{ 2346 return box_closure_on_domain(map, dom, isl_set_copy(dom), 1); 2347} 2348 2349/* Check whether app is the transitive closure of map. 2350 * In particular, check that app is acyclic and, if so, 2351 * check that 2352 * 2353 * app \subset (map \cup (map \circ app)) 2354 */ 2355static int check_exactness_omega(__isl_keep isl_map *map, 2356 __isl_keep isl_map *app) 2357{ 2358 isl_set *delta; 2359 int i; 2360 int is_empty, is_exact; 2361 unsigned d; 2362 isl_map *test; 2363 2364 delta = isl_map_deltas(isl_map_copy(app)); 2365 d = isl_set_dim(delta, isl_dim_set); 2366 for (i = 0; i < d; ++i) 2367 delta = isl_set_fix_si(delta, isl_dim_set, i, 0); 2368 is_empty = isl_set_is_empty(delta); 2369 isl_set_free(delta); 2370 if (is_empty < 0) 2371 return -1; 2372 if (!is_empty) 2373 return 0; 2374 2375 test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map)); 2376 test = isl_map_union(test, isl_map_copy(map)); 2377 is_exact = isl_map_is_subset(app, test); 2378 isl_map_free(test); 2379 2380 return is_exact; 2381} 2382 2383/* Check if basic map M_i can be combined with all the other 2384 * basic maps such that 2385 * 2386 * (\cup_j M_j)^+ 2387 * 2388 * can be computed as 2389 * 2390 * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+ 2391 * 2392 * In particular, check if we can compute a compact representation 2393 * of 2394 * 2395 * M_i^* \circ M_j \circ M_i^* 2396 * 2397 * for each j != i. 2398 * Let M_i^? be an extension of M_i^+ that allows paths 2399 * of length zero, i.e., the result of box_closure(., 1). 2400 * The criterion, as proposed by Kelly et al., is that 2401 * id = M_i^? - M_i^+ can be represented as a basic map 2402 * and that 2403 * 2404 * id \circ M_j \circ id = M_j 2405 * 2406 * for each j != i. 2407 * 2408 * If this function returns 1, then tc and qc are set to 2409 * M_i^+ and M_i^?, respectively. 2410 */ 2411static int can_be_split_off(__isl_keep isl_map *map, int i, 2412 __isl_give isl_map **tc, __isl_give isl_map **qc) 2413{ 2414 isl_map *map_i, *id = NULL; 2415 int j = -1; 2416 isl_set *C; 2417 2418 *tc = NULL; 2419 *qc = NULL; 2420 2421 C = isl_set_union(isl_map_domain(isl_map_copy(map)), 2422 isl_map_range(isl_map_copy(map))); 2423 C = isl_set_from_basic_set(isl_set_simple_hull(C)); 2424 if (!C) 2425 goto error; 2426 2427 map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i])); 2428 *tc = box_closure(isl_map_copy(map_i)); 2429 *qc = box_closure_with_identity(map_i, C); 2430 id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc)); 2431 2432 if (!id || !*qc) 2433 goto error; 2434 if (id->n != 1 || (*qc)->n != 1) 2435 goto done; 2436 2437 for (j = 0; j < map->n; ++j) { 2438 isl_map *map_j, *test; 2439 int is_ok; 2440 2441 if (i == j) 2442 continue; 2443 map_j = isl_map_from_basic_map( 2444 isl_basic_map_copy(map->p[j])); 2445 test = isl_map_apply_range(isl_map_copy(id), 2446 isl_map_copy(map_j)); 2447 test = isl_map_apply_range(test, isl_map_copy(id)); 2448 is_ok = isl_map_is_equal(test, map_j); 2449 isl_map_free(map_j); 2450 isl_map_free(test); 2451 if (is_ok < 0) 2452 goto error; 2453 if (!is_ok) 2454 break; 2455 } 2456 2457done: 2458 isl_map_free(id); 2459 if (j == map->n) 2460 return 1; 2461 2462 isl_map_free(*qc); 2463 isl_map_free(*tc); 2464 *qc = NULL; 2465 *tc = NULL; 2466 2467 return 0; 2468error: 2469 isl_map_free(id); 2470 isl_map_free(*qc); 2471 isl_map_free(*tc); 2472 *qc = NULL; 2473 *tc = NULL; 2474 return -1; 2475} 2476 2477static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map, 2478 int *exact) 2479{ 2480 isl_map *app; 2481 2482 app = box_closure(isl_map_copy(map)); 2483 if (exact) 2484 *exact = check_exactness_omega(map, app); 2485 2486 isl_map_free(map); 2487 return app; 2488} 2489 2490/* Compute an overapproximation of the transitive closure of "map" 2491 * using a variation of the algorithm from 2492 * "Transitive Closure of Infinite Graphs and its Applications" 2493 * by Kelly et al. 2494 * 2495 * We first check whether we can can split of any basic map M_i and 2496 * compute 2497 * 2498 * (\cup_j M_j)^+ 2499 * 2500 * as 2501 * 2502 * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+ 2503 * 2504 * using a recursive call on the remaining map. 2505 * 2506 * If not, we simply call box_closure on the whole map. 2507 */ 2508static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map, 2509 int *exact) 2510{ 2511 int i, j; 2512 int exact_i; 2513 isl_map *app; 2514 2515 if (!map) 2516 return NULL; 2517 if (map->n == 1) 2518 return box_closure_with_check(map, exact); 2519 2520 for (i = 0; i < map->n; ++i) { 2521 int ok; 2522 isl_map *qc, *tc; 2523 ok = can_be_split_off(map, i, &tc, &qc); 2524 if (ok < 0) 2525 goto error; 2526 if (!ok) 2527 continue; 2528 2529 app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0); 2530 2531 for (j = 0; j < map->n; ++j) { 2532 if (j == i) 2533 continue; 2534 app = isl_map_add_basic_map(app, 2535 isl_basic_map_copy(map->p[j])); 2536 } 2537 2538 app = isl_map_apply_range(isl_map_copy(qc), app); 2539 app = isl_map_apply_range(app, qc); 2540 2541 app = isl_map_union(tc, transitive_closure_omega(app, NULL)); 2542 exact_i = check_exactness_omega(map, app); 2543 if (exact_i == 1) { 2544 if (exact) 2545 *exact = exact_i; 2546 isl_map_free(map); 2547 return app; 2548 } 2549 isl_map_free(app); 2550 if (exact_i < 0) 2551 goto error; 2552 } 2553 2554 return box_closure_with_check(map, exact); 2555error: 2556 isl_map_free(map); 2557 return NULL; 2558} 2559 2560/* Compute the transitive closure of "map", or an overapproximation. 2561 * If the result is exact, then *exact is set to 1. 2562 * Simply use map_power to compute the powers of map, but tell 2563 * it to project out the lengths of the paths instead of equating 2564 * the length to a parameter. 2565 */ 2566__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map, 2567 int *exact) 2568{ 2569 isl_space *target_dim; 2570 int closed; 2571 2572 if (!map) 2573 goto error; 2574 2575 if (map->ctx->opt->closure == ISL_CLOSURE_BOX) 2576 return transitive_closure_omega(map, exact); 2577 2578 map = isl_map_compute_divs(map); 2579 map = isl_map_coalesce(map); 2580 closed = isl_map_is_transitively_closed(map); 2581 if (closed < 0) 2582 goto error; 2583 if (closed) { 2584 if (exact) 2585 *exact = 1; 2586 return map; 2587 } 2588 2589 target_dim = isl_map_get_space(map); 2590 map = map_power(map, exact, 1); 2591 map = isl_map_reset_space(map, target_dim); 2592 2593 return map; 2594error: 2595 isl_map_free(map); 2596 return NULL; 2597} 2598 2599static int inc_count(__isl_take isl_map *map, void *user) 2600{ 2601 int *n = user; 2602 2603 *n += map->n; 2604 2605 isl_map_free(map); 2606 2607 return 0; 2608} 2609 2610static int collect_basic_map(__isl_take isl_map *map, void *user) 2611{ 2612 int i; 2613 isl_basic_map ***next = user; 2614 2615 for (i = 0; i < map->n; ++i) { 2616 **next = isl_basic_map_copy(map->p[i]); 2617 if (!**next) 2618 goto error; 2619 (*next)++; 2620 } 2621 2622 isl_map_free(map); 2623 return 0; 2624error: 2625 isl_map_free(map); 2626 return -1; 2627} 2628 2629/* Perform Floyd-Warshall on the given list of basic relations. 2630 * The basic relations may live in different dimensions, 2631 * but basic relations that get assigned to the diagonal of the 2632 * grid have domains and ranges of the same dimension and so 2633 * the standard algorithm can be used because the nested transitive 2634 * closures are only applied to diagonal elements and because all 2635 * compositions are peformed on relations with compatible domains and ranges. 2636 */ 2637static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx, 2638 __isl_keep isl_basic_map **list, int n, int *exact) 2639{ 2640 int i, j, k; 2641 int n_group; 2642 int *group = NULL; 2643 isl_set **set = NULL; 2644 isl_map ***grid = NULL; 2645 isl_union_map *app; 2646 2647 group = setup_groups(ctx, list, n, &set, &n_group); 2648 if (!group) 2649 goto error; 2650 2651 grid = isl_calloc_array(ctx, isl_map **, n_group); 2652 if (!grid) 2653 goto error; 2654 for (i = 0; i < n_group; ++i) { 2655 grid[i] = isl_calloc_array(ctx, isl_map *, n_group); 2656 if (!grid[i]) 2657 goto error; 2658 for (j = 0; j < n_group; ++j) { 2659 isl_space *dim1, *dim2, *dim; 2660 dim1 = isl_space_reverse(isl_set_get_space(set[i])); 2661 dim2 = isl_set_get_space(set[j]); 2662 dim = isl_space_join(dim1, dim2); 2663 grid[i][j] = isl_map_empty(dim); 2664 } 2665 } 2666 2667 for (k = 0; k < n; ++k) { 2668 i = group[2 * k]; 2669 j = group[2 * k + 1]; 2670 grid[i][j] = isl_map_union(grid[i][j], 2671 isl_map_from_basic_map( 2672 isl_basic_map_copy(list[k]))); 2673 } 2674 2675 floyd_warshall_iterate(grid, n_group, exact); 2676 2677 app = isl_union_map_empty(isl_map_get_space(grid[0][0])); 2678 2679 for (i = 0; i < n_group; ++i) { 2680 for (j = 0; j < n_group; ++j) 2681 app = isl_union_map_add_map(app, grid[i][j]); 2682 free(grid[i]); 2683 } 2684 free(grid); 2685 2686 for (i = 0; i < 2 * n; ++i) 2687 isl_set_free(set[i]); 2688 free(set); 2689 2690 free(group); 2691 return app; 2692error: 2693 if (grid) 2694 for (i = 0; i < n_group; ++i) { 2695 if (!grid[i]) 2696 continue; 2697 for (j = 0; j < n_group; ++j) 2698 isl_map_free(grid[i][j]); 2699 free(grid[i]); 2700 } 2701 free(grid); 2702 if (set) { 2703 for (i = 0; i < 2 * n; ++i) 2704 isl_set_free(set[i]); 2705 free(set); 2706 } 2707 free(group); 2708 return NULL; 2709} 2710 2711/* Perform Floyd-Warshall on the given union relation. 2712 * The implementation is very similar to that for non-unions. 2713 * The main difference is that it is applied unconditionally. 2714 * We first extract a list of basic maps from the union map 2715 * and then perform the algorithm on this list. 2716 */ 2717static __isl_give isl_union_map *union_floyd_warshall( 2718 __isl_take isl_union_map *umap, int *exact) 2719{ 2720 int i, n; 2721 isl_ctx *ctx; 2722 isl_basic_map **list = NULL; 2723 isl_basic_map **next; 2724 isl_union_map *res; 2725 2726 n = 0; 2727 if (isl_union_map_foreach_map(umap, inc_count, &n) < 0) 2728 goto error; 2729 2730 ctx = isl_union_map_get_ctx(umap); 2731 list = isl_calloc_array(ctx, isl_basic_map *, n); 2732 if (!list) 2733 goto error; 2734 2735 next = list; 2736 if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0) 2737 goto error; 2738 2739 res = union_floyd_warshall_on_list(ctx, list, n, exact); 2740 2741 if (list) { 2742 for (i = 0; i < n; ++i) 2743 isl_basic_map_free(list[i]); 2744 free(list); 2745 } 2746 2747 isl_union_map_free(umap); 2748 return res; 2749error: 2750 if (list) { 2751 for (i = 0; i < n; ++i) 2752 isl_basic_map_free(list[i]); 2753 free(list); 2754 } 2755 isl_union_map_free(umap); 2756 return NULL; 2757} 2758 2759/* Decompose the give union relation into strongly connected components. 2760 * The implementation is essentially the same as that of 2761 * construct_power_components with the major difference that all 2762 * operations are performed on union maps. 2763 */ 2764static __isl_give isl_union_map *union_components( 2765 __isl_take isl_union_map *umap, int *exact) 2766{ 2767 int i; 2768 int n; 2769 isl_ctx *ctx; 2770 isl_basic_map **list = NULL; 2771 isl_basic_map **next; 2772 isl_union_map *path = NULL; 2773 struct isl_tc_follows_data data; 2774 struct isl_tarjan_graph *g = NULL; 2775 int c, l; 2776 int recheck = 0; 2777 2778 n = 0; 2779 if (isl_union_map_foreach_map(umap, inc_count, &n) < 0) 2780 goto error; 2781 2782 if (n == 0) 2783 return umap; 2784 if (n <= 1) 2785 return union_floyd_warshall(umap, exact); 2786 2787 ctx = isl_union_map_get_ctx(umap); 2788 list = isl_calloc_array(ctx, isl_basic_map *, n); 2789 if (!list) 2790 goto error; 2791 2792 next = list; 2793 if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0) 2794 goto error; 2795 2796 data.list = list; 2797 data.check_closed = 0; 2798 g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data); 2799 if (!g) 2800 goto error; 2801 2802 c = 0; 2803 i = 0; 2804 l = n; 2805 path = isl_union_map_empty(isl_union_map_get_space(umap)); 2806 while (l) { 2807 isl_union_map *comp; 2808 isl_union_map *path_comp, *path_comb; 2809 comp = isl_union_map_empty(isl_union_map_get_space(umap)); 2810 while (g->order[i] != -1) { 2811 comp = isl_union_map_add_map(comp, 2812 isl_map_from_basic_map( 2813 isl_basic_map_copy(list[g->order[i]]))); 2814 --l; 2815 ++i; 2816 } 2817 path_comp = union_floyd_warshall(comp, exact); 2818 path_comb = isl_union_map_apply_range(isl_union_map_copy(path), 2819 isl_union_map_copy(path_comp)); 2820 path = isl_union_map_union(path, path_comp); 2821 path = isl_union_map_union(path, path_comb); 2822 ++i; 2823 ++c; 2824 } 2825 2826 if (c > 1 && data.check_closed && !*exact) { 2827 int closed; 2828 2829 closed = isl_union_map_is_transitively_closed(path); 2830 if (closed < 0) 2831 goto error; 2832 recheck = !closed; 2833 } 2834 2835 isl_tarjan_graph_free(g); 2836 2837 for (i = 0; i < n; ++i) 2838 isl_basic_map_free(list[i]); 2839 free(list); 2840 2841 if (recheck) { 2842 isl_union_map_free(path); 2843 return union_floyd_warshall(umap, exact); 2844 } 2845 2846 isl_union_map_free(umap); 2847 2848 return path; 2849error: 2850 isl_tarjan_graph_free(g); 2851 if (list) { 2852 for (i = 0; i < n; ++i) 2853 isl_basic_map_free(list[i]); 2854 free(list); 2855 } 2856 isl_union_map_free(umap); 2857 isl_union_map_free(path); 2858 return NULL; 2859} 2860 2861/* Compute the transitive closure of "umap", or an overapproximation. 2862 * If the result is exact, then *exact is set to 1. 2863 */ 2864__isl_give isl_union_map *isl_union_map_transitive_closure( 2865 __isl_take isl_union_map *umap, int *exact) 2866{ 2867 int closed; 2868 2869 if (!umap) 2870 return NULL; 2871 2872 if (exact) 2873 *exact = 1; 2874 2875 umap = isl_union_map_compute_divs(umap); 2876 umap = isl_union_map_coalesce(umap); 2877 closed = isl_union_map_is_transitively_closed(umap); 2878 if (closed < 0) 2879 goto error; 2880 if (closed) 2881 return umap; 2882 umap = union_components(umap, exact); 2883 return umap; 2884error: 2885 isl_union_map_free(umap); 2886 return NULL; 2887} 2888 2889struct isl_union_power { 2890 isl_union_map *pow; 2891 int *exact; 2892}; 2893 2894static int power(__isl_take isl_map *map, void *user) 2895{ 2896 struct isl_union_power *up = user; 2897 2898 map = isl_map_power(map, up->exact); 2899 up->pow = isl_union_map_from_map(map); 2900 2901 return -1; 2902} 2903 2904/* Construct a map [x] -> [x+1], with parameters prescribed by "dim". 2905 */ 2906static __isl_give isl_union_map *increment(__isl_take isl_space *dim) 2907{ 2908 int k; 2909 isl_basic_map *bmap; 2910 2911 dim = isl_space_add_dims(dim, isl_dim_in, 1); 2912 dim = isl_space_add_dims(dim, isl_dim_out, 1); 2913 bmap = isl_basic_map_alloc_space(dim, 0, 1, 0); 2914 k = isl_basic_map_alloc_equality(bmap); 2915 if (k < 0) 2916 goto error; 2917 isl_seq_clr(bmap->eq[k], isl_basic_map_total_dim(bmap)); 2918 isl_int_set_si(bmap->eq[k][0], 1); 2919 isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1); 2920 isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1); 2921 return isl_union_map_from_map(isl_map_from_basic_map(bmap)); 2922error: 2923 isl_basic_map_free(bmap); 2924 return NULL; 2925} 2926 2927/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim". 2928 */ 2929static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim) 2930{ 2931 isl_basic_map *bmap; 2932 2933 dim = isl_space_add_dims(dim, isl_dim_in, 1); 2934 dim = isl_space_add_dims(dim, isl_dim_out, 1); 2935 bmap = isl_basic_map_universe(dim); 2936 bmap = isl_basic_map_deltas_map(bmap); 2937 2938 return isl_union_map_from_map(isl_map_from_basic_map(bmap)); 2939} 2940 2941/* Compute the positive powers of "map", or an overapproximation. 2942 * The result maps the exponent to a nested copy of the corresponding power. 2943 * If the result is exact, then *exact is set to 1. 2944 */ 2945__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap, 2946 int *exact) 2947{ 2948 int n; 2949 isl_union_map *inc; 2950 isl_union_map *dm; 2951 2952 if (!umap) 2953 return NULL; 2954 n = isl_union_map_n_map(umap); 2955 if (n == 0) 2956 return umap; 2957 if (n == 1) { 2958 struct isl_union_power up = { NULL, exact }; 2959 isl_union_map_foreach_map(umap, &power, &up); 2960 isl_union_map_free(umap); 2961 return up.pow; 2962 } 2963 inc = increment(isl_union_map_get_space(umap)); 2964 umap = isl_union_map_product(inc, umap); 2965 umap = isl_union_map_transitive_closure(umap, exact); 2966 umap = isl_union_map_zip(umap); 2967 dm = deltas_map(isl_union_map_get_space(umap)); 2968 umap = isl_union_map_apply_domain(umap, dm); 2969 2970 return umap; 2971} 2972 2973#undef TYPE 2974#define TYPE isl_map 2975#include "isl_power_templ.c" 2976 2977#undef TYPE 2978#define TYPE isl_union_map 2979#include "isl_power_templ.c" 2980