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 /* Check whether any path of length at least one along "steps" is acyclic.
156 * That is, check whether
158 * \sum_i k_i \delta_i = 0
162 * with \delta_i the rows of "steps", is infeasible.
164 static int is_acyclic(__isl_keep isl_mat
*steps
)
168 struct isl_basic_set
*test
;
173 test
= isl_basic_set_alloc(steps
->ctx
, 0, steps
->n_row
, 0,
174 steps
->n_col
, steps
->n_row
+ 1);
176 for (i
= 0; i
< steps
->n_col
; ++i
) {
177 k
= isl_basic_set_alloc_equality(test
);
180 isl_int_set_si(test
->eq
[k
][0], 0);
181 for (j
= 0; j
< steps
->n_row
; ++j
)
182 isl_int_set(test
->eq
[k
][1 + j
], steps
->row
[j
][i
]);
184 for (j
= 0; j
< steps
->n_row
; ++j
) {
185 k
= isl_basic_set_alloc_inequality(test
);
188 isl_seq_clr(test
->ineq
[k
], 1 + steps
->n_row
);
189 isl_int_set_si(test
->ineq
[k
][1 + j
], 1);
192 k
= isl_basic_set_alloc_inequality(test
);
195 isl_int_set_si(test
->ineq
[k
][0], -1);
196 for (j
= 0; j
< steps
->n_row
; ++j
)
197 isl_int_set_si(test
->ineq
[k
][1 + j
], 1);
199 acyclic
= isl_basic_set_is_empty(test
);
201 isl_basic_set_free(test
);
204 isl_basic_set_free(test
);
208 /* Shift variable at position "pos" up by one.
209 * That is, replace the corresponding variable v by v - 1.
211 static __isl_give isl_basic_map
*basic_map_shift_pos(
212 __isl_take isl_basic_map
*bmap
, unsigned pos
)
216 bmap
= isl_basic_map_cow(bmap
);
220 for (i
= 0; i
< bmap
->n_eq
; ++i
)
221 isl_int_sub(bmap
->eq
[i
][0], bmap
->eq
[i
][0], bmap
->eq
[i
][pos
]);
223 for (i
= 0; i
< bmap
->n_ineq
; ++i
)
224 isl_int_sub(bmap
->ineq
[i
][0],
225 bmap
->ineq
[i
][0], bmap
->ineq
[i
][pos
]);
227 for (i
= 0; i
< bmap
->n_div
; ++i
) {
228 if (isl_int_is_zero(bmap
->div
[i
][0]))
230 isl_int_sub(bmap
->div
[i
][1],
231 bmap
->div
[i
][1], bmap
->div
[i
][1 + pos
]);
237 /* Shift variable at position "pos" up by one.
238 * That is, replace the corresponding variable v by v - 1.
240 static __isl_give isl_map
*map_shift_pos(__isl_take isl_map
*map
, unsigned pos
)
244 map
= isl_map_cow(map
);
248 for (i
= 0; i
< map
->n
; ++i
) {
249 map
->p
[i
] = basic_map_shift_pos(map
->p
[i
], pos
);
253 ISL_F_CLR(map
, ISL_MAP_NORMALIZED
);
260 /* Check whether the overapproximation of the power of "map" is exactly
261 * the power of "map". Let R be "map" and A_k the overapproximation.
262 * The approximation is exact if
265 * A_k = A_{k-1} \circ R k >= 2
267 * Since A_k is known to be an overapproximation, we only need to check
270 * A_k \subset A_{k-1} \circ R k >= 2
273 static int check_power_exactness(__isl_take isl_map
*map
,
274 __isl_take isl_map
*app
, unsigned param
)
280 app_1
= isl_map_fix_si(isl_map_copy(app
), isl_dim_param
, param
, 1);
282 exact
= isl_map_is_subset(app_1
, map
);
285 if (!exact
|| exact
< 0) {
291 app_2
= isl_map_lower_bound_si(isl_map_copy(app
),
292 isl_dim_param
, param
, 2);
293 app_1
= map_shift_pos(app
, 1 + param
);
294 app_1
= isl_map_apply_range(map
, app_1
);
296 exact
= isl_map_is_subset(app_2
, app_1
);
304 /* Check whether the overapproximation of the power of "map" is exactly
305 * the power of "map", possibly after projecting out the power (if "project"
308 * If "project" is set and if "steps" can only result in acyclic paths,
311 * A = R \cup (A \circ R)
313 * where A is the overapproximation with the power projected out, i.e.,
314 * an overapproximation of the transitive closure.
315 * More specifically, since A is known to be an overapproximation, we check
317 * A \subset R \cup (A \circ R)
319 * Otherwise, we check if the power is exact.
321 static int check_exactness(__isl_take isl_map
*map
, __isl_take isl_map
*app
,
322 __isl_keep isl_mat
*steps
, unsigned param
, int project
)
328 project
= is_acyclic(steps
);
334 return check_power_exactness(map
, app
, param
);
336 map
= isl_map_project_out(map
, isl_dim_param
, param
, 1);
337 app
= isl_map_project_out(app
, isl_dim_param
, param
, 1);
339 test
= isl_map_apply_range(isl_map_copy(map
), isl_map_copy(app
));
340 test
= isl_map_union(test
, isl_map_copy(map
));
342 exact
= isl_map_is_subset(app
, test
);
356 /* Compute the positive powers of "map", or an overapproximation.
357 * The power is given by parameter "param". If the result is exact,
358 * then *exact is set to 1.
359 * If project is set, then we are actually interested in the transitive
360 * closure, so we can use a more relaxed exactness check.
362 static __isl_give isl_map
*map_power(__isl_take isl_map
*map
, unsigned param
,
363 int *exact
, int project
)
365 struct isl_mat
*steps
= NULL
;
366 struct isl_set
*domain
= NULL
;
367 struct isl_set
*range
= NULL
;
368 struct isl_map
*app
= NULL
;
369 struct isl_map
*path
= NULL
;
375 map
= isl_map_remove_empty_parts(map
);
379 if (isl_map_fast_is_empty(map
))
382 isl_assert(map
->ctx
, param
< isl_map_dim(map
, isl_dim_param
), goto error
);
384 isl_map_dim(map
, isl_dim_in
) == isl_map_dim(map
, isl_dim_out
),
387 domain
= isl_map_domain(isl_map_copy(map
));
388 domain
= isl_set_coalesce(domain
);
389 range
= isl_map_range(isl_map_copy(map
));
390 range
= isl_set_coalesce(range
);
391 app
= isl_map_from_domain_and_range(isl_set_copy(domain
),
392 isl_set_copy(range
));
394 steps
= extract_steps(map
, &ok
);
398 path
= path_along_steps(isl_map_get_dim(map
), steps
, param
);
399 app
= isl_map_intersect(app
, isl_map_copy(path
));
402 (*exact
= check_exactness(isl_map_copy(map
), isl_map_copy(app
),
403 steps
, param
, project
)) < 0)
411 isl_set_free(domain
);
418 isl_set_free(domain
);
427 /* Compute the positive powers of "map", or an overapproximation.
428 * The power is given by parameter "param". If the result is exact,
429 * then *exact is set to 1.
431 __isl_give isl_map
*isl_map_power(__isl_take isl_map
*map
, unsigned param
,
434 return map_power(map
, param
, exact
, 0);
437 /* Compute the transitive closure of "map", or an overapproximation.
438 * If the result is exact, then *exact is set to 1.
439 * Simply compute the powers of map and then project out the parameter
440 * describing the power.
442 __isl_give isl_map
*isl_map_transitive_closure(__isl_take isl_map
*map
,
450 param
= isl_map_dim(map
, isl_dim_param
);
451 map
= isl_map_add(map
, isl_dim_param
, 1);
452 map
= map_power(map
, param
, exact
, 1);
453 map
= isl_map_project_out(map
, isl_dim_param
, param
, 1);