Bug Summary

File:polly/lib/External/isl/isl_transitive_closure.c
Warning:line 2817, column 37
Dereference of null pointer (loaded from variable 'exact')

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-pc-linux-gnu -analyze -disable-free -disable-llvm-verifier -discard-value-names -main-file-name isl_transitive_closure.c -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -analyzer-config-compatibility-mode=true -mrelocation-model pic -pic-level 2 -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -ffunction-sections -fdata-sections -fcoverage-compilation-dir=/build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/tools/polly/lib/External -resource-dir /usr/lib/llvm-14/lib/clang/14.0.0 -D _DEBUG -D _GNU_SOURCE -D __STDC_CONSTANT_MACROS -D __STDC_FORMAT_MACROS -D __STDC_LIMIT_MACROS -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/tools/polly/lib/External -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/isl -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/isl/include -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/isl/imath -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/tools/polly/lib/External/isl -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/tools/polly/include -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/pet/include -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/tools/polly/lib/External/isl/include -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/include -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/include -I /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/llvm/include -D NDEBUG -U NDEBUG -internal-isystem /usr/lib/llvm-14/lib/clang/14.0.0/include -internal-isystem /usr/local/include -internal-isystem /usr/lib/gcc/x86_64-linux-gnu/10/../../../../x86_64-linux-gnu/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -Wno-unused-parameter -Wwrite-strings -Wno-missing-field-initializers -Wno-long-long -Wno-comment -std=gnu99 -fconst-strings -fdebug-compilation-dir=/build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/build-llvm/tools/polly/lib/External -fdebug-prefix-map=/build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0=. -ferror-limit 19 -stack-protector 2 -fgnuc-version=4.2.1 -vectorize-loops -vectorize-slp -analyzer-output=html -analyzer-config stable-report-filename=true -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /tmp/scan-build-2021-08-28-193554-24367-1 -x c /build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/isl/isl_transitive_closure.c
1/*
2 * Copyright 2010 INRIA Saclay
3 *
4 * Use of this software is governed by the MIT license
5 *
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8 * 91893 Orsay, France
9 */
10
11#include <isl_ctx_private.h>
12#include <isl_map_private.h>
13#include <isl/map.h>
14#include <isl_seq.h>
15#include <isl_space_private.h>
16#include <isl_lp_private.h>
17#include <isl/union_map.h>
18#include <isl_mat_private.h>
19#include <isl_vec_private.h>
20#include <isl_options_private.h>
21#include <isl_tarjan.h>
22
23isl_bool isl_map_is_transitively_closed(__isl_keep isl_map *map)
24{
25 isl_map *map2;
26 isl_bool closed;
27
28 map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));
29 closed = isl_map_is_subset(map2, map);
30 isl_map_free(map2);
31
32 return closed;
33}
34
35isl_bool isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)
36{
37 isl_union_map *umap2;
38 isl_bool closed;
39
40 umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),
41 isl_union_map_copy(umap));
42 closed = isl_union_map_is_subset(umap2, umap);
43 isl_union_map_free(umap2);
44
45 return closed;
46}
47
48/* Given a map that represents a path with the length of the path
49 * encoded as the difference between the last output coordindate
50 * and the last input coordinate, set this length to either
51 * exactly "length" (if "exactly" is set) or at least "length"
52 * (if "exactly" is not set).
53 */
54static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
55 int exactly, int length)
56{
57 isl_space *space;
58 struct isl_basic_map *bmap;
59 isl_size d;
60 isl_size nparam;
61 isl_size total;
62 int k;
63 isl_int *c;
64
65 if (!map)
66 return NULL((void*)0);
67
68 space = isl_map_get_space(map);
69 d = isl_space_dim(space, isl_dim_in);
70 nparam = isl_space_dim(space, isl_dim_param);
71 total = isl_space_dim(space, isl_dim_all);
72 if (d < 0 || nparam < 0 || total < 0)
73 space = isl_space_free(space);
74 bmap = isl_basic_map_alloc_space(space, 0, 1, 1);
75 if (exactly) {
76 k = isl_basic_map_alloc_equality(bmap);
77 if (k < 0)
78 goto error;
79 c = bmap->eq[k];
80 } else {
81 k = isl_basic_map_alloc_inequality(bmap);
82 if (k < 0)
83 goto error;
84 c = bmap->ineq[k];
85 }
86 isl_seq_clr(c, 1 + total);
87 isl_int_set_si(c[0], -length)isl_sioimath_set_si((c[0]), -length);
88 isl_int_set_si(c[1 + nparam + d - 1], -1)isl_sioimath_set_si((c[1 + nparam + d - 1]), -1);
89 isl_int_set_si(c[1 + nparam + d + d - 1], 1)isl_sioimath_set_si((c[1 + nparam + d + d - 1]), 1);
90
91 bmap = isl_basic_map_finalize(bmap);
92 map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
93
94 return map;
95error:
96 isl_basic_map_free(bmap);
97 isl_map_free(map);
98 return NULL((void*)0);
99}
100
101/* Check whether the overapproximation of the power of "map" is exactly
102 * the power of "map". Let R be "map" and A_k the overapproximation.
103 * The approximation is exact if
104 *
105 * A_1 = R
106 * A_k = A_{k-1} \circ R k >= 2
107 *
108 * Since A_k is known to be an overapproximation, we only need to check
109 *
110 * A_1 \subset R
111 * A_k \subset A_{k-1} \circ R k >= 2
112 *
113 * In practice, "app" has an extra input and output coordinate
114 * to encode the length of the path. So, we first need to add
115 * this coordinate to "map" and set the length of the path to
116 * one.
117 */
118static isl_bool check_power_exactness(__isl_take isl_map *map,
119 __isl_take isl_map *app)
120{
121 isl_bool exact;
122 isl_map *app_1;
123 isl_map *app_2;
124
125 map = isl_map_add_dims(map, isl_dim_in, 1);
126 map = isl_map_add_dims(map, isl_dim_out, 1);
127 map = set_path_length(map, 1, 1);
128
129 app_1 = set_path_length(isl_map_copy(app), 1, 1);
130
131 exact = isl_map_is_subset(app_1, map);
132 isl_map_free(app_1);
133
134 if (!exact || exact < 0) {
135 isl_map_free(app);
136 isl_map_free(map);
137 return exact;
138 }
139
140 app_1 = set_path_length(isl_map_copy(app), 0, 1);
141 app_2 = set_path_length(app, 0, 2);
142 app_1 = isl_map_apply_range(map, app_1);
143
144 exact = isl_map_is_subset(app_2, app_1);
145
146 isl_map_free(app_1);
147 isl_map_free(app_2);
148
149 return exact;
150}
151
152/* Check whether the overapproximation of the power of "map" is exactly
153 * the power of "map", possibly after projecting out the power (if "project"
154 * is set).
155 *
156 * If "project" is set and if "steps" can only result in acyclic paths,
157 * then we check
158 *
159 * A = R \cup (A \circ R)
160 *
161 * where A is the overapproximation with the power projected out, i.e.,
162 * an overapproximation of the transitive closure.
163 * More specifically, since A is known to be an overapproximation, we check
164 *
165 * A \subset R \cup (A \circ R)
166 *
167 * Otherwise, we check if the power is exact.
168 *
169 * Note that "app" has an extra input and output coordinate to encode
170 * the length of the part. If we are only interested in the transitive
171 * closure, then we can simply project out these coordinates first.
172 */
173static isl_bool check_exactness(__isl_take isl_map *map,
174 __isl_take isl_map *app, int project)
175{
176 isl_map *test;
177 isl_bool exact;
178 isl_size d;
179
180 if (!project)
181 return check_power_exactness(map, app);
182
183 d = isl_map_dim(map, isl_dim_in);
184 if (d < 0)
185 app = isl_map_free(app);
186 app = set_path_length(app, 0, 1);
187 app = isl_map_project_out(app, isl_dim_in, d, 1);
188 app = isl_map_project_out(app, isl_dim_out, d, 1);
189
190 app = isl_map_reset_space(app, isl_map_get_space(map));
191
192 test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
193 test = isl_map_union(test, isl_map_copy(map));
194
195 exact = isl_map_is_subset(app, test);
196
197 isl_map_free(app);
198 isl_map_free(test);
199
200 isl_map_free(map);
201
202 return exact;
203}
204
205/*
206 * The transitive closure implementation is based on the paper
207 * "Computing the Transitive Closure of a Union of Affine Integer
208 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
209 * Albert Cohen.
210 */
211
212/* Given a set of n offsets v_i (the rows of "steps"), construct a relation
213 * of the given dimension specification (Z^{n+1} -> Z^{n+1})
214 * that maps an element x to any element that can be reached
215 * by taking a non-negative number of steps along any of
216 * the extended offsets v'_i = [v_i 1].
217 * That is, construct
218 *
219 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
220 *
221 * For any element in this relation, the number of steps taken
222 * is equal to the difference in the final coordinates.
223 */
224static __isl_give isl_map *path_along_steps(__isl_take isl_space *space,
225 __isl_keep isl_mat *steps)
226{
227 int i, j, k;
228 struct isl_basic_map *path = NULL((void*)0);
229 isl_size d;
230 unsigned n;
231 isl_size nparam;
232 isl_size total;
233
234 d = isl_space_dim(space, isl_dim_in);
235 nparam = isl_space_dim(space, isl_dim_param);
236 if (d < 0 || nparam < 0 || !steps)
237 goto error;
238
239 n = steps->n_row;
240
241 path = isl_basic_map_alloc_space(isl_space_copy(space), n, d, n);
242
243 for (i = 0; i < n; ++i) {
244 k = isl_basic_map_alloc_div(path);
245 if (k < 0)
246 goto error;
247 isl_assert(steps->ctx, i == k, goto error)do { if (i == k) break; do { isl_handle_error(steps->ctx, isl_error_unknown
, "Assertion \"" "i == k" "\" failed", "/build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/isl/isl_transitive_closure.c"
, 247); goto error; } while (0); } while (0)
;
248 isl_int_set_si(path->div[k][0], 0)isl_sioimath_set_si((path->div[k][0]), 0);
249 }
250
251 total = isl_basic_map_dim(path, isl_dim_all);
252 if (total < 0)
253 goto error;
254 for (i = 0; i < d; ++i) {
255 k = isl_basic_map_alloc_equality(path);
256 if (k < 0)
257 goto error;
258 isl_seq_clr(path->eq[k], 1 + total);
259 isl_int_set_si(path->eq[k][1 + nparam + i], 1)isl_sioimath_set_si((path->eq[k][1 + nparam + i]), 1);
260 isl_int_set_si(path->eq[k][1 + nparam + d + i], -1)isl_sioimath_set_si((path->eq[k][1 + nparam + d + i]), -1);
261 if (i == d - 1)
262 for (j = 0; j < n; ++j)
263 isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1)isl_sioimath_set_si((path->eq[k][1 + nparam + 2 * d + j]),
1)
;
264 else
265 for (j = 0; j < n; ++j)
266 isl_int_set(path->eq[k][1 + nparam + 2 * d + j],isl_sioimath_set((path->eq[k][1 + nparam + 2 * d + j]), *(
steps->row[j][i]))
267 steps->row[j][i])isl_sioimath_set((path->eq[k][1 + nparam + 2 * d + j]), *(
steps->row[j][i]))
;
268 }
269
270 for (i = 0; i < n; ++i) {
271 k = isl_basic_map_alloc_inequality(path);
272 if (k < 0)
273 goto error;
274 isl_seq_clr(path->ineq[k], 1 + total);
275 isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1)isl_sioimath_set_si((path->ineq[k][1 + nparam + 2 * d + i]
), 1)
;
276 }
277
278 isl_space_free(space);
279
280 path = isl_basic_map_simplify(path);
281 path = isl_basic_map_finalize(path);
282 return isl_map_from_basic_map(path);
283error:
284 isl_space_free(space);
285 isl_basic_map_free(path);
286 return NULL((void*)0);
287}
288
289#define IMPURE0 0
290#define PURE_PARAM1 1
291#define PURE_VAR2 2
292#define MIXED3 3
293
294/* Check whether the parametric constant term of constraint c is never
295 * positive in "bset".
296 */
297static isl_bool parametric_constant_never_positive(
298 __isl_keep isl_basic_setisl_basic_map *bset, isl_int *c, int *div_purity)
299{
300 isl_size d;
301 isl_size n_div;
302 isl_size nparam;
303 isl_size total;
304 int i;
305 int k;
306 isl_bool empty;
307
308 n_div = isl_basic_set_dim(bset, isl_dim_div);
309 d = isl_basic_set_dim(bset, isl_dim_set);
310 nparam = isl_basic_set_dim(bset, isl_dim_param);
311 total = isl_basic_set_dim(bset, isl_dim_all);
312 if (n_div < 0 || d < 0 || nparam < 0 || total < 0)
313 return isl_bool_error;
314
315 bset = isl_basic_set_copy(bset);
316 bset = isl_basic_set_cow(bset);
317 bset = isl_basic_set_extend_constraints(bset, 0, 1);
318 k = isl_basic_set_alloc_inequality(bset);
319 if (k < 0)
320 goto error;
321 isl_seq_clr(bset->ineq[k], 1 + total);
322 isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
323 for (i = 0; i < n_div; ++i) {
324 if (div_purity[i] != PURE_PARAM1)
325 continue;
326 isl_int_set(bset->ineq[k][1 + nparam + d + i],isl_sioimath_set((bset->ineq[k][1 + nparam + d + i]), *(c[
1 + nparam + d + i]))
327 c[1 + nparam + d + i])isl_sioimath_set((bset->ineq[k][1 + nparam + d + i]), *(c[
1 + nparam + d + i]))
;
328 }
329 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1)isl_sioimath_sub_ui((bset->ineq[k][0]), *(bset->ineq[k]
[0]), 1)
;
330 empty = isl_basic_set_is_empty(bset);
331 isl_basic_set_free(bset);
332
333 return empty;
334error:
335 isl_basic_set_free(bset);
336 return isl_bool_error;
337}
338
339/* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
340 * Return PURE_VAR if only the coefficients of the set variables are non-zero.
341 * Return MIXED if only the coefficients of the parameters and the set
342 * variables are non-zero and if moreover the parametric constant
343 * can never attain positive values.
344 * Return IMPURE otherwise.
345 */
346static int purity(__isl_keep isl_basic_setisl_basic_map *bset, isl_int *c, int *div_purity,
347 int eq)
348{
349 isl_size d;
350 isl_size n_div;
351 isl_size nparam;
352 isl_bool empty;
353 int i;
354 int p = 0, v = 0;
355
356 n_div = isl_basic_set_dim(bset, isl_dim_div);
357 d = isl_basic_set_dim(bset, isl_dim_set);
358 nparam = isl_basic_set_dim(bset, isl_dim_param);
359 if (n_div < 0 || d < 0 || nparam < 0)
360 return -1;
361
362 for (i = 0; i < n_div; ++i) {
363 if (isl_int_is_zero(c[1 + nparam + d + i])(isl_sioimath_sgn(*(c[1 + nparam + d + i])) == 0))
364 continue;
365 switch (div_purity[i]) {
366 case PURE_PARAM1: p = 1; break;
367 case PURE_VAR2: v = 1; break;
368 default: return IMPURE0;
369 }
370 }
371 if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1)
372 return PURE_VAR2;
373 if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
374 return PURE_PARAM1;
375
376 empty = parametric_constant_never_positive(bset, c, div_purity);
377 if (eq && empty >= 0 && !empty) {
378 isl_seq_neg(c, c, 1 + nparam + d + n_div);
379 empty = parametric_constant_never_positive(bset, c, div_purity);
380 }
381
382 return empty < 0 ? -1 : empty ? MIXED3 : IMPURE0;
383}
384
385/* Return an array of integers indicating the type of each div in bset.
386 * If the div is (recursively) defined in terms of only the parameters,
387 * then the type is PURE_PARAM.
388 * If the div is (recursively) defined in terms of only the set variables,
389 * then the type is PURE_VAR.
390 * Otherwise, the type is IMPURE.
391 */
392static __isl_give int *get_div_purity(__isl_keep isl_basic_setisl_basic_map *bset)
393{
394 int i, j;
395 int *div_purity;
396 isl_size d;
397 isl_size n_div;
398 isl_size nparam;
399
400 n_div = isl_basic_set_dim(bset, isl_dim_div);
401 d = isl_basic_set_dim(bset, isl_dim_set);
402 nparam = isl_basic_set_dim(bset, isl_dim_param);
403 if (n_div < 0 || d < 0 || nparam < 0)
404 return NULL((void*)0);
405
406 div_purity = isl_alloc_array(bset->ctx, int, n_div)((int *)isl_malloc_or_die(bset->ctx, (n_div)*sizeof(int)));
407 if (n_div && !div_purity)
408 return NULL((void*)0);
409
410 for (i = 0; i < bset->n_div; ++i) {
411 int p = 0, v = 0;
412 if (isl_int_is_zero(bset->div[i][0])(isl_sioimath_sgn(*(bset->div[i][0])) == 0)) {
413 div_purity[i] = IMPURE0;
414 continue;
415 }
416 if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1)
417 p = 1;
418 if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)
419 v = 1;
420 for (j = 0; j < i; ++j) {
421 if (isl_int_is_zero(bset->div[i][2 + nparam + d + j])(isl_sioimath_sgn(*(bset->div[i][2 + nparam + d + j])) == 0
)
)
422 continue;
423 switch (div_purity[j]) {
424 case PURE_PARAM1: p = 1; break;
425 case PURE_VAR2: v = 1; break;
426 default: p = v = 1; break;
427 }
428 }
429 div_purity[i] = v ? p ? IMPURE0 : PURE_VAR2 : PURE_PARAM1;
430 }
431
432 return div_purity;
433}
434
435/* Given a path with the as yet unconstrained length at div position "pos",
436 * check if setting the length to zero results in only the identity
437 * mapping.
438 */
439static isl_bool empty_path_is_identity(__isl_keep isl_basic_map *path,
440 unsigned pos)
441{
442 isl_basic_map *test = NULL((void*)0);
443 isl_basic_map *id = NULL((void*)0);
444 isl_bool is_id;
445
446 test = isl_basic_map_copy(path);
447 test = isl_basic_map_fix_si(test, isl_dim_div, pos, 0);
448 id = isl_basic_map_identity(isl_basic_map_get_space(path));
449 is_id = isl_basic_map_is_equal(test, id);
450 isl_basic_map_free(test);
451 isl_basic_map_free(id);
452 return is_id;
453}
454
455/* If any of the constraints is found to be impure then this function
456 * sets *impurity to 1.
457 *
458 * If impurity is NULL then we are dealing with a non-parametric set
459 * and so the constraints are obviously PURE_VAR.
460 */
461static __isl_give isl_basic_map *add_delta_constraints(
462 __isl_take isl_basic_map *path,
463 __isl_keep isl_basic_setisl_basic_map *delta, unsigned off, unsigned nparam,
464 unsigned d, int *div_purity, int eq, int *impurity)
465{
466 int i, k;
467 int n = eq ? delta->n_eq : delta->n_ineq;
468 isl_int **delta_c = eq ? delta->eq : delta->ineq;
469 isl_size n_div, total;
470
471 n_div = isl_basic_set_dim(delta, isl_dim_div);
472 total = isl_basic_map_dim(path, isl_dim_all);
473 if (n_div < 0 || total < 0)
474 return isl_basic_map_free(path);
475
476 for (i = 0; i < n; ++i) {
477 isl_int *path_c;
478 int p = PURE_VAR2;
479 if (impurity)
480 p = purity(delta, delta_c[i], div_purity, eq);
481 if (p < 0)
482 goto error;
483 if (p != PURE_VAR2 && p != PURE_PARAM1 && !*impurity)
484 *impurity = 1;
485 if (p == IMPURE0)
486 continue;
487 if (eq && p != MIXED3) {
488 k = isl_basic_map_alloc_equality(path);
489 if (k < 0)
490 goto error;
491 path_c = path->eq[k];
492 } else {
493 k = isl_basic_map_alloc_inequality(path);
494 if (k < 0)
495 goto error;
496 path_c = path->ineq[k];
497 }
498 isl_seq_clr(path_c, 1 + total);
499 if (p == PURE_VAR2) {
500 isl_seq_cpy(path_c + off,
501 delta_c[i] + 1 + nparam, d);
502 isl_int_set(path_c[off + d], delta_c[i][0])isl_sioimath_set((path_c[off + d]), *(delta_c[i][0]));
503 } else if (p == PURE_PARAM1) {
504 isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
505 } else {
506 isl_seq_cpy(path_c + off,
507 delta_c[i] + 1 + nparam, d);
508 isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
509 }
510 isl_seq_cpy(path_c + off - n_div,
511 delta_c[i] + 1 + nparam + d, n_div);
512 }
513
514 return path;
515error:
516 isl_basic_map_free(path);
517 return NULL((void*)0);
518}
519
520/* Given a set of offsets "delta", construct a relation of the
521 * given dimension specification (Z^{n+1} -> Z^{n+1}) that
522 * is an overapproximation of the relations that
523 * maps an element x to any element that can be reached
524 * by taking a non-negative number of steps along any of
525 * the elements in "delta".
526 * That is, construct an approximation of
527 *
528 * { [x] -> [y] : exists f \in \delta, k \in Z :
529 * y = x + k [f, 1] and k >= 0 }
530 *
531 * For any element in this relation, the number of steps taken
532 * is equal to the difference in the final coordinates.
533 *
534 * In particular, let delta be defined as
535 *
536 * \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and
537 * C x + C'p + c >= 0 and
538 * D x + D'p + d >= 0 }
539 *
540 * where the constraints C x + C'p + c >= 0 are such that the parametric
541 * constant term of each constraint j, "C_j x + C'_j p + c_j",
542 * can never attain positive values, then the relation is constructed as
543 *
544 * { [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_setisl_basic_map *delta)
577{
578 isl_basic_map *path = NULL((void*)0);
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((void*)0);
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)isl_sioimath_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)isl_sioimath_set_si((path->eq[k][1 + nparam + i]), 1);
614 isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1)isl_sioimath_set_si((path->eq[k][1 + nparam + d + 1 + i]),
-1)
;
615 isl_int_set_si(path->eq[k][off + i], 1)isl_sioimath_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((void*)0), 1, NULL((void*)0));
638 path = add_delta_constraints(path, delta, off, nparam, d,
639 NULL((void*)0), 0, NULL((void*)0));
640 path = isl_basic_map_gauss(path, NULL((void*)0));
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)isl_sioimath_set_si((path->ineq[k][0]), -1);
653 isl_int_set_si(path->ineq[k][off + d], 1)isl_sioimath_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((void*)0);
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)isl_sioimath_set_si((bmap->eq[k][1 + param]), -1);
698 isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1)isl_sioimath_set_si((bmap->eq[k][1 + nparam + d - 1]), -1);
699 isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1)isl_sioimath_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)isl_sioimath_set_si((bmap->ineq[k][1 + param]), 1);
706 isl_int_set_si(bmap->ineq[k][0], -1)isl_sioimath_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((void*)0);
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_setisl_map *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((void*)0);
777 struct isl_map *path = NULL((void*)0);
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_setisl_basic_map *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((void*)0);
841}
842
843static isl_bool isl_set_overlaps(__isl_keep isl_setisl_map *set1,
844 __isl_keep isl_setisl_map *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_setisl_map *domain = NULL((void*)0);
872 struct isl_setisl_map *range = NULL((void*)0);
873 struct isl_map *app = NULL((void*)0);
874 struct isl_map *path = NULL((void*)0);
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((void*)0);
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((void*)0));
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((void*)0);
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((void*)0);
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_setisl_map *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((void*)0);
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_setisl_map *dom, __isl_keep isl_setisl_map *ran)
975{
976 isl_setisl_map *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 LEFT2 2
1010#define RIGHT1 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_setisl_map *C, int i,
1042 isl_setisl_map **dom, isl_setisl_map **ran, int *left, int *right,
1043 __isl_keep isl_map *map)
1044{
1045 int j;
1046 int ok;
1047
1048 ok = LEFT2 | RIGHT1;
1049 for (j = 0; j < map->n && ok; ++j) {
1050 isl_bool overlaps, subset;
1051 if (j == i)
1052 continue;
1053
1054 if (ok & RIGHT1) {
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 &= ~RIGHT1;
1074 }
1075 }
1076
1077 if (ok & LEFT2) {
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 &= ~LEFT2;
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((void*)0);
1167
1168 if (!map)
1169 goto error;
1170 isl_assert(map->ctx, left || right, goto error)do { if (left || right) break; do { isl_handle_error(map->
ctx, isl_error_unknown, "Assertion \"" "left || right" "\" failed"
, "/build/llvm-toolchain-snapshot-14~++20210828111110+16086d47c0d0/polly/lib/External/isl/isl_transitive_closure.c"
, 1170); goto error; } while (0); } while (0)
;
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((void*)0);
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_setisl_map **dom, isl_setisl_map **ran, int *left, int *right,
1223 __isl_give isl_map **res)
1224{
1225 int i;
1226 isl_setisl_map *C;
1227 isl_size d;
1228
1229 *res = NULL((void*)0);
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((void*)0);
1289 }
1290
1291 isl_set_free(C);
1292
1293 return isl_bool_ok(*res != NULL((void*)0));
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_setisl_map **dom = NULL((void*)0);
1312 isl_setisl_map **ran = NULL((void*)0);
1313 int *left = NULL((void*)0);
1314 int *right = NULL((void*)0);
1315 isl_setisl_map *C;
1316 isl_size d;
1317 isl_map *res = NULL((void*)0);
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)((isl_map * *)isl_calloc_or_die(map->ctx, map->n, sizeof
(isl_map *)))
;
1334 ran = isl_calloc_array(map->ctx, isl_set *, map->n)((isl_map * *)isl_calloc_or_die(map->ctx, map->n, sizeof
(isl_map *)))
;
1335 left = isl_calloc_array(map->ctx, int, map->n)((int *)isl_calloc_or_die(map->ctx, map->n, sizeof(int)
))
;
1336 right = isl_calloc_array(map->ctx, int, map->n)((int *)isl_calloc_or_die(map->ctx, map->n, sizeof(int)
))
;
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 & LEFT2) ? left : NULL((void*)0),
1394 (comp & RIGHT1) ? right : NULL((void*)0));
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 & LEFT2) ? left : NULL((void*)0),
1403 (comp & RIGHT1) ? right : NULL((void*)0), &exact_i);
1404 if (!res)
1405 goto error;
1406 if (exact_i)
1407 break;
1408 isl_map_free(res);
1409 res = NULL((void*)0);
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((void*)0);
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_setisl_map **set, int *group, __isl_take isl_setisl_map *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((void*)0);
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)isl_sioimath_set_si((bmap->eq[k][0]), 1);
1503 isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1)isl_sioimath_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)isl_sioimath_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((void*)0);
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((void*)0));
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((void*)0);
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)((isl_map ** *)isl_calloc_or_die(map->ctx, n, sizeof(isl_map
**)))
;
1615 if (!grid)
1616 goto error;
1617 for (i = 0; i < n; ++i) {
1618 grid[i] = isl_calloc_array(map->ctx, isl_map *, n)((isl_map * *)isl_calloc_or_die(map->ctx, n, sizeof(isl_map
*)))
;
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((void*)0);
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_setisl_map ***set, int *n_group)
1682{
1683 int i;
1684 int *group = NULL((void*)0);
1685 int g;
1686
1687 *set = isl_calloc_array(ctx, isl_set *, 2 * n)((isl_map * *)isl_calloc_or_die(ctx, 2 * n, sizeof(isl_map *)
))
;
1688 group = isl_alloc_array(ctx, int, 2 * n)((int *)isl_malloc_or_die(ctx, (2 * n)*sizeof(int)));
1689
1690 if (!*set || !group)
1691 goto error;
1692
1693 for (i = 0; i < n; ++i) {
1694 isl_setisl_map *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((void*)0);
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((void*)0);
1725 }
1726 free(group);
1727 return NULL((void*)0);
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_setisl_map **set = NULL((void*)0);
1743 int *group = NULL((void*)0);
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((void*)0);
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((void*)0);
1800 struct isl_map *map21 = NULL((void*)0);
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((void*)0);
1885 struct isl_tc_follows_data data;
1886 struct isl_tarjan_graph *g = NULL((void*)0);
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((void*)0);
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((void*)0);
1992 isl_space *space = NULL((void*)0);
1993
1994 if (!map)
1995 return NULL((void*)0);
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((void*)0);
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_setisl_map *dom, __isl_take isl_setisl_map *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_setisl_map *delta;
2174 isl_map *app = NULL((void*)0);
2175 isl_basic_setisl_basic_map *aff = NULL((void*)0);
2176 isl_basic_map *bmap = NULL((void*)0);
2177 isl_vec *obj = NULL((void*)0);
2178 isl_int opt;
2179
2180 isl_int_init(opt)isl_sioimath_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)isl_sioimath_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)isl_sioimath_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)isl_sioimath_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((void*)0), NULL((void*)0));
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)isl_sioimath_set_si((bmap->ineq[k][1 + nparam + i]), -1);
2234 isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1)isl_sioimath_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)isl_sioimath_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((void*)0), NULL((void*)0));
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)isl_sioimath_set_si((bmap->ineq[k][1 + nparam + i]), 1);
2249 isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1)isl_sioimath_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)isl_sioimath_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)isl_sioimath_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)isl_sioimath_set_si((bmap->ineq[k][0]), -1);
2262 isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1)isl_sioimath_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)isl_sioimath_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)isl_sioimath_clear((opt));
2286 return NULL((void*)0);
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_setisl_map *domain;
2305 isl_setisl_map *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_setisl_map *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_setisl_map *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((void*)0);
2400 int j = -1;
2401 isl_setisl_map *C;
2402
2403 *tc = NULL((void*)0);
2404 *qc = NULL((void*)0);
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((void*)0);
2450 *tc = NULL((void*)0);
2451
2452 return 0;
2453error:
2454 isl_map_free(id);
2455 isl_map_free(*qc);
2456 isl_map_free(*tc);
2457 *qc = NULL((void*)0);
2458 *tc = NULL((void*)0);
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((void*)0);
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((void*)0)));
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((void*)0);
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_BOX1)
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((void*)0);
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((void*)0);
2634 isl_setisl_map **set = NULL((void*)0);
2635 isl_map ***grid = NULL((void*)0);
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)((isl_map ** *)isl_calloc_or_die(ctx, n_group, sizeof(isl_map
**)))
;
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)((isl_map * *)isl_calloc_or_die(ctx, n_group, sizeof(isl_map *
)))
;
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((void*)0);
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((void*)0);
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)((isl_basic_map * *)isl_calloc_or_die(ctx, n, sizeof(isl_basic_map
*)))
;
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((void*)0);
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((void*)0);
2762 isl_basic_map **next;
2763 isl_union_map *path = NULL((void*)0);
2764 struct isl_tc_follows_data data;
2765 struct isl_tarjan_graph *g = NULL((void*)0);
2766 int c, l;
2767 int recheck = 0;
2768
2769 n = 0;
2770 if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
19
Assuming the condition is false
20
Taking false branch
2771 goto error;
2772
2773 if (n == 0)
21
Assuming 'n' is not equal to 0
22
Taking false branch
2774 return umap;
2775 if (n <= 1)
23
Assuming 'n' is > 1
24
Taking false branch
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)((isl_basic_map * *)isl_calloc_or_die(ctx, n, sizeof(isl_basic_map
*)))
;
2780 if (!list)
25
Assuming 'list' is non-null
26
Taking false branch
2781 goto error;
2782
2783 next = list;
2784 if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
27
Assuming the condition is false
28
Taking false branch
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)
29
Assuming 'g' is non-null
30
Taking false branch
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) {
31
Loop condition is true. Entering loop body
34
Loop condition is true. Entering loop body
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) {
32
Assuming the condition is false
33
Loop condition is false. Execution continues on line 2808
35
Assuming the condition is true
36
Loop condition is true. Entering loop body
37
Assuming the condition is false
38
Loop condition is false. Execution continues on line 2808
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
38.1
'c' is > 1
> 1 && data.check_closed && !*exact) {
39
Assuming field 'check_closed' is not equal to 0
40
Dereference of null pointer (loaded from variable '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((void*)0);
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)
9
Assuming 'umap' is non-null
10
Taking false branch
2861 return NULL((void*)0);
2862
2863 if (exact)
11
Assuming 'exact' is null
12
Taking false branch
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)
13
Assuming 'closed' is >= 0
14
Taking false branch
2870 goto error;
2871 if (closed)
15
Assuming 'closed' is 0
16
Taking false branch
2872 return umap;
2873 umap = union_components(umap, exact);
17
Passing null pointer value via 2nd parameter 'exact'
18
Calling 'union_components'
2874 return umap;
2875error:
2876 isl_union_map_free(umap);
2877 return NULL((void*)0);
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)
1
Assuming 'n' is >= 0
2
Taking false branch
2922 return isl_union_map_free(umap);
2923 if (n == 0)
3
Assuming 'n' is not equal to 0
4
Taking false branch
2924 return umap;
2925 if (n == 1) {
5
Assuming 'n' is not equal to 1
6
Taking false branch
2926 struct isl_union_power up = { NULL((void*)0), 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);
7
Passing value via 2nd parameter 'exact'
8
Calling 'isl_union_map_transitive_closure'
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 TYPEisl_union_map
2942#define TYPEisl_union_map isl_map
2943#include "isl_power_templ.c"
2944
2945#undef TYPEisl_union_map
2946#define TYPEisl_union_map isl_union_map
2947#include "isl_power_templ.c"