Bug Summary

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

Annotated Source Code

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