154359Sroberto/* 254359Sroberto * Copyright 2010 INRIA Saclay 354359Sroberto * 454359Sroberto * Use of this software is governed by the MIT license 554359Sroberto * 654359Sroberto * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France, 754359Sroberto * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod, 854359Sroberto * 91893 Orsay, France 954359Sroberto */ 1054359Sroberto 1154359Sroberto#include <isl_ctx_private.h> 1254359Sroberto#include <isl_map_private.h> 1354359Sroberto#include <isl/map.h> 1454359Sroberto#include <isl_seq.h> 1554359Sroberto#include <isl_space_private.h> 1654359Sroberto#include <isl_lp_private.h> 1754359Sroberto#include <isl/union_map.h> 1854359Sroberto#include <isl_mat_private.h> 1954359Sroberto#include <isl_vec_private.h> 2054359Sroberto#include <isl_options_private.h> 2154359Sroberto#include <isl_tarjan.h> 2254359Sroberto 2354359Srobertoisl_bool isl_map_is_transitively_closed(__isl_keep isl_map *map) 2454359Sroberto{ 2554359Sroberto isl_map *map2; 2654359Sroberto isl_bool closed; 2754359Sroberto 2854359Sroberto map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map)); 2954359Sroberto closed = isl_map_is_subset(map2, map); 3054359Sroberto isl_map_free(map2); 3154359Sroberto 3254359Sroberto return closed; 3354359Sroberto} 3454359Sroberto 3554359Srobertoisl_bool isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap) 3654359Sroberto{ 3754359Sroberto isl_union_map *umap2; 3854359Sroberto isl_bool closed; 3954359Sroberto 4054359Sroberto umap2 = isl_union_map_apply_range(isl_union_map_copy(umap), 4154359Sroberto isl_union_map_copy(umap)); 4254359Sroberto closed = isl_union_map_is_subset(umap2, umap); 4354359Sroberto isl_union_map_free(umap2); 4454359Sroberto 4554359Sroberto return closed; 4654359Sroberto} 4754359Sroberto 4854359Sroberto/* Given a map that represents a path with the length of the path 4954359Sroberto * encoded as the difference between the last output coordindate 5054359Sroberto * and the last input coordinate, set this length to either 5154359Sroberto * exactly "length" (if "exactly" is set) or at least "length" 5254359Sroberto * (if "exactly" is not set). 5354359Sroberto */ 5454359Srobertostatic __isl_give isl_map *set_path_length(__isl_take isl_map *map, 5554359Sroberto int exactly, int length) 5654359Sroberto{ 5754359Sroberto isl_space *space; 5854359Sroberto struct isl_basic_map *bmap; 5954359Sroberto isl_size d; 6054359Sroberto isl_size nparam; 6154359Sroberto isl_size total; 6254359Sroberto int k; 6354359Sroberto isl_int *c; 6454359Sroberto 6554359Sroberto if (!map) 6654359Sroberto return NULL; 6754359Sroberto 6854359Sroberto space = isl_map_get_space(map); 6954359Sroberto d = isl_space_dim(space, isl_dim_in); 7054359Sroberto nparam = isl_space_dim(space, isl_dim_param); 7154359Sroberto total = isl_space_dim(space, isl_dim_all); 7254359Sroberto if (d < 0 || nparam < 0 || total < 0) 7354359Sroberto space = isl_space_free(space); 7454359Sroberto bmap = isl_basic_map_alloc_space(space, 0, 1, 1); 7554359Sroberto if (exactly) { 7654359Sroberto k = isl_basic_map_alloc_equality(bmap); 7754359Sroberto if (k < 0) 7854359Sroberto goto error; 7954359Sroberto c = bmap->eq[k]; 8054359Sroberto } else { 8154359Sroberto k = isl_basic_map_alloc_inequality(bmap); 8254359Sroberto if (k < 0) 8354359Sroberto goto error; 8454359Sroberto c = bmap->ineq[k]; 8554359Sroberto } 8654359Sroberto isl_seq_clr(c, 1 + total); 8754359Sroberto isl_int_set_si(c[0], -length); 8854359Sroberto isl_int_set_si(c[1 + nparam + d - 1], -1); 8954359Sroberto isl_int_set_si(c[1 + nparam + d + d - 1], 1); 9054359Sroberto 9154359Sroberto bmap = isl_basic_map_finalize(bmap); 9254359Sroberto map = isl_map_intersect(map, isl_map_from_basic_map(bmap)); 9354359Sroberto 9454359Sroberto return map; 9554359Srobertoerror: 9654359Sroberto isl_basic_map_free(bmap); 9754359Sroberto isl_map_free(map); 9854359Sroberto return NULL; 9954359Sroberto} 10054359Sroberto 10154359Sroberto/* Check whether the overapproximation of the power of "map" is exactly 10254359Sroberto * the power of "map". Let R be "map" and A_k the overapproximation. 10354359Sroberto * The approximation is exact if 10454359Sroberto * 10554359Sroberto * A_1 = R 10654359Sroberto * A_k = A_{k-1} \circ R k >= 2 10754359Sroberto * 10854359Sroberto * Since A_k is known to be an overapproximation, we only need to check 10954359Sroberto * 11054359Sroberto * A_1 \subset R 11154359Sroberto * A_k \subset A_{k-1} \circ R k >= 2 11254359Sroberto * 11354359Sroberto * In practice, "app" has an extra input and output coordinate 11454359Sroberto * to encode the length of the path. So, we first need to add 11554359Sroberto * this coordinate to "map" and set the length of the path to 11654359Sroberto * one. 11754359Sroberto */ 11854359Srobertostatic isl_bool check_power_exactness(__isl_take isl_map *map, 11954359Sroberto __isl_take isl_map *app) 12054359Sroberto{ 12154359Sroberto isl_bool exact; 12254359Sroberto isl_map *app_1; 12354359Sroberto isl_map *app_2; 12454359Sroberto 12554359Sroberto map = isl_map_add_dims(map, isl_dim_in, 1); 12654359Sroberto map = isl_map_add_dims(map, isl_dim_out, 1); 12754359Sroberto map = set_path_length(map, 1, 1); 12854359Sroberto 12954359Sroberto app_1 = set_path_length(isl_map_copy(app), 1, 1); 13054359Sroberto 13154359Sroberto exact = isl_map_is_subset(app_1, map); 13254359Sroberto isl_map_free(app_1); 13354359Sroberto 13454359Sroberto if (!exact || exact < 0) { 13554359Sroberto isl_map_free(app); 13654359Sroberto isl_map_free(map); 13754359Sroberto return exact; 13854359Sroberto } 13954359Sroberto 14054359Sroberto app_1 = set_path_length(isl_map_copy(app), 0, 1); 14154359Sroberto app_2 = set_path_length(app, 0, 2); 14254359Sroberto app_1 = isl_map_apply_range(map, app_1); 14354359Sroberto 14454359Sroberto exact = isl_map_is_subset(app_2, app_1); 14554359Sroberto 14654359Sroberto isl_map_free(app_1); 14754359Sroberto isl_map_free(app_2); 14854359Sroberto 14954359Sroberto return exact; 15054359Sroberto} 15154359Sroberto 15254359Sroberto/* Check whether the overapproximation of the power of "map" is exactly 15354359Sroberto * the power of "map", possibly after projecting out the power (if "project" 15454359Sroberto * is set). 15554359Sroberto * 15654359Sroberto * If "project" is set and if "steps" can only result in acyclic paths, 15754359Sroberto * then we check 15854359Sroberto * 15954359Sroberto * A = R \cup (A \circ R) 16054359Sroberto * 16154359Sroberto * where A is the overapproximation with the power projected out, i.e., 16254359Sroberto * an overapproximation of the transitive closure. 16354359Sroberto * More specifically, since A is known to be an overapproximation, we check 16454359Sroberto * 16554359Sroberto * A \subset R \cup (A \circ R) 16654359Sroberto * 16754359Sroberto * Otherwise, we check if the power is exact. 16854359Sroberto * 16954359Sroberto * Note that "app" has an extra input and output coordinate to encode 17054359Sroberto * the length of the part. If we are only interested in the transitive 17154359Sroberto * closure, then we can simply project out these coordinates first. 17254359Sroberto */ 17354359Srobertostatic isl_bool check_exactness(__isl_take isl_map *map, 17454359Sroberto __isl_take isl_map *app, int project) 17554359Sroberto{ 17654359Sroberto isl_map *test; 17754359Sroberto isl_bool exact; 17854359Sroberto isl_size d; 17954359Sroberto 18054359Sroberto if (!project) 18154359Sroberto return check_power_exactness(map, app); 18254359Sroberto 18354359Sroberto d = isl_map_dim(map, isl_dim_in); 18454359Sroberto if (d < 0) 18554359Sroberto app = isl_map_free(app); 18654359Sroberto app = set_path_length(app, 0, 1); 18754359Sroberto app = isl_map_project_out(app, isl_dim_in, d, 1); 18854359Sroberto app = isl_map_project_out(app, isl_dim_out, d, 1); 18954359Sroberto 19054359Sroberto app = isl_map_reset_space(app, isl_map_get_space(map)); 19154359Sroberto 19254359Sroberto test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app)); 19354359Sroberto test = isl_map_union(test, isl_map_copy(map)); 19454359Sroberto 19554359Sroberto exact = isl_map_is_subset(app, test); 19654359Sroberto 19754359Sroberto isl_map_free(app); 19854359Sroberto isl_map_free(test); 19954359Sroberto 20054359Sroberto isl_map_free(map); 20154359Sroberto 20254359Sroberto return exact; 20354359Sroberto} 20454359Sroberto 20554359Sroberto/* 20654359Sroberto * The transitive closure implementation is based on the paper 20754359Sroberto * "Computing the Transitive Closure of a Union of Affine Integer 20854359Sroberto * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and 20954359Sroberto * Albert Cohen. 21054359Sroberto */ 21154359Sroberto 21254359Sroberto/* Given a set of n offsets v_i (the rows of "steps"), construct a relation 21354359Sroberto * of the given dimension specification (Z^{n+1} -> Z^{n+1}) 21454359Sroberto * that maps an element x to any element that can be reached 21554359Sroberto * by taking a non-negative number of steps along any of 21654359Sroberto * the extended offsets v'_i = [v_i 1]. 21754359Sroberto * That is, construct 21854359Sroberto * 21954359Sroberto * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i } 22054359Sroberto * 22154359Sroberto * For any element in this relation, the number of steps taken 22254359Sroberto * is equal to the difference in the final coordinates. 22354359Sroberto */ 22454359Srobertostatic __isl_give isl_map *path_along_steps(__isl_take isl_space *space, 22554359Sroberto __isl_keep isl_mat *steps) 22654359Sroberto{ 22754359Sroberto int i, j, k; 22854359Sroberto struct isl_basic_map *path = NULL; 22954359Sroberto isl_size d; 23054359Sroberto unsigned n; 23154359Sroberto isl_size nparam; 23254359Sroberto isl_size total; 23354359Sroberto 23454359Sroberto d = isl_space_dim(space, isl_dim_in); 23554359Sroberto nparam = isl_space_dim(space, isl_dim_param); 23654359Sroberto if (d < 0 || nparam < 0 || !steps) 23754359Sroberto goto error; 23854359Sroberto 23954359Sroberto n = steps->n_row; 24054359Sroberto 24154359Sroberto path = isl_basic_map_alloc_space(isl_space_copy(space), n, d, n); 24254359Sroberto 24354359Sroberto for (i = 0; i < n; ++i) { 24454359Sroberto k = isl_basic_map_alloc_div(path); 24554359Sroberto if (k < 0) 24654359Sroberto goto error; 24754359Sroberto isl_assert(steps->ctx, i == k, goto error); 24854359Sroberto isl_int_set_si(path->div[k][0], 0); 24954359Sroberto } 25054359Sroberto 25154359Sroberto total = isl_basic_map_dim(path, isl_dim_all); 25254359Sroberto if (total < 0) 25354359Sroberto goto error; 25454359Sroberto for (i = 0; i < d; ++i) { 25554359Sroberto k = isl_basic_map_alloc_equality(path); 25654359Sroberto if (k < 0) 25754359Sroberto goto error; 25854359Sroberto isl_seq_clr(path->eq[k], 1 + total); 25954359Sroberto isl_int_set_si(path->eq[k][1 + nparam + i], 1); 26054359Sroberto isl_int_set_si(path->eq[k][1 + nparam + d + i], -1); 26154359Sroberto if (i == d - 1) 26254359Sroberto for (j = 0; j < n; ++j) 26354359Sroberto isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1); 26454359Sroberto else 26554359Sroberto for (j = 0; j < n; ++j) 26654359Sroberto isl_int_set(path->eq[k][1 + nparam + 2 * d + j], 26754359Sroberto steps->row[j][i]); 26854359Sroberto } 26954359Sroberto 27054359Sroberto for (i = 0; i < n; ++i) { 27154359Sroberto k = isl_basic_map_alloc_inequality(path); 27254359Sroberto if (k < 0) 27354359Sroberto goto error; 27454359Sroberto isl_seq_clr(path->ineq[k], 1 + total); 27554359Sroberto isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1); 27654359Sroberto } 27754359Sroberto 27854359Sroberto isl_space_free(space); 27954359Sroberto 28054359Sroberto path = isl_basic_map_simplify(path); 28154359Sroberto path = isl_basic_map_finalize(path); 28254359Sroberto return isl_map_from_basic_map(path); 28354359Srobertoerror: 28454359Sroberto isl_space_free(space); 28554359Sroberto isl_basic_map_free(path); 28654359Sroberto return NULL; 28754359Sroberto} 28854359Sroberto 28954359Sroberto#define IMPURE 0 29054359Sroberto#define PURE_PARAM 1 29154359Sroberto#define PURE_VAR 2 29254359Sroberto#define MIXED 3 29354359Sroberto 29454359Sroberto/* Check whether the parametric constant term of constraint c is never 29554359Sroberto * positive in "bset". 29654359Sroberto */ 29754359Srobertostatic isl_bool parametric_constant_never_positive( 29854359Sroberto __isl_keep isl_basic_set *bset, isl_int *c, int *div_purity) 29954359Sroberto{ 30054359Sroberto isl_size d; 30154359Sroberto isl_size n_div; 30254359Sroberto isl_size nparam; 30354359Sroberto isl_size total; 30454359Sroberto int i; 30554359Sroberto int k; 30654359Sroberto isl_bool empty; 30754359Sroberto 30854359Sroberto n_div = isl_basic_set_dim(bset, isl_dim_div); 30954359Sroberto d = isl_basic_set_dim(bset, isl_dim_set); 31054359Sroberto nparam = isl_basic_set_dim(bset, isl_dim_param); 31154359Sroberto total = isl_basic_set_dim(bset, isl_dim_all); 31254359Sroberto if (n_div < 0 || d < 0 || nparam < 0 || total < 0) 31354359Sroberto return isl_bool_error; 31454359Sroberto 31554359Sroberto bset = isl_basic_set_copy(bset); 31654359Sroberto bset = isl_basic_set_cow(bset); 31754359Sroberto bset = isl_basic_set_extend_constraints(bset, 0, 1); 31854359Sroberto k = isl_basic_set_alloc_inequality(bset); 31954359Sroberto if (k < 0) 32054359Sroberto goto error; 32154359Sroberto isl_seq_clr(bset->ineq[k], 1 + total); 32254359Sroberto isl_seq_cpy(bset->ineq[k], c, 1 + nparam); 32354359Sroberto for (i = 0; i < n_div; ++i) { 32454359Sroberto if (div_purity[i] != PURE_PARAM) 32554359Sroberto continue; 32654359Sroberto isl_int_set(bset->ineq[k][1 + nparam + d + i], 32754359Sroberto c[1 + nparam + d + i]); 32854359Sroberto } 32954359Sroberto isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1); 33054359Sroberto empty = isl_basic_set_is_empty(bset); 33154359Sroberto isl_basic_set_free(bset); 33254359Sroberto 33354359Sroberto return empty; 33454359Srobertoerror: 33554359Sroberto isl_basic_set_free(bset); 33654359Sroberto return isl_bool_error; 33754359Sroberto} 33854359Sroberto 33954359Sroberto/* Return PURE_PARAM if only the coefficients of the parameters are non-zero. 34054359Sroberto * Return PURE_VAR if only the coefficients of the set variables are non-zero. 34154359Sroberto * Return MIXED if only the coefficients of the parameters and the set 34254359Sroberto * variables are non-zero and if moreover the parametric constant 34354359Sroberto * can never attain positive values. 34454359Sroberto * Return IMPURE otherwise. 34554359Sroberto */ 34654359Srobertostatic int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity, 34754359Sroberto int eq) 34854359Sroberto{ 34954359Sroberto isl_size d; 35054359Sroberto isl_size n_div; 35154359Sroberto isl_size nparam; 35254359Sroberto isl_bool empty; 35354359Sroberto int i; 35454359Sroberto int p = 0, v = 0; 35554359Sroberto 35654359Sroberto n_div = isl_basic_set_dim(bset, isl_dim_div); 35754359Sroberto d = isl_basic_set_dim(bset, isl_dim_set); 35854359Sroberto nparam = isl_basic_set_dim(bset, isl_dim_param); 35954359Sroberto if (n_div < 0 || d < 0 || nparam < 0) 36054359Sroberto return -1; 36154359Sroberto 36254359Sroberto for (i = 0; i < n_div; ++i) { 36354359Sroberto if (isl_int_is_zero(c[1 + nparam + d + i])) 36454359Sroberto continue; 36554359Sroberto switch (div_purity[i]) { 36654359Sroberto case PURE_PARAM: p = 1; break; 36754359Sroberto case PURE_VAR: v = 1; break; 36854359Sroberto default: return IMPURE; 36954359Sroberto } 37054359Sroberto } 37154359Sroberto if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1) 37254359Sroberto return PURE_VAR; 37354359Sroberto if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1) 37454359Sroberto return PURE_PARAM; 37554359Sroberto 37654359Sroberto empty = parametric_constant_never_positive(bset, c, div_purity); 37754359Sroberto if (eq && empty >= 0 && !empty) { 37854359Sroberto isl_seq_neg(c, c, 1 + nparam + d + n_div); 37954359Sroberto empty = parametric_constant_never_positive(bset, c, div_purity); 38054359Sroberto } 38154359Sroberto 38254359Sroberto return empty < 0 ? -1 : empty ? MIXED : IMPURE; 38354359Sroberto} 38454359Sroberto 38554359Sroberto/* Return an array of integers indicating the type of each div in bset. 38654359Sroberto * If the div is (recursively) defined in terms of only the parameters, 38754359Sroberto * then the type is PURE_PARAM. 38854359Sroberto * If the div is (recursively) defined in terms of only the set variables, 38954359Sroberto * then the type is PURE_VAR. 39054359Sroberto * Otherwise, the type is IMPURE. 39154359Sroberto */ 39254359Srobertostatic __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset) 39354359Sroberto{ 39454359Sroberto int i, j; 39554359Sroberto int *div_purity; 39654359Sroberto isl_size d; 39754359Sroberto isl_size n_div; 39854359Sroberto isl_size nparam; 39954359Sroberto 40054359Sroberto n_div = isl_basic_set_dim(bset, isl_dim_div); 40154359Sroberto d = isl_basic_set_dim(bset, isl_dim_set); 40254359Sroberto nparam = isl_basic_set_dim(bset, isl_dim_param); 40354359Sroberto if (n_div < 0 || d < 0 || nparam < 0) 40454359Sroberto return NULL; 40554359Sroberto 40654359Sroberto div_purity = isl_alloc_array(bset->ctx, int, n_div); 40754359Sroberto if (n_div && !div_purity) 40854359Sroberto return NULL; 40954359Sroberto 41054359Sroberto for (i = 0; i < bset->n_div; ++i) { 41154359Sroberto int p = 0, v = 0; 41254359Sroberto if (isl_int_is_zero(bset->div[i][0])) { 41354359Sroberto div_purity[i] = IMPURE; 41454359Sroberto continue; 41554359Sroberto } 41654359Sroberto if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1) 41754359Sroberto p = 1; 41854359Sroberto if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1) 41954359Sroberto v = 1; 42054359Sroberto for (j = 0; j < i; ++j) { 42154359Sroberto if (isl_int_is_zero(bset->div[i][2 + nparam + d + j])) 42254359Sroberto continue; 42354359Sroberto switch (div_purity[j]) { 42454359Sroberto case PURE_PARAM: p = 1; break; 42554359Sroberto case PURE_VAR: v = 1; break; 42654359Sroberto default: p = v = 1; break; 42754359Sroberto } 42854359Sroberto } 42954359Sroberto div_purity[i] = v ? p ? IMPURE : PURE_VAR : PURE_PARAM; 43054359Sroberto } 43154359Sroberto 43254359Sroberto return div_purity; 43354359Sroberto} 43454359Sroberto 43554359Sroberto/* Given a path with the as yet unconstrained length at div position "pos", 43654359Sroberto * check if setting the length to zero results in only the identity 43754359Sroberto * mapping. 43854359Sroberto */ 43954359Srobertostatic isl_bool empty_path_is_identity(__isl_keep isl_basic_map *path, 44054359Sroberto unsigned pos) 44154359Sroberto{ 44254359Sroberto isl_basic_map *test = NULL; 44354359Sroberto isl_basic_map *id = NULL; 44454359Sroberto isl_bool is_id; 44554359Sroberto 44654359Sroberto test = isl_basic_map_copy(path); 44754359Sroberto test = isl_basic_map_fix_si(test, isl_dim_div, pos, 0); 44854359Sroberto id = isl_basic_map_identity(isl_basic_map_get_space(path)); 44954359Sroberto is_id = isl_basic_map_is_equal(test, id); 45054359Sroberto isl_basic_map_free(test); 45154359Sroberto isl_basic_map_free(id); 45254359Sroberto return is_id; 45354359Sroberto} 45454359Sroberto 45554359Sroberto/* If any of the constraints is found to be impure then this function 45654359Sroberto * sets *impurity to 1. 45754359Sroberto * 45854359Sroberto * If impurity is NULL then we are dealing with a non-parametric set 45954359Sroberto * and so the constraints are obviously PURE_VAR. 46054359Sroberto */ 46154359Srobertostatic __isl_give isl_basic_map *add_delta_constraints( 46254359Sroberto __isl_take isl_basic_map *path, 46354359Sroberto __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam, 46454359Sroberto unsigned d, int *div_purity, int eq, int *impurity) 46554359Sroberto{ 46654359Sroberto int i, k; 46754359Sroberto int n = eq ? delta->n_eq : delta->n_ineq; 46854359Sroberto isl_int **delta_c = eq ? delta->eq : delta->ineq; 46954359Sroberto isl_size n_div, total; 47054359Sroberto 47154359Sroberto n_div = isl_basic_set_dim(delta, isl_dim_div); 47254359Sroberto total = isl_basic_map_dim(path, isl_dim_all); 47354359Sroberto if (n_div < 0 || total < 0) 47454359Sroberto return isl_basic_map_free(path); 47554359Sroberto 47654359Sroberto for (i = 0; i < n; ++i) { 47754359Sroberto isl_int *path_c; 47854359Sroberto int p = PURE_VAR; 47954359Sroberto if (impurity) 48054359Sroberto p = purity(delta, delta_c[i], div_purity, eq); 48154359Sroberto if (p < 0) 48254359Sroberto goto error; 48354359Sroberto if (p != PURE_VAR && p != PURE_PARAM && !*impurity) 48454359Sroberto *impurity = 1; 48554359Sroberto if (p == IMPURE) 48654359Sroberto continue; 48754359Sroberto if (eq && p != MIXED) { 48854359Sroberto k = isl_basic_map_alloc_equality(path); 48954359Sroberto if (k < 0) 49054359Sroberto goto error; 49154359Sroberto path_c = path->eq[k]; 49254359Sroberto } else { 49354359Sroberto k = isl_basic_map_alloc_inequality(path); 49454359Sroberto if (k < 0) 49554359Sroberto goto error; 49654359Sroberto path_c = path->ineq[k]; 49754359Sroberto } 49854359Sroberto isl_seq_clr(path_c, 1 + total); 49954359Sroberto if (p == PURE_VAR) { 50054359Sroberto isl_seq_cpy(path_c + off, 50154359Sroberto delta_c[i] + 1 + nparam, d); 50254359Sroberto isl_int_set(path_c[off + d], delta_c[i][0]); 50354359Sroberto } else if (p == PURE_PARAM) { 50454359Sroberto isl_seq_cpy(path_c, delta_c[i], 1 + nparam); 50554359Sroberto } else { 50654359Sroberto isl_seq_cpy(path_c + off, 50754359Sroberto delta_c[i] + 1 + nparam, d); 50854359Sroberto isl_seq_cpy(path_c, delta_c[i], 1 + nparam); 50954359Sroberto } 51054359Sroberto isl_seq_cpy(path_c + off - n_div, 51154359Sroberto delta_c[i] + 1 + nparam + d, n_div); 51254359Sroberto } 51354359Sroberto 51454359Sroberto return path; 51554359Srobertoerror: 51654359Sroberto isl_basic_map_free(path); 51754359Sroberto return NULL; 51854359Sroberto} 51954359Sroberto 52054359Sroberto/* Given a set of offsets "delta", construct a relation of the 52154359Sroberto * given dimension specification (Z^{n+1} -> Z^{n+1}) that 52254359Sroberto * is an overapproximation of the relations that 52354359Sroberto * maps an element x to any element that can be reached 52454359Sroberto * by taking a non-negative number of steps along any of 52554359Sroberto * the elements in "delta". 52654359Sroberto * That is, construct an approximation of 52754359Sroberto * 52854359Sroberto * { [x] -> [y] : exists f \in \delta, k \in Z : 52954359Sroberto * y = x + k [f, 1] and k >= 0 } 53054359Sroberto * 53154359Sroberto * For any element in this relation, the number of steps taken 53254359Sroberto * is equal to the difference in the final coordinates. 53354359Sroberto * 53454359Sroberto * In particular, let delta be defined as 53554359Sroberto * 53654359Sroberto * \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and 53754359Sroberto * C x + C'p + c >= 0 and 53854359Sroberto * D x + D'p + d >= 0 } 53954359Sroberto * 54054359Sroberto * where the constraints C x + C'p + c >= 0 are such that the parametric 54154359Sroberto * constant term of each constraint j, "C_j x + C'_j p + c_j", 54254359Sroberto * can never attain positive values, then the relation is constructed as 54354359Sroberto * 54454359Sroberto * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and 545 * A f + k a >= 0 and B p + b >= 0 and 546 * C f + C'p + c >= 0 and k >= 1 } 547 * union { [x] -> [x] } 548 * 549 * If the zero-length paths happen to correspond exactly to the identity 550 * mapping, then we return 551 * 552 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and 553 * A f + k a >= 0 and B p + b >= 0 and 554 * C f + C'p + c >= 0 and k >= 0 } 555 * 556 * instead. 557 * 558 * Existentially quantified variables in \delta are handled by 559 * classifying them as independent of the parameters, purely 560 * parameter dependent and others. Constraints containing 561 * any of the other existentially quantified variables are removed. 562 * This is safe, but leads to an additional overapproximation. 563 * 564 * If there are any impure constraints, then we also eliminate 565 * the parameters from \delta, resulting in a set 566 * 567 * \delta' = { [x] : E x + e >= 0 } 568 * 569 * and add the constraints 570 * 571 * E f + k e >= 0 572 * 573 * to the constructed relation. 574 */ 575static __isl_give isl_map *path_along_delta(__isl_take isl_space *space, 576 __isl_take isl_basic_set *delta) 577{ 578 isl_basic_map *path = NULL; 579 isl_size d; 580 isl_size n_div; 581 isl_size nparam; 582 isl_size total; 583 unsigned off; 584 int i, k; 585 isl_bool is_id; 586 int *div_purity = NULL; 587 int impurity = 0; 588 589 n_div = isl_basic_set_dim(delta, isl_dim_div); 590 d = isl_basic_set_dim(delta, isl_dim_set); 591 nparam = isl_basic_set_dim(delta, isl_dim_param); 592 if (n_div < 0 || d < 0 || nparam < 0) 593 goto error; 594 path = isl_basic_map_alloc_space(isl_space_copy(space), n_div + d + 1, 595 d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1); 596 off = 1 + nparam + 2 * (d + 1) + n_div; 597 598 for (i = 0; i < n_div + d + 1; ++i) { 599 k = isl_basic_map_alloc_div(path); 600 if (k < 0) 601 goto error; 602 isl_int_set_si(path->div[k][0], 0); 603 } 604 605 total = isl_basic_map_dim(path, isl_dim_all); 606 if (total < 0) 607 goto error; 608 for (i = 0; i < d + 1; ++i) { 609 k = isl_basic_map_alloc_equality(path); 610 if (k < 0) 611 goto error; 612 isl_seq_clr(path->eq[k], 1 + total); 613 isl_int_set_si(path->eq[k][1 + nparam + i], 1); 614 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1); 615 isl_int_set_si(path->eq[k][off + i], 1); 616 } 617 618 div_purity = get_div_purity(delta); 619 if (n_div && !div_purity) 620 goto error; 621 622 path = add_delta_constraints(path, delta, off, nparam, d, 623 div_purity, 1, &impurity); 624 path = add_delta_constraints(path, delta, off, nparam, d, 625 div_purity, 0, &impurity); 626 if (impurity) { 627 isl_space *space = isl_basic_set_get_space(delta); 628 delta = isl_basic_set_project_out(delta, 629 isl_dim_param, 0, nparam); 630 delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam); 631 delta = isl_basic_set_reset_space(delta, space); 632 if (!delta) 633 goto error; 634 path = isl_basic_map_extend_constraints(path, delta->n_eq, 635 delta->n_ineq + 1); 636 path = add_delta_constraints(path, delta, off, nparam, d, 637 NULL, 1, NULL); 638 path = add_delta_constraints(path, delta, off, nparam, d, 639 NULL, 0, NULL); 640 path = isl_basic_map_gauss(path, NULL); 641 } 642 643 is_id = empty_path_is_identity(path, n_div + d); 644 if (is_id < 0) 645 goto error; 646 647 k = isl_basic_map_alloc_inequality(path); 648 if (k < 0) 649 goto error; 650 isl_seq_clr(path->ineq[k], 1 + total); 651 if (!is_id) 652 isl_int_set_si(path->ineq[k][0], -1); 653 isl_int_set_si(path->ineq[k][off + d], 1); 654 655 free(div_purity); 656 isl_basic_set_free(delta); 657 path = isl_basic_map_finalize(path); 658 if (is_id) { 659 isl_space_free(space); 660 return isl_map_from_basic_map(path); 661 } 662 return isl_basic_map_union(path, isl_basic_map_identity(space)); 663error: 664 free(div_purity); 665 isl_space_free(space); 666 isl_basic_set_free(delta); 667 isl_basic_map_free(path); 668 return NULL; 669} 670 671/* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param", 672 * construct a map that equates the parameter to the difference 673 * in the final coordinates and imposes that this difference is positive. 674 * That is, construct 675 * 676 * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 } 677 */ 678static __isl_give isl_map *equate_parameter_to_length( 679 __isl_take isl_space *space, unsigned param) 680{ 681 struct isl_basic_map *bmap; 682 isl_size d; 683 isl_size nparam; 684 isl_size total; 685 int k; 686 687 d = isl_space_dim(space, isl_dim_in); 688 nparam = isl_space_dim(space, isl_dim_param); 689 total = isl_space_dim(space, isl_dim_all); 690 if (d < 0 || nparam < 0 || total < 0) 691 space = isl_space_free(space); 692 bmap = isl_basic_map_alloc_space(space, 0, 1, 1); 693 k = isl_basic_map_alloc_equality(bmap); 694 if (k < 0) 695 goto error; 696 isl_seq_clr(bmap->eq[k], 1 + total); 697 isl_int_set_si(bmap->eq[k][1 + param], -1); 698 isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1); 699 isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1); 700 701 k = isl_basic_map_alloc_inequality(bmap); 702 if (k < 0) 703 goto error; 704 isl_seq_clr(bmap->ineq[k], 1 + total); 705 isl_int_set_si(bmap->ineq[k][1 + param], 1); 706 isl_int_set_si(bmap->ineq[k][0], -1); 707 708 bmap = isl_basic_map_finalize(bmap); 709 return isl_map_from_basic_map(bmap); 710error: 711 isl_basic_map_free(bmap); 712 return NULL; 713} 714 715/* Check whether "path" is acyclic, where the last coordinates of domain 716 * and range of path encode the number of steps taken. 717 * That is, check whether 718 * 719 * { d | d = y - x and (x,y) in path } 720 * 721 * does not contain any element with positive last coordinate (positive length) 722 * and zero remaining coordinates (cycle). 723 */ 724static isl_bool is_acyclic(__isl_take isl_map *path) 725{ 726 int i; 727 isl_bool acyclic; 728 isl_size dim; 729 struct isl_set *delta; 730 731 delta = isl_map_deltas(path); 732 dim = isl_set_dim(delta, isl_dim_set); 733 if (dim < 0) 734 delta = isl_set_free(delta); 735 for (i = 0; i < dim; ++i) { 736 if (i == dim -1) 737 delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1); 738 else 739 delta = isl_set_fix_si(delta, isl_dim_set, i, 0); 740 } 741 742 acyclic = isl_set_is_empty(delta); 743 isl_set_free(delta); 744 745 return acyclic; 746} 747 748/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D 749 * and a dimension specification (Z^{n+1} -> Z^{n+1}), 750 * construct a map that is an overapproximation of the map 751 * that takes an element from the space D \times Z to another 752 * element from the same space, such that the first n coordinates of the 753 * difference between them is a sum of differences between images 754 * and pre-images in one of the R_i and such that the last coordinate 755 * is equal to the number of steps taken. 756 * That is, let 757 * 758 * \Delta_i = { y - x | (x, y) in R_i } 759 * 760 * then the constructed map is an overapproximation of 761 * 762 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 763 * d = (\sum_i k_i \delta_i, \sum_i k_i) } 764 * 765 * The elements of the singleton \Delta_i's are collected as the 766 * rows of the steps matrix. For all these \Delta_i's together, 767 * a single path is constructed. 768 * For each of the other \Delta_i's, we compute an overapproximation 769 * of the paths along elements of \Delta_i. 770 * Since each of these paths performs an addition, composition is 771 * symmetric and we can simply compose all resulting paths in any order. 772 */ 773static __isl_give isl_map *construct_extended_path(__isl_take isl_space *space, 774 __isl_keep isl_map *map, int *project) 775{ 776 struct isl_mat *steps = NULL; 777 struct isl_map *path = NULL; 778 isl_size d; 779 int i, j, n; 780 781 d = isl_map_dim(map, isl_dim_in); 782 if (d < 0) 783 goto error; 784 785 path = isl_map_identity(isl_space_copy(space)); 786 787 steps = isl_mat_alloc(map->ctx, map->n, d); 788 if (!steps) 789 goto error; 790 791 n = 0; 792 for (i = 0; i < map->n; ++i) { 793 struct isl_basic_set *delta; 794 795 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i])); 796 797 for (j = 0; j < d; ++j) { 798 isl_bool fixed; 799 800 fixed = isl_basic_set_plain_dim_is_fixed(delta, j, 801 &steps->row[n][j]); 802 if (fixed < 0) { 803 isl_basic_set_free(delta); 804 goto error; 805 } 806 if (!fixed) 807 break; 808 } 809 810 811 if (j < d) { 812 path = isl_map_apply_range(path, 813 path_along_delta(isl_space_copy(space), delta)); 814 path = isl_map_coalesce(path); 815 } else { 816 isl_basic_set_free(delta); 817 ++n; 818 } 819 } 820 821 if (n > 0) { 822 steps->n_row = n; 823 path = isl_map_apply_range(path, 824 path_along_steps(isl_space_copy(space), steps)); 825 } 826 827 if (project && *project) { 828 *project = is_acyclic(isl_map_copy(path)); 829 if (*project < 0) 830 goto error; 831 } 832 833 isl_space_free(space); 834 isl_mat_free(steps); 835 return path; 836error: 837 isl_space_free(space); 838 isl_mat_free(steps); 839 isl_map_free(path); 840 return NULL; 841} 842 843static isl_bool isl_set_overlaps(__isl_keep isl_set *set1, 844 __isl_keep isl_set *set2) 845{ 846 return isl_bool_not(isl_set_is_disjoint(set1, set2)); 847} 848 849/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D 850 * and a dimension specification (Z^{n+1} -> Z^{n+1}), 851 * construct a map that is an overapproximation of the map 852 * that takes an element from the dom R \times Z to an 853 * element from ran R \times Z, such that the first n coordinates of the 854 * difference between them is a sum of differences between images 855 * and pre-images in one of the R_i and such that the last coordinate 856 * is equal to the number of steps taken. 857 * That is, let 858 * 859 * \Delta_i = { y - x | (x, y) in R_i } 860 * 861 * then the constructed map is an overapproximation of 862 * 863 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 864 * d = (\sum_i k_i \delta_i, \sum_i k_i) and 865 * x in dom R and x + d in ran R and 866 * \sum_i k_i >= 1 } 867 */ 868static __isl_give isl_map *construct_component(__isl_take isl_space *space, 869 __isl_keep isl_map *map, isl_bool *exact, int project) 870{ 871 struct isl_set *domain = NULL; 872 struct isl_set *range = NULL; 873 struct isl_map *app = NULL; 874 struct isl_map *path = NULL; 875 isl_bool overlaps; 876 int check; 877 878 domain = isl_map_domain(isl_map_copy(map)); 879 domain = isl_set_coalesce(domain); 880 range = isl_map_range(isl_map_copy(map)); 881 range = isl_set_coalesce(range); 882 overlaps = isl_set_overlaps(domain, range); 883 if (overlaps < 0 || !overlaps) { 884 isl_set_free(domain); 885 isl_set_free(range); 886 isl_space_free(space); 887 888 if (overlaps < 0) 889 map = NULL; 890 map = isl_map_copy(map); 891 map = isl_map_add_dims(map, isl_dim_in, 1); 892 map = isl_map_add_dims(map, isl_dim_out, 1); 893 map = set_path_length(map, 1, 1); 894 return map; 895 } 896 app = isl_map_from_domain_and_range(domain, range); 897 app = isl_map_add_dims(app, isl_dim_in, 1); 898 app = isl_map_add_dims(app, isl_dim_out, 1); 899 900 check = exact && *exact == isl_bool_true; 901 path = construct_extended_path(isl_space_copy(space), map, 902 check ? &project : NULL); 903 app = isl_map_intersect(app, path); 904 905 if (check && 906 (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app), 907 project)) < 0) 908 goto error; 909 910 isl_space_free(space); 911 app = set_path_length(app, 0, 1); 912 return app; 913error: 914 isl_space_free(space); 915 isl_map_free(app); 916 return NULL; 917} 918 919/* Call construct_component and, if "project" is set, project out 920 * the final coordinates. 921 */ 922static __isl_give isl_map *construct_projected_component( 923 __isl_take isl_space *space, 924 __isl_keep isl_map *map, isl_bool *exact, int project) 925{ 926 isl_map *app; 927 unsigned d; 928 929 if (!space) 930 return NULL; 931 d = isl_space_dim(space, isl_dim_in); 932 933 app = construct_component(space, map, exact, project); 934 if (project) { 935 app = isl_map_project_out(app, isl_dim_in, d - 1, 1); 936 app = isl_map_project_out(app, isl_dim_out, d - 1, 1); 937 } 938 return app; 939} 940 941/* Compute an extended version, i.e., with path lengths, of 942 * an overapproximation of the transitive closure of "bmap" 943 * with path lengths greater than or equal to zero and with 944 * domain and range equal to "dom". 945 */ 946static __isl_give isl_map *q_closure(__isl_take isl_space *space, 947 __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, 948 isl_bool *exact) 949{ 950 int project = 1; 951 isl_map *path; 952 isl_map *map; 953 isl_map *app; 954 955 dom = isl_set_add_dims(dom, isl_dim_set, 1); 956 app = isl_map_from_domain_and_range(dom, isl_set_copy(dom)); 957 map = isl_map_from_basic_map(isl_basic_map_copy(bmap)); 958 path = construct_extended_path(space, map, &project); 959 app = isl_map_intersect(app, path); 960 961 if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0) 962 goto error; 963 964 return app; 965error: 966 isl_map_free(app); 967 return NULL; 968} 969 970/* Check whether qc has any elements of length at least one 971 * with domain and/or range outside of dom and ran. 972 */ 973static isl_bool has_spurious_elements(__isl_keep isl_map *qc, 974 __isl_keep isl_set *dom, __isl_keep isl_set *ran) 975{ 976 isl_set *s; 977 isl_bool subset; 978 isl_size d; 979 980 d = isl_map_dim(qc, isl_dim_in); 981 if (d < 0 || !dom || !ran) 982 return isl_bool_error; 983 984 qc = isl_map_copy(qc); 985 qc = set_path_length(qc, 0, 1); 986 qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1); 987 qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1); 988 989 s = isl_map_domain(isl_map_copy(qc)); 990 subset = isl_set_is_subset(s, dom); 991 isl_set_free(s); 992 if (subset < 0) 993 goto error; 994 if (!subset) { 995 isl_map_free(qc); 996 return isl_bool_true; 997 } 998 999 s = isl_map_range(qc); 1000 subset = isl_set_is_subset(s, ran); 1001 isl_set_free(s); 1002 1003 return isl_bool_not(subset); 1004error: 1005 isl_map_free(qc); 1006 return isl_bool_error; 1007} 1008 1009#define LEFT 2 1010#define RIGHT 1 1011 1012/* For each basic map in "map", except i, check whether it combines 1013 * with the transitive closure that is reflexive on C combines 1014 * to the left and to the right. 1015 * 1016 * In particular, if 1017 * 1018 * dom map_j \subseteq C 1019 * 1020 * then right[j] is set to 1. Otherwise, if 1021 * 1022 * ran map_i \cap dom map_j = \emptyset 1023 * 1024 * then right[j] is set to 0. Otherwise, composing to the right 1025 * is impossible. 1026 * 1027 * Similar, for composing to the left, we have if 1028 * 1029 * ran map_j \subseteq C 1030 * 1031 * then left[j] is set to 1. Otherwise, if 1032 * 1033 * dom map_i \cap ran map_j = \emptyset 1034 * 1035 * then left[j] is set to 0. Otherwise, composing to the left 1036 * is impossible. 1037 * 1038 * The return value is or'd with LEFT if composing to the left 1039 * is possible and with RIGHT if composing to the right is possible. 1040 */ 1041static int composability(__isl_keep isl_set *C, int i, 1042 isl_set **dom, isl_set **ran, int *left, int *right, 1043 __isl_keep isl_map *map) 1044{ 1045 int j; 1046 int ok; 1047 1048 ok = LEFT | RIGHT; 1049 for (j = 0; j < map->n && ok; ++j) { 1050 isl_bool overlaps, subset; 1051 if (j == i) 1052 continue; 1053 1054 if (ok & RIGHT) { 1055 if (!dom[j]) 1056 dom[j] = isl_set_from_basic_set( 1057 isl_basic_map_domain( 1058 isl_basic_map_copy(map->p[j]))); 1059 if (!dom[j]) 1060 return -1; 1061 overlaps = isl_set_overlaps(ran[i], dom[j]); 1062 if (overlaps < 0) 1063 return -1; 1064 if (!overlaps) 1065 right[j] = 0; 1066 else { 1067 subset = isl_set_is_subset(dom[j], C); 1068 if (subset < 0) 1069 return -1; 1070 if (subset) 1071 right[j] = 1; 1072 else 1073 ok &= ~RIGHT; 1074 } 1075 } 1076 1077 if (ok & LEFT) { 1078 if (!ran[j]) 1079 ran[j] = isl_set_from_basic_set( 1080 isl_basic_map_range( 1081 isl_basic_map_copy(map->p[j]))); 1082 if (!ran[j]) 1083 return -1; 1084 overlaps = isl_set_overlaps(dom[i], ran[j]); 1085 if (overlaps < 0) 1086 return -1; 1087 if (!overlaps) 1088 left[j] = 0; 1089 else { 1090 subset = isl_set_is_subset(ran[j], C); 1091 if (subset < 0) 1092 return -1; 1093 if (subset) 1094 left[j] = 1; 1095 else 1096 ok &= ~LEFT; 1097 } 1098 } 1099 } 1100 1101 return ok; 1102} 1103 1104static __isl_give isl_map *anonymize(__isl_take isl_map *map) 1105{ 1106 map = isl_map_reset(map, isl_dim_in); 1107 map = isl_map_reset(map, isl_dim_out); 1108 return map; 1109} 1110 1111/* Return a map that is a union of the basic maps in "map", except i, 1112 * composed to left and right with qc based on the entries of "left" 1113 * and "right". 1114 */ 1115static __isl_give isl_map *compose(__isl_keep isl_map *map, int i, 1116 __isl_take isl_map *qc, int *left, int *right) 1117{ 1118 int j; 1119 isl_map *comp; 1120 1121 comp = isl_map_empty(isl_map_get_space(map)); 1122 for (j = 0; j < map->n; ++j) { 1123 isl_map *map_j; 1124 1125 if (j == i) 1126 continue; 1127 1128 map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j])); 1129 map_j = anonymize(map_j); 1130 if (left && left[j]) 1131 map_j = isl_map_apply_range(map_j, isl_map_copy(qc)); 1132 if (right && right[j]) 1133 map_j = isl_map_apply_range(isl_map_copy(qc), map_j); 1134 comp = isl_map_union(comp, map_j); 1135 } 1136 1137 comp = isl_map_compute_divs(comp); 1138 comp = isl_map_coalesce(comp); 1139 1140 isl_map_free(qc); 1141 1142 return comp; 1143} 1144 1145/* Compute the transitive closure of "map" incrementally by 1146 * computing 1147 * 1148 * map_i^+ \cup qc^+ 1149 * 1150 * or 1151 * 1152 * map_i^+ \cup ((id \cup map_i^) \circ qc^+) 1153 * 1154 * or 1155 * 1156 * map_i^+ \cup (qc^+ \circ (id \cup map_i^)) 1157 * 1158 * depending on whether left or right are NULL. 1159 */ 1160static __isl_give isl_map *compute_incremental( 1161 __isl_take isl_space *space, __isl_keep isl_map *map, 1162 int i, __isl_take isl_map *qc, int *left, int *right, isl_bool *exact) 1163{ 1164 isl_map *map_i; 1165 isl_map *tc; 1166 isl_map *rtc = NULL; 1167 1168 if (!map) 1169 goto error; 1170 isl_assert(map->ctx, left || right, goto error); 1171 1172 map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i])); 1173 tc = construct_projected_component(isl_space_copy(space), map_i, 1174 exact, 1); 1175 isl_map_free(map_i); 1176 1177 if (*exact) 1178 qc = isl_map_transitive_closure(qc, exact); 1179 1180 if (!*exact) { 1181 isl_space_free(space); 1182 isl_map_free(tc); 1183 isl_map_free(qc); 1184 return isl_map_universe(isl_map_get_space(map)); 1185 } 1186 1187 if (!left || !right) 1188 rtc = isl_map_union(isl_map_copy(tc), 1189 isl_map_identity(isl_map_get_space(tc))); 1190 if (!right) 1191 qc = isl_map_apply_range(rtc, qc); 1192 if (!left) 1193 qc = isl_map_apply_range(qc, rtc); 1194 qc = isl_map_union(tc, qc); 1195 1196 isl_space_free(space); 1197 1198 return qc; 1199error: 1200 isl_space_free(space); 1201 isl_map_free(qc); 1202 return NULL; 1203} 1204 1205/* Given a map "map", try to find a basic map such that 1206 * map^+ can be computed as 1207 * 1208 * map^+ = map_i^+ \cup 1209 * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+ 1210 * 1211 * with C the simple hull of the domain and range of the input map. 1212 * map_i^ \cup Id_C is computed by allowing the path lengths to be zero 1213 * and by intersecting domain and range with C. 1214 * Of course, we need to check that this is actually equal to map_i^ \cup Id_C. 1215 * Also, we only use the incremental computation if all the transitive 1216 * closures are exact and if the number of basic maps in the union, 1217 * after computing the integer divisions, is smaller than the number 1218 * of basic maps in the input map. 1219 */ 1220static isl_bool incremental_on_entire_domain(__isl_keep isl_space *space, 1221 __isl_keep isl_map *map, 1222 isl_set **dom, isl_set **ran, int *left, int *right, 1223 __isl_give isl_map **res) 1224{ 1225 int i; 1226 isl_set *C; 1227 isl_size d; 1228 1229 *res = NULL; 1230 1231 d = isl_map_dim(map, isl_dim_in); 1232 if (d < 0) 1233 return isl_bool_error; 1234 1235 C = isl_set_union(isl_map_domain(isl_map_copy(map)), 1236 isl_map_range(isl_map_copy(map))); 1237 C = isl_set_from_basic_set(isl_set_simple_hull(C)); 1238 if (!C) 1239 return isl_bool_error; 1240 if (C->n != 1) { 1241 isl_set_free(C); 1242 return isl_bool_false; 1243 } 1244 1245 for (i = 0; i < map->n; ++i) { 1246 isl_map *qc; 1247 isl_bool exact_i; 1248 isl_bool spurious; 1249 int j; 1250 dom[i] = isl_set_from_basic_set(isl_basic_map_domain( 1251 isl_basic_map_copy(map->p[i]))); 1252 ran[i] = isl_set_from_basic_set(isl_basic_map_range( 1253 isl_basic_map_copy(map->p[i]))); 1254 qc = q_closure(isl_space_copy(space), isl_set_copy(C), 1255 map->p[i], &exact_i); 1256 if (!qc) 1257 goto error; 1258 if (!exact_i) { 1259 isl_map_free(qc); 1260 continue; 1261 } 1262 spurious = has_spurious_elements(qc, dom[i], ran[i]); 1263 if (spurious) { 1264 isl_map_free(qc); 1265 if (spurious < 0) 1266 goto error; 1267 continue; 1268 } 1269 qc = isl_map_project_out(qc, isl_dim_in, d, 1); 1270 qc = isl_map_project_out(qc, isl_dim_out, d, 1); 1271 qc = isl_map_compute_divs(qc); 1272 for (j = 0; j < map->n; ++j) 1273 left[j] = right[j] = 1; 1274 qc = compose(map, i, qc, left, right); 1275 if (!qc) 1276 goto error; 1277 if (qc->n >= map->n) { 1278 isl_map_free(qc); 1279 continue; 1280 } 1281 *res = compute_incremental(isl_space_copy(space), map, i, qc, 1282 left, right, &exact_i); 1283 if (!*res) 1284 goto error; 1285 if (exact_i) 1286 break; 1287 isl_map_free(*res); 1288 *res = NULL; 1289 } 1290 1291 isl_set_free(C); 1292 1293 return isl_bool_ok(*res != NULL); 1294error: 1295 isl_set_free(C); 1296 return isl_bool_error; 1297} 1298 1299/* Try and compute the transitive closure of "map" as 1300 * 1301 * map^+ = map_i^+ \cup 1302 * \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+ 1303 * 1304 * with C either the simple hull of the domain and range of the entire 1305 * map or the simple hull of domain and range of map_i. 1306 */ 1307static __isl_give isl_map *incremental_closure(__isl_take isl_space *space, 1308 __isl_keep isl_map *map, isl_bool *exact, int project) 1309{ 1310 int i; 1311 isl_set **dom = NULL; 1312 isl_set **ran = NULL; 1313 int *left = NULL; 1314 int *right = NULL; 1315 isl_set *C; 1316 isl_size d; 1317 isl_map *res = NULL; 1318 1319 if (!project) 1320 return construct_projected_component(space, map, exact, 1321 project); 1322 1323 if (!map) 1324 goto error; 1325 if (map->n <= 1) 1326 return construct_projected_component(space, map, exact, 1327 project); 1328 1329 d = isl_map_dim(map, isl_dim_in); 1330 if (d < 0) 1331 goto error; 1332 1333 dom = isl_calloc_array(map->ctx, isl_set *, map->n); 1334 ran = isl_calloc_array(map->ctx, isl_set *, map->n); 1335 left = isl_calloc_array(map->ctx, int, map->n); 1336 right = isl_calloc_array(map->ctx, int, map->n); 1337 if (!ran || !dom || !left || !right) 1338 goto error; 1339 1340 if (incremental_on_entire_domain(space, map, dom, ran, left, right, 1341 &res) < 0) 1342 goto error; 1343 1344 for (i = 0; !res && i < map->n; ++i) { 1345 isl_map *qc; 1346 int comp; 1347 isl_bool exact_i, spurious; 1348 if (!dom[i]) 1349 dom[i] = isl_set_from_basic_set( 1350 isl_basic_map_domain( 1351 isl_basic_map_copy(map->p[i]))); 1352 if (!dom[i]) 1353 goto error; 1354 if (!ran[i]) 1355 ran[i] = isl_set_from_basic_set( 1356 isl_basic_map_range( 1357 isl_basic_map_copy(map->p[i]))); 1358 if (!ran[i]) 1359 goto error; 1360 C = isl_set_union(isl_set_copy(dom[i]), 1361 isl_set_copy(ran[i])); 1362 C = isl_set_from_basic_set(isl_set_simple_hull(C)); 1363 if (!C) 1364 goto error; 1365 if (C->n != 1) { 1366 isl_set_free(C); 1367 continue; 1368 } 1369 comp = composability(C, i, dom, ran, left, right, map); 1370 if (!comp || comp < 0) { 1371 isl_set_free(C); 1372 if (comp < 0) 1373 goto error; 1374 continue; 1375 } 1376 qc = q_closure(isl_space_copy(space), C, map->p[i], &exact_i); 1377 if (!qc) 1378 goto error; 1379 if (!exact_i) { 1380 isl_map_free(qc); 1381 continue; 1382 } 1383 spurious = has_spurious_elements(qc, dom[i], ran[i]); 1384 if (spurious) { 1385 isl_map_free(qc); 1386 if (spurious < 0) 1387 goto error; 1388 continue; 1389 } 1390 qc = isl_map_project_out(qc, isl_dim_in, d, 1); 1391 qc = isl_map_project_out(qc, isl_dim_out, d, 1); 1392 qc = isl_map_compute_divs(qc); 1393 qc = compose(map, i, qc, (comp & LEFT) ? left : NULL, 1394 (comp & RIGHT) ? right : NULL); 1395 if (!qc) 1396 goto error; 1397 if (qc->n >= map->n) { 1398 isl_map_free(qc); 1399 continue; 1400 } 1401 res = compute_incremental(isl_space_copy(space), map, i, qc, 1402 (comp & LEFT) ? left : NULL, 1403 (comp & RIGHT) ? right : NULL, &exact_i); 1404 if (!res) 1405 goto error; 1406 if (exact_i) 1407 break; 1408 isl_map_free(res); 1409 res = NULL; 1410 } 1411 1412 for (i = 0; i < map->n; ++i) { 1413 isl_set_free(dom[i]); 1414 isl_set_free(ran[i]); 1415 } 1416 free(dom); 1417 free(ran); 1418 free(left); 1419 free(right); 1420 1421 if (res) { 1422 isl_space_free(space); 1423 return res; 1424 } 1425 1426 return construct_projected_component(space, map, exact, project); 1427error: 1428 if (dom) 1429 for (i = 0; i < map->n; ++i) 1430 isl_set_free(dom[i]); 1431 free(dom); 1432 if (ran) 1433 for (i = 0; i < map->n; ++i) 1434 isl_set_free(ran[i]); 1435 free(ran); 1436 free(left); 1437 free(right); 1438 isl_space_free(space); 1439 return NULL; 1440} 1441 1442/* Given an array of sets "set", add "dom" at position "pos" 1443 * and search for elements at earlier positions that overlap with "dom". 1444 * If any can be found, then merge all of them, together with "dom", into 1445 * a single set and assign the union to the first in the array, 1446 * which becomes the new group leader for all groups involved in the merge. 1447 * During the search, we only consider group leaders, i.e., those with 1448 * group[i] = i, as the other sets have already been combined 1449 * with one of the group leaders. 1450 */ 1451static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos) 1452{ 1453 int i; 1454 1455 group[pos] = pos; 1456 set[pos] = isl_set_copy(dom); 1457 1458 for (i = pos - 1; i >= 0; --i) { 1459 isl_bool o; 1460 1461 if (group[i] != i) 1462 continue; 1463 1464 o = isl_set_overlaps(set[i], dom); 1465 if (o < 0) 1466 goto error; 1467 if (!o) 1468 continue; 1469 1470 set[i] = isl_set_union(set[i], set[group[pos]]); 1471 set[group[pos]] = NULL; 1472 if (!set[i]) 1473 goto error; 1474 group[group[pos]] = i; 1475 group[pos] = i; 1476 } 1477 1478 isl_set_free(dom); 1479 return 0; 1480error: 1481 isl_set_free(dom); 1482 return -1; 1483} 1484 1485/* Construct a map [x] -> [x+1], with parameters prescribed by "space". 1486 */ 1487static __isl_give isl_map *increment(__isl_take isl_space *space) 1488{ 1489 int k; 1490 isl_basic_map *bmap; 1491 isl_size total; 1492 1493 space = isl_space_set_from_params(space); 1494 space = isl_space_add_dims(space, isl_dim_set, 1); 1495 space = isl_space_map_from_set(space); 1496 bmap = isl_basic_map_alloc_space(space, 0, 1, 0); 1497 total = isl_basic_map_dim(bmap, isl_dim_all); 1498 k = isl_basic_map_alloc_equality(bmap); 1499 if (total < 0 || k < 0) 1500 goto error; 1501 isl_seq_clr(bmap->eq[k], 1 + total); 1502 isl_int_set_si(bmap->eq[k][0], 1); 1503 isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1); 1504 isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1); 1505 return isl_map_from_basic_map(bmap); 1506error: 1507 isl_basic_map_free(bmap); 1508 return NULL; 1509} 1510 1511/* Replace each entry in the n by n grid of maps by the cross product 1512 * with the relation { [i] -> [i + 1] }. 1513 */ 1514static isl_stat add_length(__isl_keep isl_map *map, isl_map ***grid, int n) 1515{ 1516 int i, j; 1517 isl_space *space; 1518 isl_map *step; 1519 1520 space = isl_space_params(isl_map_get_space(map)); 1521 step = increment(space); 1522 1523 if (!step) 1524 return isl_stat_error; 1525 1526 for (i = 0; i < n; ++i) 1527 for (j = 0; j < n; ++j) 1528 grid[i][j] = isl_map_product(grid[i][j], 1529 isl_map_copy(step)); 1530 1531 isl_map_free(step); 1532 1533 return isl_stat_ok; 1534} 1535 1536/* The core of the Floyd-Warshall algorithm. 1537 * Updates the given n x x matrix of relations in place. 1538 * 1539 * The algorithm iterates over all vertices. In each step, the whole 1540 * matrix is updated to include all paths that go to the current vertex, 1541 * possibly stay there a while (including passing through earlier vertices) 1542 * and then come back. At the start of each iteration, the diagonal 1543 * element corresponding to the current vertex is replaced by its 1544 * transitive closure to account for all indirect paths that stay 1545 * in the current vertex. 1546 */ 1547static void floyd_warshall_iterate(isl_map ***grid, int n, isl_bool *exact) 1548{ 1549 int r, p, q; 1550 1551 for (r = 0; r < n; ++r) { 1552 isl_bool r_exact; 1553 int check = exact && *exact == isl_bool_true; 1554 grid[r][r] = isl_map_transitive_closure(grid[r][r], 1555 check ? &r_exact : NULL); 1556 if (check && !r_exact) 1557 *exact = isl_bool_false; 1558 1559 for (p = 0; p < n; ++p) 1560 for (q = 0; q < n; ++q) { 1561 isl_map *loop; 1562 if (p == r && q == r) 1563 continue; 1564 loop = isl_map_apply_range( 1565 isl_map_copy(grid[p][r]), 1566 isl_map_copy(grid[r][q])); 1567 grid[p][q] = isl_map_union(grid[p][q], loop); 1568 loop = isl_map_apply_range( 1569 isl_map_copy(grid[p][r]), 1570 isl_map_apply_range( 1571 isl_map_copy(grid[r][r]), 1572 isl_map_copy(grid[r][q]))); 1573 grid[p][q] = isl_map_union(grid[p][q], loop); 1574 grid[p][q] = isl_map_coalesce(grid[p][q]); 1575 } 1576 } 1577} 1578 1579/* Given a partition of the domains and ranges of the basic maps in "map", 1580 * apply the Floyd-Warshall algorithm with the elements in the partition 1581 * as vertices. 1582 * 1583 * In particular, there are "n" elements in the partition and "group" is 1584 * an array of length 2 * map->n with entries in [0,n-1]. 1585 * 1586 * We first construct a matrix of relations based on the partition information, 1587 * apply Floyd-Warshall on this matrix of relations and then take the 1588 * union of all entries in the matrix as the final result. 1589 * 1590 * If we are actually computing the power instead of the transitive closure, 1591 * i.e., when "project" is not set, then the result should have the 1592 * path lengths encoded as the difference between an extra pair of 1593 * coordinates. We therefore apply the nested transitive closures 1594 * to relations that include these lengths. In particular, we replace 1595 * the input relation by the cross product with the unit length relation 1596 * { [i] -> [i + 1] }. 1597 */ 1598static __isl_give isl_map *floyd_warshall_with_groups( 1599 __isl_take isl_space *space, __isl_keep isl_map *map, 1600 isl_bool *exact, int project, int *group, int n) 1601{ 1602 int i, j, k; 1603 isl_map ***grid = NULL; 1604 isl_map *app; 1605 1606 if (!map) 1607 goto error; 1608 1609 if (n == 1) { 1610 free(group); 1611 return incremental_closure(space, map, exact, project); 1612 } 1613 1614 grid = isl_calloc_array(map->ctx, isl_map **, n); 1615 if (!grid) 1616 goto error; 1617 for (i = 0; i < n; ++i) { 1618 grid[i] = isl_calloc_array(map->ctx, isl_map *, n); 1619 if (!grid[i]) 1620 goto error; 1621 for (j = 0; j < n; ++j) 1622 grid[i][j] = isl_map_empty(isl_map_get_space(map)); 1623 } 1624 1625 for (k = 0; k < map->n; ++k) { 1626 i = group[2 * k]; 1627 j = group[2 * k + 1]; 1628 grid[i][j] = isl_map_union(grid[i][j], 1629 isl_map_from_basic_map( 1630 isl_basic_map_copy(map->p[k]))); 1631 } 1632 1633 if (!project && add_length(map, grid, n) < 0) 1634 goto error; 1635 1636 floyd_warshall_iterate(grid, n, exact); 1637 1638 app = isl_map_empty(isl_map_get_space(grid[0][0])); 1639 1640 for (i = 0; i < n; ++i) { 1641 for (j = 0; j < n; ++j) 1642 app = isl_map_union(app, grid[i][j]); 1643 free(grid[i]); 1644 } 1645 free(grid); 1646 1647 free(group); 1648 isl_space_free(space); 1649 1650 return app; 1651error: 1652 if (grid) 1653 for (i = 0; i < n; ++i) { 1654 if (!grid[i]) 1655 continue; 1656 for (j = 0; j < n; ++j) 1657 isl_map_free(grid[i][j]); 1658 free(grid[i]); 1659 } 1660 free(grid); 1661 free(group); 1662 isl_space_free(space); 1663 return NULL; 1664} 1665 1666/* Partition the domains and ranges of the n basic relations in list 1667 * into disjoint cells. 1668 * 1669 * To find the partition, we simply consider all of the domains 1670 * and ranges in turn and combine those that overlap. 1671 * "set" contains the partition elements and "group" indicates 1672 * to which partition element a given domain or range belongs. 1673 * The domain of basic map i corresponds to element 2 * i in these arrays, 1674 * while the domain corresponds to element 2 * i + 1. 1675 * During the construction group[k] is either equal to k, 1676 * in which case set[k] contains the union of all the domains and 1677 * ranges in the corresponding group, or is equal to some l < k, 1678 * with l another domain or range in the same group. 1679 */ 1680static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n, 1681 isl_set ***set, int *n_group) 1682{ 1683 int i; 1684 int *group = NULL; 1685 int g; 1686 1687 *set = isl_calloc_array(ctx, isl_set *, 2 * n); 1688 group = isl_alloc_array(ctx, int, 2 * n); 1689 1690 if (!*set || !group) 1691 goto error; 1692 1693 for (i = 0; i < n; ++i) { 1694 isl_set *dom; 1695 dom = isl_set_from_basic_set(isl_basic_map_domain( 1696 isl_basic_map_copy(list[i]))); 1697 if (merge(*set, group, dom, 2 * i) < 0) 1698 goto error; 1699 dom = isl_set_from_basic_set(isl_basic_map_range( 1700 isl_basic_map_copy(list[i]))); 1701 if (merge(*set, group, dom, 2 * i + 1) < 0) 1702 goto error; 1703 } 1704 1705 g = 0; 1706 for (i = 0; i < 2 * n; ++i) 1707 if (group[i] == i) { 1708 if (g != i) { 1709 (*set)[g] = (*set)[i]; 1710 (*set)[i] = NULL; 1711 } 1712 group[i] = g++; 1713 } else 1714 group[i] = group[group[i]]; 1715 1716 *n_group = g; 1717 1718 return group; 1719error: 1720 if (*set) { 1721 for (i = 0; i < 2 * n; ++i) 1722 isl_set_free((*set)[i]); 1723 free(*set); 1724 *set = NULL; 1725 } 1726 free(group); 1727 return NULL; 1728} 1729 1730/* Check if the domains and ranges of the basic maps in "map" can 1731 * be partitioned, and if so, apply Floyd-Warshall on the elements 1732 * of the partition. Note that we also apply this algorithm 1733 * if we want to compute the power, i.e., when "project" is not set. 1734 * However, the results are unlikely to be exact since the recursive 1735 * calls inside the Floyd-Warshall algorithm typically result in 1736 * non-linear path lengths quite quickly. 1737 */ 1738static __isl_give isl_map *floyd_warshall(__isl_take isl_space *space, 1739 __isl_keep isl_map *map, isl_bool *exact, int project) 1740{ 1741 int i; 1742 isl_set **set = NULL; 1743 int *group = NULL; 1744 int n; 1745 1746 if (!map) 1747 goto error; 1748 if (map->n <= 1) 1749 return incremental_closure(space, map, exact, project); 1750 1751 group = setup_groups(map->ctx, map->p, map->n, &set, &n); 1752 if (!group) 1753 goto error; 1754 1755 for (i = 0; i < 2 * map->n; ++i) 1756 isl_set_free(set[i]); 1757 1758 free(set); 1759 1760 return floyd_warshall_with_groups(space, map, exact, project, group, n); 1761error: 1762 isl_space_free(space); 1763 return NULL; 1764} 1765 1766/* Structure for representing the nodes of the graph of which 1767 * strongly connected components are being computed. 1768 * 1769 * list contains the actual nodes 1770 * check_closed is set if we may have used the fact that 1771 * a pair of basic maps can be interchanged 1772 */ 1773struct isl_tc_follows_data { 1774 isl_basic_map **list; 1775 int check_closed; 1776}; 1777 1778/* Check whether in the computation of the transitive closure 1779 * "list[i]" (R_1) should follow (or be part of the same component as) 1780 * "list[j]" (R_2). 1781 * 1782 * That is check whether 1783 * 1784 * R_1 \circ R_2 1785 * 1786 * is a subset of 1787 * 1788 * R_2 \circ R_1 1789 * 1790 * If so, then there is no reason for R_1 to immediately follow R_2 1791 * in any path. 1792 * 1793 * *check_closed is set if the subset relation holds while 1794 * R_1 \circ R_2 is not empty. 1795 */ 1796static isl_bool basic_map_follows(int i, int j, void *user) 1797{ 1798 struct isl_tc_follows_data *data = user; 1799 struct isl_map *map12 = NULL; 1800 struct isl_map *map21 = NULL; 1801 isl_bool applies, subset; 1802 1803 applies = isl_basic_map_applies_range(data->list[j], data->list[i]); 1804 if (applies < 0) 1805 return isl_bool_error; 1806 if (!applies) 1807 return isl_bool_false; 1808 1809 map21 = isl_map_from_basic_map( 1810 isl_basic_map_apply_range( 1811 isl_basic_map_copy(data->list[j]), 1812 isl_basic_map_copy(data->list[i]))); 1813 subset = isl_map_is_empty(map21); 1814 if (subset < 0) 1815 goto error; 1816 if (subset) { 1817 isl_map_free(map21); 1818 return isl_bool_false; 1819 } 1820 1821 if (!isl_basic_map_is_transformation(data->list[i]) || 1822 !isl_basic_map_is_transformation(data->list[j])) { 1823 isl_map_free(map21); 1824 return isl_bool_true; 1825 } 1826 1827 map12 = isl_map_from_basic_map( 1828 isl_basic_map_apply_range( 1829 isl_basic_map_copy(data->list[i]), 1830 isl_basic_map_copy(data->list[j]))); 1831 1832 subset = isl_map_is_subset(map21, map12); 1833 1834 isl_map_free(map12); 1835 isl_map_free(map21); 1836 1837 if (subset) 1838 data->check_closed = 1; 1839 1840 return isl_bool_not(subset); 1841error: 1842 isl_map_free(map21); 1843 return isl_bool_error; 1844} 1845 1846/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D 1847 * and a dimension specification (Z^{n+1} -> Z^{n+1}), 1848 * construct a map that is an overapproximation of the map 1849 * that takes an element from the dom R \times Z to an 1850 * element from ran R \times Z, such that the first n coordinates of the 1851 * difference between them is a sum of differences between images 1852 * and pre-images in one of the R_i and such that the last coordinate 1853 * is equal to the number of steps taken. 1854 * If "project" is set, then these final coordinates are not included, 1855 * i.e., a relation of type Z^n -> Z^n is returned. 1856 * That is, let 1857 * 1858 * \Delta_i = { y - x | (x, y) in R_i } 1859 * 1860 * then the constructed map is an overapproximation of 1861 * 1862 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1863 * d = (\sum_i k_i \delta_i, \sum_i k_i) and 1864 * x in dom R and x + d in ran R } 1865 * 1866 * or 1867 * 1868 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1869 * d = (\sum_i k_i \delta_i) and 1870 * x in dom R and x + d in ran R } 1871 * 1872 * if "project" is set. 1873 * 1874 * We first split the map into strongly connected components, perform 1875 * the above on each component and then join the results in the correct 1876 * order, at each join also taking in the union of both arguments 1877 * to allow for paths that do not go through one of the two arguments. 1878 */ 1879static __isl_give isl_map *construct_power_components( 1880 __isl_take isl_space *space, __isl_keep isl_map *map, isl_bool *exact, 1881 int project) 1882{ 1883 int i, n, c; 1884 struct isl_map *path = NULL; 1885 struct isl_tc_follows_data data; 1886 struct isl_tarjan_graph *g = NULL; 1887 isl_bool *orig_exact; 1888 isl_bool local_exact; 1889 1890 if (!map) 1891 goto error; 1892 if (map->n <= 1) 1893 return floyd_warshall(space, map, exact, project); 1894 1895 data.list = map->p; 1896 data.check_closed = 0; 1897 g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data); 1898 if (!g) 1899 goto error; 1900 1901 orig_exact = exact; 1902 if (data.check_closed && !exact) 1903 exact = &local_exact; 1904 1905 c = 0; 1906 i = 0; 1907 n = map->n; 1908 if (project) 1909 path = isl_map_empty(isl_map_get_space(map)); 1910 else 1911 path = isl_map_empty(isl_space_copy(space)); 1912 path = anonymize(path); 1913 while (n) { 1914 struct isl_map *comp; 1915 isl_map *path_comp, *path_comb; 1916 comp = isl_map_alloc_space(isl_map_get_space(map), n, 0); 1917 while (g->order[i] != -1) { 1918 comp = isl_map_add_basic_map(comp, 1919 isl_basic_map_copy(map->p[g->order[i]])); 1920 --n; 1921 ++i; 1922 } 1923 path_comp = floyd_warshall(isl_space_copy(space), 1924 comp, exact, project); 1925 path_comp = anonymize(path_comp); 1926 path_comb = isl_map_apply_range(isl_map_copy(path), 1927 isl_map_copy(path_comp)); 1928 path = isl_map_union(path, path_comp); 1929 path = isl_map_union(path, path_comb); 1930 isl_map_free(comp); 1931 ++i; 1932 ++c; 1933 } 1934 1935 if (c > 1 && data.check_closed && !*exact) { 1936 isl_bool closed; 1937 1938 closed = isl_map_is_transitively_closed(path); 1939 if (closed < 0) 1940 goto error; 1941 if (!closed) { 1942 isl_tarjan_graph_free(g); 1943 isl_map_free(path); 1944 return floyd_warshall(space, map, orig_exact, project); 1945 } 1946 } 1947 1948 isl_tarjan_graph_free(g); 1949 isl_space_free(space); 1950 1951 return path; 1952error: 1953 isl_tarjan_graph_free(g); 1954 isl_space_free(space); 1955 isl_map_free(path); 1956 return NULL; 1957} 1958 1959/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D, 1960 * construct a map that is an overapproximation of the map 1961 * that takes an element from the space D to another 1962 * element from the same space, such that the difference between 1963 * them is a strictly positive sum of differences between images 1964 * and pre-images in one of the R_i. 1965 * The number of differences in the sum is equated to parameter "param". 1966 * That is, let 1967 * 1968 * \Delta_i = { y - x | (x, y) in R_i } 1969 * 1970 * then the constructed map is an overapproximation of 1971 * 1972 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1973 * d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 } 1974 * or 1975 * 1976 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i : 1977 * d = \sum_i k_i \delta_i and \sum_i k_i > 0 } 1978 * 1979 * if "project" is set. 1980 * 1981 * If "project" is not set, then 1982 * we construct an extended mapping with an extra coordinate 1983 * that indicates the number of steps taken. In particular, 1984 * the difference in the last coordinate is equal to the number 1985 * of steps taken to move from a domain element to the corresponding 1986 * image element(s). 1987 */ 1988static __isl_give isl_map *construct_power(__isl_keep isl_map *map, 1989 isl_bool *exact, int project) 1990{ 1991 struct isl_map *app = NULL; 1992 isl_space *space = NULL; 1993 1994 if (!map) 1995 return NULL; 1996 1997 space = isl_map_get_space(map); 1998 1999 space = isl_space_add_dims(space, isl_dim_in, 1); 2000 space = isl_space_add_dims(space, isl_dim_out, 1); 2001 2002 app = construct_power_components(isl_space_copy(space), map, 2003 exact, project); 2004 2005 isl_space_free(space); 2006 2007 return app; 2008} 2009 2010/* Compute the positive powers of "map", or an overapproximation. 2011 * If the result is exact, then *exact is set to 1. 2012 * 2013 * If project is set, then we are actually interested in the transitive 2014 * closure, so we can use a more relaxed exactness check. 2015 * The lengths of the paths are also projected out instead of being 2016 * encoded as the difference between an extra pair of final coordinates. 2017 */ 2018static __isl_give isl_map *map_power(__isl_take isl_map *map, 2019 isl_bool *exact, int project) 2020{ 2021 struct isl_map *app = NULL; 2022 2023 if (exact) 2024 *exact = isl_bool_true; 2025 2026 if (isl_map_check_transformation(map) < 0) 2027 return isl_map_free(map); 2028 2029 app = construct_power(map, exact, project); 2030 2031 isl_map_free(map); 2032 return app; 2033} 2034 2035/* Compute the positive powers of "map", or an overapproximation. 2036 * The result maps the exponent to a nested copy of the corresponding power. 2037 * If the result is exact, then *exact is set to 1. 2038 * map_power constructs an extended relation with the path lengths 2039 * encoded as the difference between the final coordinates. 2040 * In the final step, this difference is equated to an extra parameter 2041 * and made positive. The extra coordinates are subsequently projected out 2042 * and the parameter is turned into the domain of the result. 2043 */ 2044__isl_give isl_map *isl_map_power(__isl_take isl_map *map, isl_bool *exact) 2045{ 2046 isl_space *target_space; 2047 isl_space *space; 2048 isl_map *diff; 2049 isl_size d; 2050 isl_size param; 2051 2052 d = isl_map_dim(map, isl_dim_in); 2053 param = isl_map_dim(map, isl_dim_param); 2054 if (d < 0 || param < 0) 2055 return isl_map_free(map); 2056 2057 map = isl_map_compute_divs(map); 2058 map = isl_map_coalesce(map); 2059 2060 if (isl_map_plain_is_empty(map)) { 2061 map = isl_map_from_range(isl_map_wrap(map)); 2062 map = isl_map_add_dims(map, isl_dim_in, 1); 2063 map = isl_map_set_dim_name(map, isl_dim_in, 0, "k"); 2064 return map; 2065 } 2066 2067 target_space = isl_map_get_space(map); 2068 target_space = isl_space_from_range(isl_space_wrap(target_space)); 2069 target_space = isl_space_add_dims(target_space, isl_dim_in, 1); 2070 target_space = isl_space_set_dim_name(target_space, isl_dim_in, 0, "k"); 2071 2072 map = map_power(map, exact, 0); 2073 2074 map = isl_map_add_dims(map, isl_dim_param, 1); 2075 space = isl_map_get_space(map); 2076 diff = equate_parameter_to_length(space, param); 2077 map = isl_map_intersect(map, diff); 2078 map = isl_map_project_out(map, isl_dim_in, d, 1); 2079 map = isl_map_project_out(map, isl_dim_out, d, 1); 2080 map = isl_map_from_range(isl_map_wrap(map)); 2081 map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1); 2082 2083 map = isl_map_reset_space(map, target_space); 2084 2085 return map; 2086} 2087 2088/* Compute a relation that maps each element in the range of the input 2089 * relation to the lengths of all paths composed of edges in the input 2090 * relation that end up in the given range element. 2091 * The result may be an overapproximation, in which case *exact is set to 0. 2092 * The resulting relation is very similar to the power relation. 2093 * The difference are that the domain has been projected out, the 2094 * range has become the domain and the exponent is the range instead 2095 * of a parameter. 2096 */ 2097__isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map, 2098 isl_bool *exact) 2099{ 2100 isl_space *space; 2101 isl_map *diff; 2102 isl_size d; 2103 isl_size param; 2104 2105 d = isl_map_dim(map, isl_dim_in); 2106 param = isl_map_dim(map, isl_dim_param); 2107 if (d < 0 || param < 0) 2108 return isl_map_free(map); 2109 2110 map = isl_map_compute_divs(map); 2111 map = isl_map_coalesce(map); 2112 2113 if (isl_map_plain_is_empty(map)) { 2114 if (exact) 2115 *exact = isl_bool_true; 2116 map = isl_map_project_out(map, isl_dim_out, 0, d); 2117 map = isl_map_add_dims(map, isl_dim_out, 1); 2118 return map; 2119 } 2120 2121 map = map_power(map, exact, 0); 2122 2123 map = isl_map_add_dims(map, isl_dim_param, 1); 2124 space = isl_map_get_space(map); 2125 diff = equate_parameter_to_length(space, param); 2126 map = isl_map_intersect(map, diff); 2127 map = isl_map_project_out(map, isl_dim_in, 0, d + 1); 2128 map = isl_map_project_out(map, isl_dim_out, d, 1); 2129 map = isl_map_reverse(map); 2130 map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1); 2131 2132 return map; 2133} 2134 2135/* Given a map, compute the smallest superset of this map that is of the form 2136 * 2137 * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } 2138 * 2139 * (where p ranges over the (non-parametric) dimensions), 2140 * compute the transitive closure of this map, i.e., 2141 * 2142 * { i -> j : exists k > 0: 2143 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2144 * 2145 * and intersect domain and range of this transitive closure with 2146 * the given domain and range. 2147 * 2148 * If with_id is set, then try to include as much of the identity mapping 2149 * as possible, by computing 2150 * 2151 * { i -> j : exists k >= 0: 2152 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2153 * 2154 * instead (i.e., allow k = 0). 2155 * 2156 * In practice, we compute the difference set 2157 * 2158 * delta = { j - i | i -> j in map }, 2159 * 2160 * look for stride constraint on the individual dimensions and compute 2161 * (constant) lower and upper bounds for each individual dimension, 2162 * adding a constraint for each bound not equal to infinity. 2163 */ 2164static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map, 2165 __isl_take isl_set *dom, __isl_take isl_set *ran, int with_id) 2166{ 2167 int i; 2168 int k; 2169 unsigned d; 2170 unsigned nparam; 2171 unsigned total; 2172 isl_space *space; 2173 isl_set *delta; 2174 isl_map *app = NULL; 2175 isl_basic_set *aff = NULL; 2176 isl_basic_map *bmap = NULL; 2177 isl_vec *obj = NULL; 2178 isl_int opt; 2179 2180 isl_int_init(opt); 2181 2182 delta = isl_map_deltas(isl_map_copy(map)); 2183 2184 aff = isl_set_affine_hull(isl_set_copy(delta)); 2185 if (!aff) 2186 goto error; 2187 space = isl_map_get_space(map); 2188 d = isl_space_dim(space, isl_dim_in); 2189 nparam = isl_space_dim(space, isl_dim_param); 2190 total = isl_space_dim(space, isl_dim_all); 2191 bmap = isl_basic_map_alloc_space(space, 2192 aff->n_div + 1, aff->n_div, 2 * d + 1); 2193 for (i = 0; i < aff->n_div + 1; ++i) { 2194 k = isl_basic_map_alloc_div(bmap); 2195 if (k < 0) 2196 goto error; 2197 isl_int_set_si(bmap->div[k][0], 0); 2198 } 2199 for (i = 0; i < aff->n_eq; ++i) { 2200 if (!isl_basic_set_eq_is_stride(aff, i)) 2201 continue; 2202 k = isl_basic_map_alloc_equality(bmap); 2203 if (k < 0) 2204 goto error; 2205 isl_seq_clr(bmap->eq[k], 1 + nparam); 2206 isl_seq_cpy(bmap->eq[k] + 1 + nparam + d, 2207 aff->eq[i] + 1 + nparam, d); 2208 isl_seq_neg(bmap->eq[k] + 1 + nparam, 2209 aff->eq[i] + 1 + nparam, d); 2210 isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d, 2211 aff->eq[i] + 1 + nparam + d, aff->n_div); 2212 isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0); 2213 } 2214 obj = isl_vec_alloc(map->ctx, 1 + nparam + d); 2215 if (!obj) 2216 goto error; 2217 isl_seq_clr(obj->el, 1 + nparam + d); 2218 for (i = 0; i < d; ++ i) { 2219 enum isl_lp_result res; 2220 2221 isl_int_set_si(obj->el[1 + nparam + i], 1); 2222 2223 res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt, 2224 NULL, NULL); 2225 if (res == isl_lp_error) 2226 goto error; 2227 if (res == isl_lp_ok) { 2228 k = isl_basic_map_alloc_inequality(bmap); 2229 if (k < 0) 2230 goto error; 2231 isl_seq_clr(bmap->ineq[k], 2232 1 + nparam + 2 * d + bmap->n_div); 2233 isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1); 2234 isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1); 2235 isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt); 2236 } 2237 2238 res = isl_set_solve_lp(delta, 1, 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_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt); 2251 } 2252 2253 isl_int_set_si(obj->el[1 + nparam + i], 0); 2254 } 2255 k = isl_basic_map_alloc_inequality(bmap); 2256 if (k < 0) 2257 goto error; 2258 isl_seq_clr(bmap->ineq[k], 2259 1 + nparam + 2 * d + bmap->n_div); 2260 if (!with_id) 2261 isl_int_set_si(bmap->ineq[k][0], -1); 2262 isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1); 2263 2264 app = isl_map_from_domain_and_range(dom, ran); 2265 2266 isl_vec_free(obj); 2267 isl_basic_set_free(aff); 2268 isl_map_free(map); 2269 bmap = isl_basic_map_finalize(bmap); 2270 isl_set_free(delta); 2271 isl_int_clear(opt); 2272 2273 map = isl_map_from_basic_map(bmap); 2274 map = isl_map_intersect(map, app); 2275 2276 return map; 2277error: 2278 isl_vec_free(obj); 2279 isl_basic_map_free(bmap); 2280 isl_basic_set_free(aff); 2281 isl_set_free(dom); 2282 isl_set_free(ran); 2283 isl_map_free(map); 2284 isl_set_free(delta); 2285 isl_int_clear(opt); 2286 return NULL; 2287} 2288 2289/* Given a map, compute the smallest superset of this map that is of the form 2290 * 2291 * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } 2292 * 2293 * (where p ranges over the (non-parametric) dimensions), 2294 * compute the transitive closure of this map, i.e., 2295 * 2296 * { i -> j : exists k > 0: 2297 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2298 * 2299 * and intersect domain and range of this transitive closure with 2300 * domain and range of the original map. 2301 */ 2302static __isl_give isl_map *box_closure(__isl_take isl_map *map) 2303{ 2304 isl_set *domain; 2305 isl_set *range; 2306 2307 domain = isl_map_domain(isl_map_copy(map)); 2308 domain = isl_set_coalesce(domain); 2309 range = isl_map_range(isl_map_copy(map)); 2310 range = isl_set_coalesce(range); 2311 2312 return box_closure_on_domain(map, domain, range, 0); 2313} 2314 2315/* Given a map, compute the smallest superset of this map that is of the form 2316 * 2317 * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } 2318 * 2319 * (where p ranges over the (non-parametric) dimensions), 2320 * compute the transitive and partially reflexive closure of this map, i.e., 2321 * 2322 * { i -> j : exists k >= 0: 2323 * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } 2324 * 2325 * and intersect domain and range of this transitive closure with 2326 * the given domain. 2327 */ 2328static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map, 2329 __isl_take isl_set *dom) 2330{ 2331 return box_closure_on_domain(map, dom, isl_set_copy(dom), 1); 2332} 2333 2334/* Check whether app is the transitive closure of map. 2335 * In particular, check that app is acyclic and, if so, 2336 * check that 2337 * 2338 * app \subset (map \cup (map \circ app)) 2339 */ 2340static isl_bool check_exactness_omega(__isl_keep isl_map *map, 2341 __isl_keep isl_map *app) 2342{ 2343 isl_set *delta; 2344 int i; 2345 isl_bool is_empty, is_exact; 2346 isl_size d; 2347 isl_map *test; 2348 2349 delta = isl_map_deltas(isl_map_copy(app)); 2350 d = isl_set_dim(delta, isl_dim_set); 2351 if (d < 0) 2352 delta = isl_set_free(delta); 2353 for (i = 0; i < d; ++i) 2354 delta = isl_set_fix_si(delta, isl_dim_set, i, 0); 2355 is_empty = isl_set_is_empty(delta); 2356 isl_set_free(delta); 2357 if (is_empty < 0 || !is_empty) 2358 return is_empty; 2359 2360 test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map)); 2361 test = isl_map_union(test, isl_map_copy(map)); 2362 is_exact = isl_map_is_subset(app, test); 2363 isl_map_free(test); 2364 2365 return is_exact; 2366} 2367 2368/* Check if basic map M_i can be combined with all the other 2369 * basic maps such that 2370 * 2371 * (\cup_j M_j)^+ 2372 * 2373 * can be computed as 2374 * 2375 * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+ 2376 * 2377 * In particular, check if we can compute a compact representation 2378 * of 2379 * 2380 * M_i^* \circ M_j \circ M_i^* 2381 * 2382 * for each j != i. 2383 * Let M_i^? be an extension of M_i^+ that allows paths 2384 * of length zero, i.e., the result of box_closure(., 1). 2385 * The criterion, as proposed by Kelly et al., is that 2386 * id = M_i^? - M_i^+ can be represented as a basic map 2387 * and that 2388 * 2389 * id \circ M_j \circ id = M_j 2390 * 2391 * for each j != i. 2392 * 2393 * If this function returns 1, then tc and qc are set to 2394 * M_i^+ and M_i^?, respectively. 2395 */ 2396static int can_be_split_off(__isl_keep isl_map *map, int i, 2397 __isl_give isl_map **tc, __isl_give isl_map **qc) 2398{ 2399 isl_map *map_i, *id = NULL; 2400 int j = -1; 2401 isl_set *C; 2402 2403 *tc = NULL; 2404 *qc = NULL; 2405 2406 C = isl_set_union(isl_map_domain(isl_map_copy(map)), 2407 isl_map_range(isl_map_copy(map))); 2408 C = isl_set_from_basic_set(isl_set_simple_hull(C)); 2409 if (!C) 2410 goto error; 2411 2412 map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i])); 2413 *tc = box_closure(isl_map_copy(map_i)); 2414 *qc = box_closure_with_identity(map_i, C); 2415 id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc)); 2416 2417 if (!id || !*qc) 2418 goto error; 2419 if (id->n != 1 || (*qc)->n != 1) 2420 goto done; 2421 2422 for (j = 0; j < map->n; ++j) { 2423 isl_map *map_j, *test; 2424 int is_ok; 2425 2426 if (i == j) 2427 continue; 2428 map_j = isl_map_from_basic_map( 2429 isl_basic_map_copy(map->p[j])); 2430 test = isl_map_apply_range(isl_map_copy(id), 2431 isl_map_copy(map_j)); 2432 test = isl_map_apply_range(test, isl_map_copy(id)); 2433 is_ok = isl_map_is_equal(test, map_j); 2434 isl_map_free(map_j); 2435 isl_map_free(test); 2436 if (is_ok < 0) 2437 goto error; 2438 if (!is_ok) 2439 break; 2440 } 2441 2442done: 2443 isl_map_free(id); 2444 if (j == map->n) 2445 return 1; 2446 2447 isl_map_free(*qc); 2448 isl_map_free(*tc); 2449 *qc = NULL; 2450 *tc = NULL; 2451 2452 return 0; 2453error: 2454 isl_map_free(id); 2455 isl_map_free(*qc); 2456 isl_map_free(*tc); 2457 *qc = NULL; 2458 *tc = NULL; 2459 return -1; 2460} 2461 2462static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map, 2463 isl_bool *exact) 2464{ 2465 isl_map *app; 2466 2467 app = box_closure(isl_map_copy(map)); 2468 if (exact) { 2469 isl_bool is_exact = check_exactness_omega(map, app); 2470 2471 if (is_exact < 0) 2472 app = isl_map_free(app); 2473 else 2474 *exact = is_exact; 2475 } 2476 2477 isl_map_free(map); 2478 return app; 2479} 2480 2481/* Compute an overapproximation of the transitive closure of "map" 2482 * using a variation of the algorithm from 2483 * "Transitive Closure of Infinite Graphs and its Applications" 2484 * by Kelly et al. 2485 * 2486 * We first check whether we can can split of any basic map M_i and 2487 * compute 2488 * 2489 * (\cup_j M_j)^+ 2490 * 2491 * as 2492 * 2493 * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+ 2494 * 2495 * using a recursive call on the remaining map. 2496 * 2497 * If not, we simply call box_closure on the whole map. 2498 */ 2499static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map, 2500 isl_bool *exact) 2501{ 2502 int i, j; 2503 isl_bool exact_i; 2504 isl_map *app; 2505 2506 if (!map) 2507 return NULL; 2508 if (map->n == 1) 2509 return box_closure_with_check(map, exact); 2510 2511 for (i = 0; i < map->n; ++i) { 2512 int ok; 2513 isl_map *qc, *tc; 2514 ok = can_be_split_off(map, i, &tc, &qc); 2515 if (ok < 0) 2516 goto error; 2517 if (!ok) 2518 continue; 2519 2520 app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0); 2521 2522 for (j = 0; j < map->n; ++j) { 2523 if (j == i) 2524 continue; 2525 app = isl_map_add_basic_map(app, 2526 isl_basic_map_copy(map->p[j])); 2527 } 2528 2529 app = isl_map_apply_range(isl_map_copy(qc), app); 2530 app = isl_map_apply_range(app, qc); 2531 2532 app = isl_map_union(tc, transitive_closure_omega(app, NULL)); 2533 exact_i = check_exactness_omega(map, app); 2534 if (exact_i == isl_bool_true) { 2535 if (exact) 2536 *exact = exact_i; 2537 isl_map_free(map); 2538 return app; 2539 } 2540 isl_map_free(app); 2541 if (exact_i < 0) 2542 goto error; 2543 } 2544 2545 return box_closure_with_check(map, exact); 2546error: 2547 isl_map_free(map); 2548 return NULL; 2549} 2550 2551/* Compute the transitive closure of "map", or an overapproximation. 2552 * If the result is exact, then *exact is set to 1. 2553 * Simply use map_power to compute the powers of map, but tell 2554 * it to project out the lengths of the paths instead of equating 2555 * the length to a parameter. 2556 */ 2557__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map, 2558 isl_bool *exact) 2559{ 2560 isl_space *target_dim; 2561 isl_bool closed; 2562 2563 if (!map) 2564 goto error; 2565 2566 if (map->ctx->opt->closure == ISL_CLOSURE_BOX) 2567 return transitive_closure_omega(map, exact); 2568 2569 map = isl_map_compute_divs(map); 2570 map = isl_map_coalesce(map); 2571 closed = isl_map_is_transitively_closed(map); 2572 if (closed < 0) 2573 goto error; 2574 if (closed) { 2575 if (exact) 2576 *exact = isl_bool_true; 2577 return map; 2578 } 2579 2580 target_dim = isl_map_get_space(map); 2581 map = map_power(map, exact, 1); 2582 map = isl_map_reset_space(map, target_dim); 2583 2584 return map; 2585error: 2586 isl_map_free(map); 2587 return NULL; 2588} 2589 2590static isl_stat inc_count(__isl_take isl_map *map, void *user) 2591{ 2592 int *n = user; 2593 2594 *n += map->n; 2595 2596 isl_map_free(map); 2597 2598 return isl_stat_ok; 2599} 2600 2601static isl_stat collect_basic_map(__isl_take isl_map *map, void *user) 2602{ 2603 int i; 2604 isl_basic_map ***next = user; 2605 2606 for (i = 0; i < map->n; ++i) { 2607 **next = isl_basic_map_copy(map->p[i]); 2608 if (!**next) 2609 goto error; 2610 (*next)++; 2611 } 2612 2613 isl_map_free(map); 2614 return isl_stat_ok; 2615error: 2616 isl_map_free(map); 2617 return isl_stat_error; 2618} 2619 2620/* Perform Floyd-Warshall on the given list of basic relations. 2621 * The basic relations may live in different dimensions, 2622 * but basic relations that get assigned to the diagonal of the 2623 * grid have domains and ranges of the same dimension and so 2624 * the standard algorithm can be used because the nested transitive 2625 * closures are only applied to diagonal elements and because all 2626 * compositions are performed on relations with compatible domains and ranges. 2627 */ 2628static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx, 2629 __isl_keep isl_basic_map **list, int n, isl_bool *exact) 2630{ 2631 int i, j, k; 2632 int n_group; 2633 int *group = NULL; 2634 isl_set **set = NULL; 2635 isl_map ***grid = NULL; 2636 isl_union_map *app; 2637 2638 group = setup_groups(ctx, list, n, &set, &n_group); 2639 if (!group) 2640 goto error; 2641 2642 grid = isl_calloc_array(ctx, isl_map **, n_group); 2643 if (!grid) 2644 goto error; 2645 for (i = 0; i < n_group; ++i) { 2646 grid[i] = isl_calloc_array(ctx, isl_map *, n_group); 2647 if (!grid[i]) 2648 goto error; 2649 for (j = 0; j < n_group; ++j) { 2650 isl_space *space1, *space2, *space; 2651 space1 = isl_space_reverse(isl_set_get_space(set[i])); 2652 space2 = isl_set_get_space(set[j]); 2653 space = isl_space_join(space1, space2); 2654 grid[i][j] = isl_map_empty(space); 2655 } 2656 } 2657 2658 for (k = 0; k < n; ++k) { 2659 i = group[2 * k]; 2660 j = group[2 * k + 1]; 2661 grid[i][j] = isl_map_union(grid[i][j], 2662 isl_map_from_basic_map( 2663 isl_basic_map_copy(list[k]))); 2664 } 2665 2666 floyd_warshall_iterate(grid, n_group, exact); 2667 2668 app = isl_union_map_empty(isl_map_get_space(grid[0][0])); 2669 2670 for (i = 0; i < n_group; ++i) { 2671 for (j = 0; j < n_group; ++j) 2672 app = isl_union_map_add_map(app, grid[i][j]); 2673 free(grid[i]); 2674 } 2675 free(grid); 2676 2677 for (i = 0; i < 2 * n; ++i) 2678 isl_set_free(set[i]); 2679 free(set); 2680 2681 free(group); 2682 return app; 2683error: 2684 if (grid) 2685 for (i = 0; i < n_group; ++i) { 2686 if (!grid[i]) 2687 continue; 2688 for (j = 0; j < n_group; ++j) 2689 isl_map_free(grid[i][j]); 2690 free(grid[i]); 2691 } 2692 free(grid); 2693 if (set) { 2694 for (i = 0; i < 2 * n; ++i) 2695 isl_set_free(set[i]); 2696 free(set); 2697 } 2698 free(group); 2699 return NULL; 2700} 2701 2702/* Perform Floyd-Warshall on the given union relation. 2703 * The implementation is very similar to that for non-unions. 2704 * The main difference is that it is applied unconditionally. 2705 * We first extract a list of basic maps from the union map 2706 * and then perform the algorithm on this list. 2707 */ 2708static __isl_give isl_union_map *union_floyd_warshall( 2709 __isl_take isl_union_map *umap, isl_bool *exact) 2710{ 2711 int i, n; 2712 isl_ctx *ctx; 2713 isl_basic_map **list = NULL; 2714 isl_basic_map **next; 2715 isl_union_map *res; 2716 2717 n = 0; 2718 if (isl_union_map_foreach_map(umap, inc_count, &n) < 0) 2719 goto error; 2720 2721 ctx = isl_union_map_get_ctx(umap); 2722 list = isl_calloc_array(ctx, isl_basic_map *, n); 2723 if (!list) 2724 goto error; 2725 2726 next = list; 2727 if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0) 2728 goto error; 2729 2730 res = union_floyd_warshall_on_list(ctx, list, n, exact); 2731 2732 if (list) { 2733 for (i = 0; i < n; ++i) 2734 isl_basic_map_free(list[i]); 2735 free(list); 2736 } 2737 2738 isl_union_map_free(umap); 2739 return res; 2740error: 2741 if (list) { 2742 for (i = 0; i < n; ++i) 2743 isl_basic_map_free(list[i]); 2744 free(list); 2745 } 2746 isl_union_map_free(umap); 2747 return NULL; 2748} 2749 2750/* Decompose the give union relation into strongly connected components. 2751 * The implementation is essentially the same as that of 2752 * construct_power_components with the major difference that all 2753 * operations are performed on union maps. 2754 */ 2755static __isl_give isl_union_map *union_components( 2756 __isl_take isl_union_map *umap, isl_bool *exact) 2757{ 2758 int i; 2759 int n; 2760 isl_ctx *ctx; 2761 isl_basic_map **list = NULL; 2762 isl_basic_map **next; 2763 isl_union_map *path = NULL; 2764 struct isl_tc_follows_data data; 2765 struct isl_tarjan_graph *g = NULL; 2766 int c, l; 2767 int recheck = 0; 2768 2769 n = 0; 2770 if (isl_union_map_foreach_map(umap, inc_count, &n) < 0) 2771 goto error; 2772 2773 if (n == 0) 2774 return umap; 2775 if (n <= 1) 2776 return union_floyd_warshall(umap, exact); 2777 2778 ctx = isl_union_map_get_ctx(umap); 2779 list = isl_calloc_array(ctx, isl_basic_map *, n); 2780 if (!list) 2781 goto error; 2782 2783 next = list; 2784 if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0) 2785 goto error; 2786 2787 data.list = list; 2788 data.check_closed = 0; 2789 g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data); 2790 if (!g) 2791 goto error; 2792 2793 c = 0; 2794 i = 0; 2795 l = n; 2796 path = isl_union_map_empty(isl_union_map_get_space(umap)); 2797 while (l) { 2798 isl_union_map *comp; 2799 isl_union_map *path_comp, *path_comb; 2800 comp = isl_union_map_empty(isl_union_map_get_space(umap)); 2801 while (g->order[i] != -1) { 2802 comp = isl_union_map_add_map(comp, 2803 isl_map_from_basic_map( 2804 isl_basic_map_copy(list[g->order[i]]))); 2805 --l; 2806 ++i; 2807 } 2808 path_comp = union_floyd_warshall(comp, exact); 2809 path_comb = isl_union_map_apply_range(isl_union_map_copy(path), 2810 isl_union_map_copy(path_comp)); 2811 path = isl_union_map_union(path, path_comp); 2812 path = isl_union_map_union(path, path_comb); 2813 ++i; 2814 ++c; 2815 } 2816 2817 if (c > 1 && data.check_closed && !*exact) { 2818 isl_bool closed; 2819 2820 closed = isl_union_map_is_transitively_closed(path); 2821 if (closed < 0) 2822 goto error; 2823 recheck = !closed; 2824 } 2825 2826 isl_tarjan_graph_free(g); 2827 2828 for (i = 0; i < n; ++i) 2829 isl_basic_map_free(list[i]); 2830 free(list); 2831 2832 if (recheck) { 2833 isl_union_map_free(path); 2834 return union_floyd_warshall(umap, exact); 2835 } 2836 2837 isl_union_map_free(umap); 2838 2839 return path; 2840error: 2841 isl_tarjan_graph_free(g); 2842 if (list) { 2843 for (i = 0; i < n; ++i) 2844 isl_basic_map_free(list[i]); 2845 free(list); 2846 } 2847 isl_union_map_free(umap); 2848 isl_union_map_free(path); 2849 return NULL; 2850} 2851 2852/* Compute the transitive closure of "umap", or an overapproximation. 2853 * If the result is exact, then *exact is set to 1. 2854 */ 2855__isl_give isl_union_map *isl_union_map_transitive_closure( 2856 __isl_take isl_union_map *umap, isl_bool *exact) 2857{ 2858 isl_bool closed; 2859 2860 if (!umap) 2861 return NULL; 2862 2863 if (exact) 2864 *exact = isl_bool_true; 2865 2866 umap = isl_union_map_compute_divs(umap); 2867 umap = isl_union_map_coalesce(umap); 2868 closed = isl_union_map_is_transitively_closed(umap); 2869 if (closed < 0) 2870 goto error; 2871 if (closed) 2872 return umap; 2873 umap = union_components(umap, exact); 2874 return umap; 2875error: 2876 isl_union_map_free(umap); 2877 return NULL; 2878} 2879 2880struct isl_union_power { 2881 isl_union_map *pow; 2882 isl_bool *exact; 2883}; 2884 2885static isl_stat power(__isl_take isl_map *map, void *user) 2886{ 2887 struct isl_union_power *up = user; 2888 2889 map = isl_map_power(map, up->exact); 2890 up->pow = isl_union_map_from_map(map); 2891 2892 return isl_stat_error; 2893} 2894 2895/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "space". 2896 */ 2897static __isl_give isl_union_map *deltas_map(__isl_take isl_space *space) 2898{ 2899 isl_basic_map *bmap; 2900 2901 space = isl_space_add_dims(space, isl_dim_in, 1); 2902 space = isl_space_add_dims(space, isl_dim_out, 1); 2903 bmap = isl_basic_map_universe(space); 2904 bmap = isl_basic_map_deltas_map(bmap); 2905 2906 return isl_union_map_from_map(isl_map_from_basic_map(bmap)); 2907} 2908 2909/* Compute the positive powers of "map", or an overapproximation. 2910 * The result maps the exponent to a nested copy of the corresponding power. 2911 * If the result is exact, then *exact is set to 1. 2912 */ 2913__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap, 2914 isl_bool *exact) 2915{ 2916 isl_size n; 2917 isl_union_map *inc; 2918 isl_union_map *dm; 2919 2920 n = isl_union_map_n_map(umap); 2921 if (n < 0) 2922 return isl_union_map_free(umap); 2923 if (n == 0) 2924 return umap; 2925 if (n == 1) { 2926 struct isl_union_power up = { NULL, exact }; 2927 isl_union_map_foreach_map(umap, &power, &up); 2928 isl_union_map_free(umap); 2929 return up.pow; 2930 } 2931 inc = isl_union_map_from_map(increment(isl_union_map_get_space(umap))); 2932 umap = isl_union_map_product(inc, umap); 2933 umap = isl_union_map_transitive_closure(umap, exact); 2934 umap = isl_union_map_zip(umap); 2935 dm = deltas_map(isl_union_map_get_space(umap)); 2936 umap = isl_union_map_apply_domain(umap, dm); 2937 2938 return umap; 2939} 2940 2941#undef TYPE 2942#define TYPE isl_map 2943#include "isl_power_templ.c" 2944 2945#undef TYPE 2946#define TYPE isl_union_map 2947#include "isl_power_templ.c" 2948