2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
12 #include "isl_map_private.h"
15 * The transitive closure implementation is based on the paper
16 * "Computing the Transitive Closure of a Union of Affine Integer
17 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
21 /* Given a union of translations (uniform dependences), return a matrix
22 * with as many rows as there are disjuncts in the union and as many
23 * columns as there are input dimensions (which should be equal to
24 * the number of output dimensions).
25 * Each row contains the translation performed by the corresponding disjunct.
26 * If "map" turns out not to be a union of translations, then the contents
27 * of the returned matrix are undefined and *ok is set to 0.
29 static __isl_give isl_mat *extract_steps(__isl_keep isl_map *map, int *ok)
32 struct isl_mat *steps;
33 unsigned dim = isl_map_dim(map, isl_dim_in);
37 steps = isl_mat_alloc(map->ctx, map->n, dim);
41 for (i = 0; i < map->n; ++i) {
42 struct isl_basic_set *delta;
44 delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
46 for (j = 0; j < dim; ++j) {
49 fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
52 isl_basic_set_free(delta);
59 isl_basic_set_free(delta);
74 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
75 * of the given dimension specification that maps a element x to any
76 * element that can be reached by taking a positive number of steps
77 * along any of the offsets, where the number of steps k is equal to
78 * parameter "param". That is, construct
80 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v_i, k = \sum_i k_i > 0 }
83 static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
84 __isl_keep isl_mat *steps, unsigned param)
87 struct isl_basic_map *path = NULL;
95 d = isl_dim_size(dim, isl_dim_in);
97 nparam = isl_dim_size(dim, isl_dim_param);
99 path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d + 1, n + 1);
101 for (i = 0; i < n; ++i) {
102 k = isl_basic_map_alloc_div(path);
105 isl_assert(steps->ctx, i == k, goto error);
106 isl_int_set_si(path->div[k][0], 0);
109 for (i = 0; i < d; ++i) {
110 k = isl_basic_map_alloc_equality(path);
113 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
114 isl_int_set_si(path->eq[k][1 + nparam + i], 1);
115 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
116 for (j = 0; j < n; ++j)
117 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
121 k = isl_basic_map_alloc_equality(path);
124 isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
125 isl_int_set_si(path->eq[k][1 + param], -1);
126 for (j = 0; j < n; ++j)
127 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
129 for (i = 0; i < n; ++i) {
130 k = isl_basic_map_alloc_inequality(path);
133 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
134 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
137 k = isl_basic_map_alloc_inequality(path);
140 isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
141 isl_int_set_si(path->ineq[k][1 + param], 1);
142 isl_int_set_si(path->ineq[k][0], -1);
146 path = isl_basic_map_simplify(path);
147 path = isl_basic_map_finalize(path);
148 return isl_map_from_basic_map(path);
151 isl_basic_map_free(path);
155 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
156 * construct a map that is an overapproximation of the map
157 * that takes an element from the space D to another
158 * element from the same space, such that the difference between
159 * them is a strictly positive sum of differences between images
160 * and pre-images in one of the R_i.
161 * The number of differences in the sum is equated to parameter "param".
164 * \Delta_i = { y - x | (x, y) in R_i }
166 * then the constructed map is an overapproximation of
168 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
169 * d = \sum_i k_i and k = \sum_i k_i > 0 }
171 * If any of the \Delta_i contains more than one element, then
172 * we currently simply return a universal map { (x) -> (y) }.
174 static __isl_give isl_map *construct_path(__isl_keep isl_map *map,
177 struct isl_mat *steps = NULL;
178 struct isl_map *path = NULL;
184 steps = extract_steps(map, &ok);
189 return isl_map_universe(isl_map_get_dim(map));
192 path = path_along_steps(isl_map_get_dim(map), steps, param);
198 /* Check whether "path" is acyclic.
199 * That is, check whether
201 * { d | d = y - x and (x,y) in path }
203 * does not contain the origin.
205 static int is_acyclic(__isl_take isl_map *path)
210 struct isl_set *delta;
212 delta = isl_map_deltas(path);
213 dim = isl_set_dim(delta, isl_dim_set);
214 for (i = 0; i < dim; ++i)
215 delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
217 acyclic = isl_set_is_empty(delta);
223 /* Shift variable at position "pos" up by one.
224 * That is, replace the corresponding variable v by v - 1.
226 static __isl_give isl_basic_map *basic_map_shift_pos(
227 __isl_take isl_basic_map *bmap, unsigned pos)
231 bmap = isl_basic_map_cow(bmap);
235 for (i = 0; i < bmap->n_eq; ++i)
236 isl_int_sub(bmap->eq[i][0], bmap->eq[i][0], bmap->eq[i][pos]);
238 for (i = 0; i < bmap->n_ineq; ++i)
239 isl_int_sub(bmap->ineq[i][0],
240 bmap->ineq[i][0], bmap->ineq[i][pos]);
242 for (i = 0; i < bmap->n_div; ++i) {
243 if (isl_int_is_zero(bmap->div[i][0]))
245 isl_int_sub(bmap->div[i][1],
246 bmap->div[i][1], bmap->div[i][1 + pos]);
252 /* Shift variable at position "pos" up by one.
253 * That is, replace the corresponding variable v by v - 1.
255 static __isl_give isl_map *map_shift_pos(__isl_take isl_map *map, unsigned pos)
259 map = isl_map_cow(map);
263 for (i = 0; i < map->n; ++i) {
264 map->p[i] = basic_map_shift_pos(map->p[i], pos);
268 ISL_F_CLR(map, ISL_MAP_NORMALIZED);
275 /* Check whether the overapproximation of the power of "map" is exactly
276 * the power of "map". Let R be "map" and A_k the overapproximation.
277 * The approximation is exact if
280 * A_k = A_{k-1} \circ R k >= 2
282 * Since A_k is known to be an overapproximation, we only need to check
285 * A_k \subset A_{k-1} \circ R k >= 2
288 static int check_power_exactness(__isl_take isl_map *map,
289 __isl_take isl_map *app, unsigned param)
295 app_1 = isl_map_fix_si(isl_map_copy(app), isl_dim_param, param, 1);
297 exact = isl_map_is_subset(app_1, map);
300 if (!exact || exact < 0) {
306 app_2 = isl_map_lower_bound_si(isl_map_copy(app),
307 isl_dim_param, param, 2);
308 app_1 = map_shift_pos(app, 1 + param);
309 app_1 = isl_map_apply_range(map, app_1);
311 exact = isl_map_is_subset(app_2, app_1);
319 /* Check whether the overapproximation of the power of "map" is exactly
320 * the power of "map", possibly after projecting out the power (if "project"
323 * If "project" is set and if "steps" can only result in acyclic paths,
326 * A = R \cup (A \circ R)
328 * where A is the overapproximation with the power projected out, i.e.,
329 * an overapproximation of the transitive closure.
330 * More specifically, since A is known to be an overapproximation, we check
332 * A \subset R \cup (A \circ R)
334 * Otherwise, we check if the power is exact.
336 static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
337 __isl_take isl_map *path, unsigned param, int project)
343 project = is_acyclic(path);
350 return check_power_exactness(map, app, param);
352 map = isl_map_project_out(map, isl_dim_param, param, 1);
353 app = isl_map_project_out(app, isl_dim_param, param, 1);
355 test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
356 test = isl_map_union(test, isl_map_copy(map));
358 exact = isl_map_is_subset(app, test);
372 /* Compute the positive powers of "map", or an overapproximation.
373 * The power is given by parameter "param". If the result is exact,
374 * then *exact is set to 1.
375 * If project is set, then we are actually interested in the transitive
376 * closure, so we can use a more relaxed exactness check.
378 static __isl_give isl_map *map_power(__isl_take isl_map *map, unsigned param,
379 int *exact, int project)
381 struct isl_set *domain = NULL;
382 struct isl_set *range = NULL;
383 struct isl_map *app = NULL;
384 struct isl_map *path = NULL;
389 map = isl_map_remove_empty_parts(map);
393 if (isl_map_fast_is_empty(map))
396 isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
398 isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
401 domain = isl_map_domain(isl_map_copy(map));
402 domain = isl_set_coalesce(domain);
403 range = isl_map_range(isl_map_copy(map));
404 range = isl_set_coalesce(range);
405 app = isl_map_from_domain_and_range(isl_set_copy(domain),
406 isl_set_copy(range));
408 path = construct_path(map, param);
409 app = isl_map_intersect(app, isl_map_copy(path));
412 (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
413 isl_map_copy(path), param, project)) < 0)
416 isl_set_free(domain);
422 isl_set_free(domain);
430 /* Compute the positive powers of "map", or an overapproximation.
431 * The power is given by parameter "param". If the result is exact,
432 * then *exact is set to 1.
434 __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
437 return map_power(map, param, exact, 0);
440 /* Compute the transitive closure of "map", or an overapproximation.
441 * If the result is exact, then *exact is set to 1.
442 * Simply compute the powers of map and then project out the parameter
443 * describing the power.
445 __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
453 param = isl_map_dim(map, isl_dim_param);
454 map = isl_map_add(map, isl_dim_param, 1);
455 map = map_power(map, param, exact, 1);
456 map = isl_map_project_out(map, isl_dim_param, param, 1);