Coverage Report

Created: 2017-11-23 03:11

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/llvm/tools/polly/lib/External/isl/isl_transitive_closure.c
Line
Count
Source (jump to first uncovered line)
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
23
int isl_map_is_transitively_closed(__isl_keep isl_map *map)
24
115
{
25
115
  isl_map *map2;
26
115
  int closed;
27
115
28
115
  map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));
29
115
  closed = isl_map_is_subset(map2, map);
30
115
  isl_map_free(map2);
31
115
32
115
  return closed;
33
115
}
34
35
int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)
36
72
{
37
72
  isl_union_map *umap2;
38
72
  int closed;
39
72
40
72
  umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),
41
72
            isl_union_map_copy(umap));
42
72
  closed = isl_union_map_is_subset(umap2, umap);
43
72
  isl_union_map_free(umap2);
44
72
45
72
  return closed;
46
72
}
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
 */
54
static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
55
  int exactly, int length)
56
310
{
57
310
  isl_space *dim;
58
310
  struct isl_basic_map *bmap;
59
310
  unsigned d;
60
310
  unsigned nparam;
61
310
  int k;
62
310
  isl_int *c;
63
310
64
310
  if (!map)
65
0
    return NULL;
66
310
67
310
  dim = isl_map_get_space(map);
68
310
  d = isl_space_dim(dim, isl_dim_in);
69
310
  nparam = isl_space_dim(dim, isl_dim_param);
70
310
  bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
71
310
  if (exactly) {
72
6
    k = isl_basic_map_alloc_equality(bmap);
73
6
    if (k < 0)
74
0
      goto error;
75
6
    c = bmap->eq[k];
76
304
  } else {
77
304
    k = isl_basic_map_alloc_inequality(bmap);
78
304
    if (k < 0)
79
0
      goto error;
80
304
    c = bmap->ineq[k];
81
304
  }
82
310
  isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
83
310
  isl_int_set_si(c[0], -length);
84
310
  isl_int_set_si(c[1 + nparam + d - 1], -1);
85
310
  isl_int_set_si(c[1 + nparam + d + d - 1], 1);
86
310
87
310
  bmap = isl_basic_map_finalize(bmap);
88
310
  map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
89
310
90
310
  return map;
91
0
error:
92
0
  isl_basic_map_free(bmap);
93
0
  isl_map_free(map);
94
0
  return NULL;
95
310
}
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
 */
114
static int check_power_exactness(__isl_take isl_map *map,
115
  __isl_take isl_map *app)
116
3
{
117
3
  int exact;
118
3
  isl_map *app_1;
119
3
  isl_map *app_2;
120
3
121
3
  map = isl_map_add_dims(map, isl_dim_in, 1);
122
3
  map = isl_map_add_dims(map, isl_dim_out, 1);
123
3
  map = set_path_length(map, 1, 1);
124
3
125
3
  app_1 = set_path_length(isl_map_copy(app), 1, 1);
126
3
127
3
  exact = isl_map_is_subset(app_1, map);
128
3
  isl_map_free(app_1);
129
3
130
3
  if (!exact || exact < 0) {
131
0
    isl_map_free(app);
132
0
    isl_map_free(map);
133
0
    return exact;
134
0
  }
135
3
136
3
  app_1 = set_path_length(isl_map_copy(app), 0, 1);
137
3
  app_2 = set_path_length(app, 0, 2);
138
3
  app_1 = isl_map_apply_range(map, app_1);
139
3
140
3
  exact = isl_map_is_subset(app_2, app_1);
141
3
142
3
  isl_map_free(app_1);
143
3
  isl_map_free(app_2);
144
3
145
3
  return exact;
146
3
}
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
 */
169
static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
170
  int project)
171
105
{
172
105
  isl_map *test;
173
105
  int exact;
174
105
  unsigned d;
175
105
176
105
  if (!project)
177
3
    return check_power_exactness(map, app);
178
102
179
102
  d = isl_map_dim(map, isl_dim_in);
180
102
  app = set_path_length(app, 0, 1);
181
102
  app = isl_map_project_out(app, isl_dim_in, d, 1);
182
102
  app = isl_map_project_out(app, isl_dim_out, d, 1);
183
102
184
102
  app = isl_map_reset_space(app, isl_map_get_space(map));
185
102
186
102
  test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
187
102
  test = isl_map_union(test, isl_map_copy(map));
188
102
189
102
  exact = isl_map_is_subset(app, test);
190
102
191
102
  isl_map_free(app);
192
102
  isl_map_free(test);
193
102
194
102
  isl_map_free(map);
195
102
196
102
  return exact;
197
102
}
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
 */
218
static __isl_give isl_map *path_along_steps(__isl_take isl_space *dim,
219
  __isl_keep isl_mat *steps)
220
189
{
221
189
  int i, j, k;
222
189
  struct isl_basic_map *path = NULL;
223
189
  unsigned d;
224
189
  unsigned n;
225
189
  unsigned nparam;
226
189
227
189
  if (!dim || !steps)
228
0
    goto error;
229
189
230
189
  d = isl_space_dim(dim, isl_dim_in);
231
189
  n = steps->n_row;
232
189
  nparam = isl_space_dim(dim, isl_dim_param);
233
189
234
189
  path = isl_basic_map_alloc_space(isl_space_copy(dim), n, d, n);
235
189
236
400
  for (i = 0; i < n; 
++i211
) {
237
211
    k = isl_basic_map_alloc_div(path);
238
211
    if (k < 0)
239
0
      goto error;
240
211
    isl_assert(steps->ctx, i == k, goto error);
241
211
    isl_int_set_si(path->div[k][0], 0);
242
211
  }
243
189
244
928
  
for (i = 0; 189
i < d;
++i739
) {
245
739
    k = isl_basic_map_alloc_equality(path);
246
739
    if (k < 0)
247
0
      goto error;
248
739
    isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
249
739
    isl_int_set_si(path->eq[k][1 + nparam + i], 1);
250
739
    isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
251
739
    if (i == d - 1)
252
400
      
for (j = 0; 189
j < n;
++j211
)
253
211
        isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
254
739
    else
255
1.16k
      
for (j = 0; 550
j < n;
++j617
)
256
617
        isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
257
739
              steps->row[j][i]);
258
739
  }
259
189
260
400
  
for (i = 0; 189
i < n;
++i211
) {
261
211
    k = isl_basic_map_alloc_inequality(path);
262
211
    if (k < 0)
263
0
      goto error;
264
211
    isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
265
211
    isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
266
211
  }
267
189
268
189
  isl_space_free(dim);
269
189
270
189
  path = isl_basic_map_simplify(path);
271
189
  path = isl_basic_map_finalize(path);
272
189
  return isl_map_from_basic_map(path);
273
0
error:
274
0
  isl_space_free(dim);
275
0
  isl_basic_map_free(path);
276
0
  return NULL;
277
189
}
278
279
120
#define IMPURE    0
280
205
#define PURE_PARAM  1
281
457
#define PURE_VAR  2
282
44
#define MIXED   3
283
284
/* Check whether the parametric constant term of constraint c is never
285
 * positive in "bset".
286
 */
287
static isl_bool parametric_constant_never_positive(
288
  __isl_keep isl_basic_set *bset, isl_int *c, int *div_purity)
289
14
{
290
14
  unsigned d;
291
14
  unsigned n_div;
292
14
  unsigned nparam;
293
14
  int i;
294
14
  int k;
295
14
  isl_bool empty;
296
14
297
14
  n_div = isl_basic_set_dim(bset, isl_dim_div);
298
14
  d = isl_basic_set_dim(bset, isl_dim_set);
299
14
  nparam = isl_basic_set_dim(bset, isl_dim_param);
300
14
301
14
  bset = isl_basic_set_copy(bset);
302
14
  bset = isl_basic_set_cow(bset);
303
14
  bset = isl_basic_set_extend_constraints(bset, 0, 1);
304
14
  k = isl_basic_set_alloc_inequality(bset);
305
14
  if (k < 0)
306
0
    goto error;
307
14
  isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
308
14
  isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
309
27
  for (i = 0; i < n_div; 
++i13
) {
310
13
    if (div_purity[i] != PURE_PARAM)
311
13
      
continue1
;
312
12
    isl_int_set(bset->ineq[k][1 + nparam + d + i],
313
12
          c[1 + nparam + d + i]);
314
12
  }
315
14
  isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
316
14
  empty = isl_basic_set_is_empty(bset);
317
14
  isl_basic_set_free(bset);
318
14
319
14
  return empty;
320
0
error:
321
0
  isl_basic_set_free(bset);
322
0
  return isl_bool_error;
323
14
}
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
 */
332
static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,
333
  int eq)
334
69
{
335
69
  unsigned d;
336
69
  unsigned n_div;
337
69
  unsigned nparam;
338
69
  isl_bool empty;
339
69
  int i;
340
69
  int p = 0, v = 0;
341
69
342
69
  n_div = isl_basic_set_dim(bset, isl_dim_div);
343
69
  d = isl_basic_set_dim(bset, isl_dim_set);
344
69
  nparam = isl_basic_set_dim(bset, isl_dim_param);
345
69
346
157
  for (i = 0; i < n_div; 
++i88
) {
347
88
    if (isl_int_is_zero(c[1 + nparam + d + i]))
348
88
      
continue75
;
349
13
    switch (div_purity[i]) {
350
13
    
case 12
PURE_PARAM12
: p = 1; break;
351
13
    
case 1
PURE_VAR1
: v = 1; break;
352
13
    
default: return 0
IMPURE0
;
353
88
    }
354
88
  }
355
69
  if (!p && 
isl_seq_first_non_zero(c + 1, nparam) == -157
)
356
40
    return PURE_VAR;
357
29
  if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1)
358
15
    return PURE_PARAM;
359
14
360
14
  empty = parametric_constant_never_positive(bset, c, div_purity);
361
14
  if (eq && 
empty >= 00
&&
!empty0
) {
362
0
    isl_seq_neg(c, c, 1 + nparam + d + n_div);
363
0
    empty = parametric_constant_never_positive(bset, c, div_purity);
364
0
  }
365
14
366
14
  return empty < 0 ? 
-10
: empty ?
MIXED1
:
IMPURE13
;
367
69
}
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
 */
376
static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)
377
10
{
378
10
  int i, j;
379
10
  int *div_purity;
380
10
  unsigned d;
381
10
  unsigned n_div;
382
10
  unsigned nparam;
383
10
384
10
  if (!bset)
385
0
    return NULL;
386
10
387
10
  n_div = isl_basic_set_dim(bset, isl_dim_div);
388
10
  d = isl_basic_set_dim(bset, isl_dim_set);
389
10
  nparam = isl_basic_set_dim(bset, isl_dim_param);
390
10
391
10
  div_purity = isl_alloc_array(bset->ctx, int, n_div);
392
10
  if (n_div && 
!div_purity4
)
393
0
    return NULL;
394
10
395
17
  
for (i = 0; 10
i < bset->n_div;
++i7
) {
396
7
    int p = 0, v = 0;
397
7
    if (isl_int_is_zero(bset->div[i][0])) {
398
0
      div_purity[i] = IMPURE;
399
0
      continue;
400
0
    }
401
7
    if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1)
402
6
      p = 1;
403
7
    if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)
404
1
      v = 1;
405
10
    for (j = 0; j < i; 
++j3
) {
406
3
      if (isl_int_is_zero(bset->div[i][2 + nparam + d + j]))
407
3
        continue;
408
0
      switch (div_purity[j]) {
409
0
      case PURE_PARAM: p = 1; break;
410
0
      case PURE_VAR: v = 1; break;
411
0
      default: p = v = 1; break;
412
3
      }
413
3
    }
414
7
    div_purity[i] = v ? 
p ? 1
IMPURE0
:
PURE_VAR1
:
PURE_PARAM6
;
415
7
  }
416
10
417
10
  return div_purity;
418
10
}
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
 */
424
static int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
425
10
{
426
10
  isl_basic_map *test = NULL;
427
10
  isl_basic_map *id = NULL;
428
10
  int k;
429
10
  int is_id;
430
10
431
10
  test = isl_basic_map_copy(path);
432
10
  test = isl_basic_map_extend_constraints(test, 1, 0);
433
10
  k = isl_basic_map_alloc_equality(test);
434
10
  if (k < 0)
435
0
    goto error;
436
10
  isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
437
10
  isl_int_set_si(test->eq[k][pos], 1);
438
10
  test = isl_basic_map_gauss(test, NULL);
439
10
  id = isl_basic_map_identity(isl_basic_map_get_space(path));
440
10
  is_id = isl_basic_map_is_equal(test, id);
441
10
  isl_basic_map_free(test);
442
10
  isl_basic_map_free(id);
443
10
  return is_id;
444
0
error:
445
0
  isl_basic_map_free(test);
446
0
  return -1;
447
10
}
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
 */
455
static __isl_give isl_basic_map *add_delta_constraints(
456
  __isl_take isl_basic_map *path,
457
  __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
458
  unsigned d, int *div_purity, int eq, int *impurity)
459
38
{
460
38
  int i, k;
461
38
  int n = eq ? 
delta->n_eq19
:
delta->n_ineq19
;
462
38
  isl_int **delta_c = eq ? 
delta->eq19
:
delta->ineq19
;
463
38
  unsigned n_div;
464
38
465
38
  n_div = isl_basic_set_dim(delta, isl_dim_div);
466
38
467
145
  for (i = 0; i < n; 
++i107
) {
468
107
    isl_int *path_c;
469
107
    int p = PURE_VAR;
470
107
    if (impurity)
471
69
      p = purity(delta, delta_c[i], div_purity, eq);
472
107
    if (p < 0)
473
0
      goto error;
474
107
    if (p != PURE_VAR && 
p != 29
PURE_PARAM29
&&
!*impurity14
)
475
9
      *impurity = 1;
476
107
    if (p == IMPURE)
477
107
      
continue13
;
478
94
    if (eq && 
p != 43
MIXED43
) {
479
43
      k = isl_basic_map_alloc_equality(path);
480
43
      if (k < 0)
481
0
        goto error;
482
43
      path_c = path->eq[k];
483
51
    } else {
484
51
      k = isl_basic_map_alloc_inequality(path);
485
51
      if (k < 0)
486
0
        goto error;
487
51
      path_c = path->ineq[k];
488
51
    }
489
94
    isl_seq_clr(path_c, 1 + isl_basic_map_total_dim(path));
490
94
    if (p == PURE_VAR) {
491
78
      isl_seq_cpy(path_c + off,
492
78
            delta_c[i] + 1 + nparam, d);
493
78
      isl_int_set(path_c[off + d], delta_c[i][0]);
494
78
    } else 
if (16
p == 16
PURE_PARAM16
) {
495
15
      isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
496
15
    } else {
497
1
      isl_seq_cpy(path_c + off,
498
1
            delta_c[i] + 1 + nparam, d);
499
1
      isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
500
1
    }
501
107
    isl_seq_cpy(path_c + off - n_div,
502
107
          delta_c[i] + 1 + nparam + d, n_div);
503
107
  }
504
38
505
38
  return path;
506
0
error:
507
0
  isl_basic_map_free(path);
508
0
  return NULL;
509
38
}
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
 */
566
static __isl_give isl_map *path_along_delta(__isl_take isl_space *dim,
567
  __isl_take isl_basic_set *delta)
568
10
{
569
10
  isl_basic_map *path = NULL;
570
10
  unsigned d;
571
10
  unsigned n_div;
572
10
  unsigned nparam;
573
10
  unsigned off;
574
10
  int i, k;
575
10
  int is_id;
576
10
  int *div_purity = NULL;
577
10
  int impurity = 0;
578
10
579
10
  if (!delta)
580
0
    goto error;
581
10
  n_div = isl_basic_set_dim(delta, isl_dim_div);
582
10
  d = isl_basic_set_dim(delta, isl_dim_set);
583
10
  nparam = isl_basic_set_dim(delta, isl_dim_param);
584
10
  path = isl_basic_map_alloc_space(isl_space_copy(dim), n_div + d + 1,
585
10
      d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);
586
10
  off = 1 + nparam + 2 * (d + 1) + n_div;
587
10
588
61
  for (i = 0; i < n_div + d + 1; 
++i51
) {
589
51
    k = isl_basic_map_alloc_div(path);
590
51
    if (k < 0)
591
0
      goto error;
592
51
    isl_int_set_si(path->div[k][0], 0);
593
51
  }
594
10
595
54
  
for (i = 0; 10
i < d + 1;
++i44
) {
596
44
    k = isl_basic_map_alloc_equality(path);
597
44
    if (k < 0)
598
0
      goto error;
599
44
    isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
600
44
    isl_int_set_si(path->eq[k][1 + nparam + i], 1);
601
44
    isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
602
44
    isl_int_set_si(path->eq[k][off + i], 1);
603
44
  }
604
10
605
10
  div_purity = get_div_purity(delta);
606
10
  if (n_div && 
!div_purity4
)
607
0
    goto error;
608
10
609
10
  path = add_delta_constraints(path, delta, off, nparam, d,
610
10
             div_purity, 1, &impurity);
611
10
  path = add_delta_constraints(path, delta, off, nparam, d,
612
10
             div_purity, 0, &impurity);
613
10
  if (impurity) {
614
9
    isl_space *dim = isl_basic_set_get_space(delta);
615
9
    delta = isl_basic_set_project_out(delta,
616
9
              isl_dim_param, 0, nparam);
617
9
    delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam);
618
9
    delta = isl_basic_set_reset_space(delta, dim);
619
9
    if (!delta)
620
0
      goto error;
621
9
    path = isl_basic_map_extend_constraints(path, delta->n_eq,
622
9
              delta->n_ineq + 1);
623
9
    path = add_delta_constraints(path, delta, off, nparam, d,
624
9
               NULL, 1, NULL);
625
9
    path = add_delta_constraints(path, delta, off, nparam, d,
626
9
               NULL, 0, NULL);
627
9
    path = isl_basic_map_gauss(path, NULL);
628
9
  }
629
10
630
10
  is_id = empty_path_is_identity(path, off + d);
631
10
  if (is_id < 0)
632
0
    goto error;
633
10
634
10
  k = isl_basic_map_alloc_inequality(path);
635
10
  if (k < 0)
636
0
    goto error;
637
10
  isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
638
10
  if (!is_id)
639
10
    
isl_int_set_si8
(path->ineq[k][0], -1);
640
10
  isl_int_set_si(path->ineq[k][off + d], 1);
641
10
      
642
10
  free(div_purity);
643
10
  isl_basic_set_free(delta);
644
10
  path = isl_basic_map_finalize(path);
645
10
  if (is_id) {
646
2
    isl_space_free(dim);
647
2
    return isl_map_from_basic_map(path);
648
2
  }
649
8
  return isl_basic_map_union(path, isl_basic_map_identity(dim));
650
0
error:
651
0
  free(div_purity);
652
0
  isl_space_free(dim);
653
0
  isl_basic_set_free(delta);
654
0
  isl_basic_map_free(path);
655
0
  return NULL;
656
10
}
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
 */
665
static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_space *dim,
666
  unsigned param)
667
3
{
668
3
  struct isl_basic_map *bmap;
669
3
  unsigned d;
670
3
  unsigned nparam;
671
3
  int k;
672
3
673
3
  d = isl_space_dim(dim, isl_dim_in);
674
3
  nparam = isl_space_dim(dim, isl_dim_param);
675
3
  bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
676
3
  k = isl_basic_map_alloc_equality(bmap);
677
3
  if (k < 0)
678
0
    goto error;
679
3
  isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
680
3
  isl_int_set_si(bmap->eq[k][1 + param], -1);
681
3
  isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
682
3
  isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
683
3
684
3
  k = isl_basic_map_alloc_inequality(bmap);
685
3
  if (k < 0)
686
0
    goto error;
687
3
  isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
688
3
  isl_int_set_si(bmap->ineq[k][1 + param], 1);
689
3
  isl_int_set_si(bmap->ineq[k][0], -1);
690
3
691
3
  bmap = isl_basic_map_finalize(bmap);
692
3
  return isl_map_from_basic_map(bmap);
693
0
error:
694
0
  isl_basic_map_free(bmap);
695
0
  return NULL;
696
3
}
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
 */
707
static int is_acyclic(__isl_take isl_map *path)
708
102
{
709
102
  int i;
710
102
  int acyclic;
711
102
  unsigned dim;
712
102
  struct isl_set *delta;
713
102
714
102
  delta = isl_map_deltas(path);
715
102
  dim = isl_set_dim(delta, isl_dim_set);
716
513
  for (i = 0; i < dim; 
++i411
) {
717
411
    if (i == dim -1)
718
102
      delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
719
309
    else
720
309
      delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
721
411
  }
722
102
723
102
  acyclic = isl_set_is_empty(delta);
724
102
  isl_set_free(delta);
725
102
726
102
  return acyclic;
727
102
}
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
 */
754
static __isl_give isl_map *construct_extended_path(__isl_take isl_space *dim,
755
  __isl_keep isl_map *map, int *project)
756
198
{
757
198
  struct isl_mat *steps = NULL;
758
198
  struct isl_map *path = NULL;
759
198
  unsigned d;
760
198
  int i, j, n;
761
198
762
198
  if (!map)
763
0
    goto error;
764
198
765
198
  d = isl_map_dim(map, isl_dim_in);
766
198
767
198
  path = isl_map_identity(isl_space_copy(dim));
768
198
769
198
  steps = isl_mat_alloc(map->ctx, map->n, d);
770
198
  if (!steps)
771
0
    goto error;
772
198
773
198
  n = 0;
774
419
  for (i = 0; i < map->n; 
++i221
) {
775
221
    struct isl_basic_set *delta;
776
221
777
221
    delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
778
221
779
843
    for (j = 0; j < d; 
++j622
) {
780
632
      isl_bool fixed;
781
632
782
632
      fixed = isl_basic_set_plain_dim_is_fixed(delta, j,
783
632
                  &steps->row[n][j]);
784
632
      if (fixed < 0) {
785
0
        isl_basic_set_free(delta);
786
0
        goto error;
787
0
      }
788
632
      if (!fixed)
789
10
        break;
790
632
    }
791
221
792
221
793
221
    if (j < d) {
794
10
      path = isl_map_apply_range(path,
795
10
        path_along_delta(isl_space_copy(dim), delta));
796
10
      path = isl_map_coalesce(path);
797
211
    } else {
798
211
      isl_basic_set_free(delta);
799
211
      ++n;
800
211
    }
801
221
  }
802
198
803
198
  if (n > 0) {
804
189
    steps->n_row = n;
805
189
    path = isl_map_apply_range(path,
806
189
        path_along_steps(isl_space_copy(dim), steps));
807
189
  }
808
198
809
198
  if (project && 
*project105
) {
810
102
    *project = is_acyclic(isl_map_copy(path));
811
102
    if (*project < 0)
812
0
      goto error;
813
198
  }
814
198
815
198
  isl_space_free(dim);
816
198
  isl_mat_free(steps);
817
198
  return path;
818
0
error:
819
0
  isl_space_free(dim);
820
0
  isl_mat_free(steps);
821
0
  isl_map_free(path);
822
0
  return NULL;
823
198
}
824
825
static isl_bool isl_set_overlaps(__isl_keep isl_set *set1,
826
  __isl_keep isl_set *set2)
827
417
{
828
417
  isl_set *i;
829
417
  isl_bool no_overlap;
830
417
831
417
  if (!set1 || !set2)
832
0
    return isl_bool_error;
833
417
834
417
  if (!isl_space_tuple_is_equal(set1->dim, isl_dim_set,
835
417
          set2->dim, isl_dim_set))
836
0
    return isl_bool_false;
837
417
838
417
  i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2));
839
417
  no_overlap = isl_set_is_empty(i);
840
417
  isl_set_free(i);
841
417
842
417
  return isl_bool_not(no_overlap);
843
417
}
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
 */
864
static __isl_give isl_map *construct_component(__isl_take isl_space *dim,
865
  __isl_keep isl_map *map, int *exact, int project)
866
112
{
867
112
  struct isl_set *domain = NULL;
868
112
  struct isl_set *range = NULL;
869
112
  struct isl_map *app = NULL;
870
112
  struct isl_map *path = NULL;
871
112
  isl_bool overlaps;
872
112
873
112
  domain = isl_map_domain(isl_map_copy(map));
874
112
  domain = isl_set_coalesce(domain);
875
112
  range = isl_map_range(isl_map_copy(map));
876
112
  range = isl_set_coalesce(range);
877
112
  overlaps = isl_set_overlaps(domain, range);
878
112
  if (overlaps < 0 || !overlaps) {
879
0
    isl_set_free(domain);
880
0
    isl_set_free(range);
881
0
    isl_space_free(dim);
882
0
883
0
    if (overlaps < 0)
884
0
      map = NULL;
885
0
    map = isl_map_copy(map);
886
0
    map = isl_map_add_dims(map, isl_dim_in, 1);
887
0
    map = isl_map_add_dims(map, isl_dim_out, 1);
888
0
    map = set_path_length(map, 1, 1);
889
0
    return map;
890
0
  }
891
112
  app = isl_map_from_domain_and_range(domain, range);
892
112
  app = isl_map_add_dims(app, isl_dim_in, 1);
893
112
  app = isl_map_add_dims(app, isl_dim_out, 1);
894
112
895
112
  path = construct_extended_path(isl_space_copy(dim), map,
896
112
          exact && 
*exact19
?
&project19
: NULL);
897
112
  app = isl_map_intersect(app, path);
898
112
899
112
  if (exact && 
*exact19
&&
900
112
      (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
901
19
              project)) < 0)
902
0
    goto error;
903
112
904
112
  isl_space_free(dim);
905
112
  app = set_path_length(app, 0, 1);
906
112
  return app;
907
0
error:
908
0
  isl_space_free(dim);
909
0
  isl_map_free(app);
910
0
  return NULL;
911
112
}
912
913
/* Call construct_component and, if "project" is set, project out
914
 * the final coordinates.
915
 */
916
static __isl_give isl_map *construct_projected_component(
917
  __isl_take isl_space *dim,
918
  __isl_keep isl_map *map, int *exact, int project)
919
112
{
920
112
  isl_map *app;
921
112
  unsigned d;
922
112
923
112
  if (!dim)
924
0
    return NULL;
925
112
  d = isl_space_dim(dim, isl_dim_in);
926
112
927
112
  app = construct_component(dim, map, exact, project);
928
112
  if (project) {
929
109
    app = isl_map_project_out(app, isl_dim_in, d - 1, 1);
930
109
    app = isl_map_project_out(app, isl_dim_out, d - 1, 1);
931
109
  }
932
112
  return app;
933
112
}
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
 */
940
static __isl_give isl_map *q_closure(__isl_take isl_space *dim,
941
  __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact)
942
86
{
943
86
  int project = 1;
944
86
  isl_map *path;
945
86
  isl_map *map;
946
86
  isl_map *app;
947
86
948
86
  dom = isl_set_add_dims(dom, isl_dim_set, 1);
949
86
  app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));
950
86
  map = isl_map_from_basic_map(isl_basic_map_copy(bmap));
951
86
  path = construct_extended_path(dim, map, &project);
952
86
  app = isl_map_intersect(app, path);
953
86
954
86
  if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0)
955
0
    goto error;
956
86
957
86
  return app;
958
0
error:
959
0
  isl_map_free(app);
960
0
  return NULL;
961
86
}
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
 */
966
static int has_spurious_elements(__isl_keep isl_map *qc,
967
  __isl_keep isl_set *dom, __isl_keep isl_set *ran)
968
84
{
969
84
  isl_set *s;
970
84
  int subset;
971
84
  unsigned d;
972
84
973
84
  if (!qc || !dom || !ran)
974
0
    return -1;
975
84
976
84
  d = isl_map_dim(qc, isl_dim_in);
977
84
978
84
  qc = isl_map_copy(qc);
979
84
  qc = set_path_length(qc, 0, 1);
980
84
  qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1);
981
84
  qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1);
982
84
983
84
  s = isl_map_domain(isl_map_copy(qc));
984
84
  subset = isl_set_is_subset(s, dom);
985
84
  isl_set_free(s);
986
84
  if (subset < 0)
987
0
    goto error;
988
84
  if (!subset) {
989
77
    isl_map_free(qc);
990
77
    return 1;
991
77
  }
992
7
993
7
  s = isl_map_range(qc);
994
7
  subset = isl_set_is_subset(s, ran);
995
7
  isl_set_free(s);
996
7
997
7
  return subset < 0 ? 
-10
: !subset;
998
0
error:
999
0
  isl_map_free(qc);
1000
0
  return -1;
1001
84
}
1002
1003
87
#define LEFT  2
1004
87
#define RIGHT 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
 */
1035
static int composability(__isl_keep isl_set *C, int i,
1036
  isl_set **dom, isl_set **ran, int *left, int *right,
1037
  __isl_keep isl_map *map)
1038
42
{
1039
42
  int j;
1040
42
  int ok;
1041
42
1042
42
  ok = LEFT | RIGHT;
1043
126
  for (j = 0; j < map->n && 
ok84
;
++j84
) {
1044
84
    isl_bool overlaps, subset;
1045
84
    if (j == i)
1046
42
      continue;
1047
42
1048
42
    if (ok & RIGHT) {
1049
42
      if (!dom[j])
1050
0
        dom[j] = isl_set_from_basic_set(
1051
0
          isl_basic_map_domain(
1052
0
            isl_basic_map_copy(map->p[j])));
1053
42
      if (!dom[j])
1054
0
        return -1;
1055
42
      overlaps = isl_set_overlaps(ran[i], dom[j]);
1056
42
      if (overlaps < 0)
1057
0
        return -1;
1058
42
      if (!overlaps)
1059
0
        right[j] = 0;
1060
42
      else {
1061
42
        subset = isl_set_is_subset(dom[j], C);
1062
42
        if (subset < 0)
1063
0
          return -1;
1064
42
        if (subset)
1065
41
          right[j] = 1;
1066
1
        else
1067
1
          ok &= ~RIGHT;
1068
42
      }
1069
42
    }
1070
42
1071
42
    if (ok & LEFT) {
1072
42
      if (!ran[j])
1073
0
        ran[j] = isl_set_from_basic_set(
1074
0
          isl_basic_map_range(
1075
0
            isl_basic_map_copy(map->p[j])));
1076
42
      if (!ran[j])
1077
0
        return -1;
1078
42
      overlaps = isl_set_overlaps(dom[i], ran[j]);
1079
42
      if (overlaps < 0)
1080
0
        return -1;
1081
42
      if (!overlaps)
1082
0
        left[j] = 0;
1083
42
      else {
1084
42
        subset = isl_set_is_subset(ran[j], C);
1085
42
        if (subset < 0)
1086
0
          return -1;
1087
42
        if (subset)
1088
41
          left[j] = 1;
1089
1
        else
1090
1
          ok &= ~LEFT;
1091
42
      }
1092
42
    }
1093
84
  }
1094
42
1095
42
  return ok;
1096
42
}
1097
1098
static __isl_give isl_map *anonymize(__isl_take isl_map *map)
1099
71
{
1100
71
  map = isl_map_reset(map, isl_dim_in);
1101
71
  map = isl_map_reset(map, isl_dim_out);
1102
71
  return map;
1103
71
}
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
 */
1109
static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
1110
  __isl_take isl_map *qc, int *left, int *right)
1111
7
{
1112
7
  int j;
1113
7
  isl_map *comp;
1114
7
1115
7
  comp = isl_map_empty(isl_map_get_space(map));
1116
21
  for (j = 0; j < map->n; 
++j14
) {
1117
14
    isl_map *map_j;
1118
14
1119
14
    if (j == i)
1120
7
      continue;
1121
7
1122
7
    map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));
1123
7
    map_j = anonymize(map_j);
1124
7
    if (left && left[j])
1125
7
      map_j = isl_map_apply_range(map_j, isl_map_copy(qc));
1126
7
    if (right && right[j])
1127
7
      map_j = isl_map_apply_range(isl_map_copy(qc), map_j);
1128
14
    comp = isl_map_union(comp, map_j);
1129
14
  }
1130
7
1131
7
  comp = isl_map_compute_divs(comp);
1132
7
  comp = isl_map_coalesce(comp);
1133
7
1134
7
  isl_map_free(qc);
1135
7
1136
7
  return comp;
1137
7
}
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
 */
1154
static __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
3
{
1158
3
  isl_map *map_i;
1159
3
  isl_map *tc;
1160
3
  isl_map *rtc = NULL;
1161
3
1162
3
  if (!map)
1163
0
    goto error;
1164
3
  isl_assert(map->ctx, left || right, goto error);
1165
3
1166
3
  map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
1167
3
  tc = construct_projected_component(isl_space_copy(dim), map_i,
1168
3
            exact, 1);
1169
3
  isl_map_free(map_i);
1170
3
1171
3
  if (*exact)
1172
3
    qc = isl_map_transitive_closure(qc, exact);
1173
3
1174
3
  if (!*exact) {
1175
0
    isl_space_free(dim);
1176
0
    isl_map_free(tc);
1177
0
    isl_map_free(qc);
1178
0
    return isl_map_universe(isl_map_get_space(map));
1179
0
  }
1180
3
1181
3
  if (!left || !right)
1182
0
    rtc = isl_map_union(isl_map_copy(tc),
1183
0
            isl_map_identity(isl_map_get_space(tc)));
1184
3
  if (!right)
1185
0
    qc = isl_map_apply_range(rtc, qc);
1186
3
  if (!left)
1187
0
    qc = isl_map_apply_range(qc, rtc);
1188
3
  qc = isl_map_union(tc, qc);
1189
3
1190
3
  isl_space_free(dim);
1191
3
1192
3
  return qc;
1193
0
error:
1194
0
  isl_space_free(dim);
1195
0
  isl_map_free(qc);
1196
0
  return NULL;
1197
3
}
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
 */
1214
static int incemental_on_entire_domain(__isl_keep isl_space *dim,
1215
  __isl_keep isl_map *map,
1216
  isl_set **dom, isl_set **ran, int *left, int *right,
1217
  __isl_give isl_map **res)
1218
24
{
1219
24
  int i;
1220
24
  isl_set *C;
1221
24
  unsigned d;
1222
24
1223
24
  *res = NULL;
1224
24
1225
24
  C = isl_set_union(isl_map_domain(isl_map_copy(map)),
1226
24
        isl_map_range(isl_map_copy(map)));
1227
24
  C = isl_set_from_basic_set(isl_set_simple_hull(C));
1228
24
  if (!C)
1229
0
    return -1;
1230
24
  if (C->n != 1) {
1231
0
    isl_set_free(C);
1232
0
    return 0;
1233
0
  }
1234
24
1235
24
  d = isl_map_dim(map, isl_dim_in);
1236
24
1237
66
  for (i = 0; i < map->n; 
++i42
) {
1238
45
    isl_map *qc;
1239
45
    int exact_i, spurious;
1240
45
    int j;
1241
45
    dom[i] = isl_set_from_basic_set(isl_basic_map_domain(
1242
45
          isl_basic_map_copy(map->p[i])));
1243
45
    ran[i] = isl_set_from_basic_set(isl_basic_map_range(
1244
45
          isl_basic_map_copy(map->p[i])));
1245
45
    qc = q_closure(isl_space_copy(dim), isl_set_copy(C),
1246
45
        map->p[i], &exact_i);
1247
45
    if (!qc)
1248
0
      goto error;
1249
45
    if (!exact_i) {
1250
1
      isl_map_free(qc);
1251
1
      continue;
1252
1
    }
1253
44
    spurious = has_spurious_elements(qc, dom[i], ran[i]);
1254
44
    if (spurious) {
1255
39
      isl_map_free(qc);
1256
39
      if (spurious < 0)
1257
0
        goto error;
1258
39
      continue;
1259
39
    }
1260
5
    qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1261
5
    qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1262
5
    qc = isl_map_compute_divs(qc);
1263
15
    for (j = 0; j < map->n; 
++j10
)
1264
10
      left[j] = right[j] = 1;
1265
5
    qc = compose(map, i, qc, left, right);
1266
5
    if (!qc)
1267
0
      goto error;
1268
5
    if (qc->n >= map->n) {
1269
2
      isl_map_free(qc);
1270
2
      continue;
1271
2
    }
1272
3
    *res = compute_incremental(isl_space_copy(dim), map, i, qc,
1273
3
        left, right, &exact_i);
1274
3
    if (!*res)
1275
0
      goto error;
1276
3
    if (exact_i)
1277
3
      break;
1278
0
    isl_map_free(*res);
1279
0
    *res = NULL;
1280
0
  }
1281
24
1282
24
  isl_set_free(C);
1283
24
1284
24
  return *res != NULL;
1285
24
error:
1286
0
  isl_set_free(C);
1287
0
  return -1;
1288
24
}
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
 */
1298
static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim,
1299
  __isl_keep isl_map *map, int *exact, int project)
1300
112
{
1301
112
  int i;
1302
112
  isl_set **dom = NULL;
1303
112
  isl_set **ran = NULL;
1304
112
  int *left = NULL;
1305
112
  int *right = NULL;
1306
112
  isl_set *C;
1307
112
  unsigned d;
1308
112
  isl_map *res = NULL;
1309
112
1310
112
  if (!project)
1311
3
    return construct_projected_component(dim, map, exact, project);
1312
109
1313
109
  if (!map)
1314
0
    goto error;
1315
109
  if (map->n <= 1)
1316
85
    return construct_projected_component(dim, map, exact, project);
1317
24
1318
24
  d = isl_map_dim(map, isl_dim_in);
1319
24
1320
24
  dom = isl_calloc_array(map->ctx, isl_set *, map->n);
1321
24
  ran = isl_calloc_array(map->ctx, isl_set *, map->n);
1322
24
  left = isl_calloc_array(map->ctx, int, map->n);
1323
24
  right = isl_calloc_array(map->ctx, int, map->n);
1324
24
  if (!ran || !dom || !left || !right)
1325
0
    goto error;
1326
24
1327
24
  if (incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 0)
1328
0
    goto error;
1329
24
1330
66
  
for (i = 0; 24
!res &&
i < map->n63
;
++i42
) {
1331
42
    isl_map *qc;
1332
42
    int exact_i, spurious, comp;
1333
42
    if (!dom[i])
1334
0
      dom[i] = isl_set_from_basic_set(
1335
0
          isl_basic_map_domain(
1336
0
            isl_basic_map_copy(map->p[i])));
1337
42
    if (!dom[i])
1338
0
      goto error;
1339
42
    if (!ran[i])
1340
0
      ran[i] = isl_set_from_basic_set(
1341
0
          isl_basic_map_range(
1342
0
            isl_basic_map_copy(map->p[i])));
1343
42
    if (!ran[i])
1344
0
      goto error;
1345
42
    C = isl_set_union(isl_set_copy(dom[i]),
1346
42
              isl_set_copy(ran[i]));
1347
42
    C = isl_set_from_basic_set(isl_set_simple_hull(C));
1348
42
    if (!C)
1349
0
      goto error;
1350
42
    if (C->n != 1) {
1351
0
      isl_set_free(C);
1352
0
      continue;
1353
0
    }
1354
42
    comp = composability(C, i, dom, ran, left, right, map);
1355
42
    if (!comp || 
comp < 041
) {
1356
1
      isl_set_free(C);
1357
1
      if (comp < 0)
1358
0
        goto error;
1359
1
      continue;
1360
1
    }
1361
41
    qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i);
1362
41
    if (!qc)
1363
0
      goto error;
1364
41
    if (!exact_i) {
1365
1
      isl_map_free(qc);
1366
1
      continue;
1367
1
    }
1368
40
    spurious = has_spurious_elements(qc, dom[i], ran[i]);
1369
40
    if (spurious) {
1370
38
      isl_map_free(qc);
1371
38
      if (spurious < 0)
1372
0
        goto error;
1373
38
      continue;
1374
38
    }
1375
2
    qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1376
2
    qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1377
2
    qc = isl_map_compute_divs(qc);
1378
2
    qc = compose(map, i, qc, (comp & LEFT) ? left : NULL,
1379
2
        (comp & RIGHT) ? right : NULL);
1380
2
    if (!qc)
1381
0
      goto error;
1382
2
    if (qc->n >= map->n) {
1383
2
      isl_map_free(qc);
1384
2
      continue;
1385
2
    }
1386
0
    res = compute_incremental(isl_space_copy(dim), map, i, qc,
1387
0
        (comp & LEFT) ? left : NULL,
1388
0
        (comp & RIGHT) ? right : NULL, &exact_i);
1389
0
    if (!res)
1390
0
      goto error;
1391
0
    if (exact_i)
1392
0
      break;
1393
0
    isl_map_free(res);
1394
0
    res = NULL;
1395
0
  }
1396
24
1397
72
  
for (i = 0; 24
i < map->n;
++i48
) {
1398
48
    isl_set_free(dom[i]);
1399
48
    isl_set_free(ran[i]);
1400
48
  }
1401
24
  free(dom);
1402
24
  free(ran);
1403
24
  free(left);
1404
24
  free(right);
1405
24
1406
24
  if (res) {
1407
3
    isl_space_free(dim);
1408
3
    return res;
1409
3
  }
1410
21
1411
21
  return construct_projected_component(dim, map, exact, project);
1412
0
error:
1413
0
  if (dom)
1414
0
    for (i = 0; i < map->n; ++i)
1415
0
      isl_set_free(dom[i]);
1416
0
  free(dom);
1417
0
  if (ran)
1418
0
    for (i = 0; i < map->n; ++i)
1419
0
      isl_set_free(ran[i]);
1420
0
  free(ran);
1421
0
  free(left);
1422
0
  free(right);
1423
0
  isl_space_free(dim);
1424
0
  return NULL;
1425
112
}
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
 */
1436
static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)
1437
336
{
1438
336
  int i;
1439
336
1440
336
  group[pos] = pos;
1441
336
  set[pos] = isl_set_copy(dom);
1442
336
1443
692
  for (i = pos - 1; i >= 0; 
--i356
) {
1444
356
    isl_bool o;
1445
356
1446
356
    if (group[i] != i)
1447
135
      continue;
1448
221
1449
221
    o = isl_set_overlaps(set[i], dom);
1450
221
    if (o < 0)
1451
0
      goto error;
1452
221
    if (!o)
1453
9
      continue;
1454
212
1455
212
    set[i] = isl_set_union(set[i], set[group[pos]]);
1456
212
    set[group[pos]] = NULL;
1457
212
    if (!set[i])
1458
0
      goto error;
1459
212
    group[group[pos]] = i;
1460
212
    group[pos] = i;
1461
212
  }
1462
336
1463
336
  isl_set_free(dom);
1464
336
  return 0;
1465
0
error:
1466
0
  isl_set_free(dom);
1467
0
  return -1;
1468
336
}
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
 */
1473
static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
1474
1
{
1475
1
  int i, j, k;
1476
1
  isl_space *dim;
1477
1
  isl_basic_map *bstep;
1478
1
  isl_map *step;
1479
1
  unsigned nparam;
1480
1
1481
1
  if (!map)
1482
0
    return -1;
1483
1
1484
1
  dim = isl_map_get_space(map);
1485
1
  nparam = isl_space_dim(dim, isl_dim_param);
1486
1
  dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in));
1487
1
  dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out));
1488
1
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
1489
1
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
1490
1
  bstep = isl_basic_map_alloc_space(dim, 0, 1, 0);
1491
1
  k = isl_basic_map_alloc_equality(bstep);
1492
1
  if (k < 0) {
1493
0
    isl_basic_map_free(bstep);
1494
0
    return -1;
1495
0
  }
1496
1
  isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep));
1497
1
  isl_int_set_si(bstep->eq[k][0], 1);
1498
1
  isl_int_set_si(bstep->eq[k][1 + nparam], 1);
1499
1
  isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1);
1500
1
  bstep = isl_basic_map_finalize(bstep);
1501
1
  step = isl_map_from_basic_map(bstep);
1502
1
1503
3
  for (i = 0; i < n; 
++i2
)
1504
6
    
for (j = 0; 2
j < n;
++j4
)
1505
4
      grid[i][j] = isl_map_product(grid[i][j],
1506
4
                 isl_map_copy(step));
1507
1
1508
1
  isl_map_free(step);
1509
1
1510
1
  return 0;
1511
1
}
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
 */
1524
static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)
1525
95
{
1526
95
  int r, p, q;
1527
95
1528
193
  for (r = 0; r < n; 
++r98
) {
1529
98
    int r_exact;
1530
98
    grid[r][r] = isl_map_transitive_closure(grid[r][r],
1531
98
        (exact && 
*exact4
) ?
&r_exact4
: NULL);
1532
98
    if (exact && 
*exact4
&&
!r_exact4
)
1533
0
      *exact = 0;
1534
98
1535
202
    for (p = 0; p < n; 
++p104
)
1536
220
      
for (q = 0; 104
q < n;
++q116
) {
1537
116
        isl_map *loop;
1538
116
        if (p == r && 
q == r104
)
1539
98
          continue;
1540
18
        loop = isl_map_apply_range(
1541
18
            isl_map_copy(grid[p][r]),
1542
18
            isl_map_copy(grid[r][q]));
1543
18
        grid[p][q] = isl_map_union(grid[p][q], loop);
1544
18
        loop = isl_map_apply_range(
1545
18
            isl_map_copy(grid[p][r]),
1546
18
          isl_map_apply_range(
1547
18
            isl_map_copy(grid[r][r]),
1548
18
            isl_map_copy(grid[r][q])));
1549
18
        grid[p][q] = isl_map_union(grid[p][q], loop);
1550
18
        grid[p][q] = isl_map_coalesce(grid[p][q]);
1551
18
      }
1552
98
  }
1553
95
}
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
 */
1574
static __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
28
{
1577
28
  int i, j, k;
1578
28
  isl_map ***grid = NULL;
1579
28
  isl_map *app;
1580
28
1581
28
  if (!map)
1582
0
    goto error;
1583
28
1584
28
  if (n == 1) {
1585
26
    free(group);
1586
26
    return incremental_closure(dim, map, exact, project);
1587
26
  }
1588
2
1589
2
  grid = isl_calloc_array(map->ctx, isl_map **, n);
1590
2
  if (!grid)
1591
0
    goto error;
1592
6
  
for (i = 0; 2
i < n;
++i4
) {
1593
4
    grid[i] = isl_calloc_array(map->ctx, isl_map *, n);
1594
4
    if (!grid[i])
1595
0
      goto error;
1596
12
    
for (j = 0; 4
j < n;
++j8
)
1597
8
      grid[i][j] = isl_map_empty(isl_map_get_space(map));
1598
4
  }
1599
2
1600
6
  
for (k = 0; 2
k < map->n;
++k4
) {
1601
4
    i = group[2 * k];
1602
4
    j = group[2 * k + 1];
1603
4
    grid[i][j] = isl_map_union(grid[i][j],
1604
4
        isl_map_from_basic_map(
1605
4
          isl_basic_map_copy(map->p[k])));
1606
4
  }
1607
2
1608
2
  if (!project && 
add_length(map, grid, n) < 01
)
1609
0
    goto error;
1610
2
1611
2
  floyd_warshall_iterate(grid, n, exact);
1612
2
1613
2
  app = isl_map_empty(isl_map_get_space(grid[0][0]));
1614
2
1615
6
  for (i = 0; i < n; 
++i4
) {
1616
12
    for (j = 0; j < n; 
++j8
)
1617
8
      app = isl_map_union(app, grid[i][j]);
1618
4
    free(grid[i]);
1619
4
  }
1620
2
  free(grid);
1621
2
1622
2
  free(group);
1623
2
  isl_space_free(dim);
1624
2
1625
2
  return app;
1626
0
error:
1627
0
  if (grid)
1628
0
    for (i = 0; i < n; ++i) {
1629
0
      if (!grid[i])
1630
0
        continue;
1631
0
      for (j = 0; j < n; ++j)
1632
0
        isl_map_free(grid[i][j]);
1633
0
      free(grid[i]);
1634
0
    }
1635
0
  free(grid);
1636
0
  free(group);
1637
0
  isl_space_free(dim);
1638
0
  return NULL;
1639
28
}
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
 */
1655
static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n,
1656
  isl_set ***set, int *n_group)
1657
121
{
1658
121
  int i;
1659
121
  int *group = NULL;
1660
121
  int g;
1661
121
1662
121
  *set = isl_calloc_array(ctx, isl_set *, 2 * n);
1663
121
  group = isl_alloc_array(ctx, int, 2 * n);
1664
121
1665
121
  if (!*set || !group)
1666
0
    goto error;
1667
121
1668
289
  
for (i = 0; 121
i < n;
++i168
) {
1669
168
    isl_set *dom;
1670
168
    dom = isl_set_from_basic_set(isl_basic_map_domain(
1671
168
        isl_basic_map_copy(list[i])));
1672
168
    if (merge(*set, group, dom, 2 * i) < 0)
1673
0
      goto error;
1674
168
    dom = isl_set_from_basic_set(isl_basic_map_range(
1675
168
        isl_basic_map_copy(list[i])));
1676
168
    if (merge(*set, group, dom, 2 * i + 1) < 0)
1677
0
      goto error;
1678
168
  }
1679
121
1680
121
  g = 0;
1681
457
  for (i = 0; i < 2 * n; 
++i336
)
1682
336
    if (group[i] == i) {
1683
124
      if (g != i) {
1684
0
        (*set)[g] = (*set)[i];
1685
0
        (*set)[i] = NULL;
1686
0
      }
1687
124
      group[i] = g++;
1688
124
    } else
1689
212
      group[i] = group[group[i]];
1690
121
1691
121
  *n_group = g;
1692
121
1693
121
  return group;
1694
0
error:
1695
0
  if (*set) {
1696
0
    for (i = 0; i < 2 * n; ++i)
1697
0
      isl_set_free((*set)[i]);
1698
0
    free(*set);
1699
0
    *set = NULL;
1700
0
  }
1701
0
  free(group);
1702
0
  return NULL;
1703
121
}
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
 */
1713
static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim,
1714
  __isl_keep isl_map *map, int *exact, int project)
1715
114
{
1716
114
  int i;
1717
114
  isl_set **set = NULL;
1718
114
  int *group = NULL;
1719
114
  int n;
1720
114
1721
114
  if (!map)
1722
0
    goto error;
1723
114
  if (map->n <= 1)
1724
86
    return incremental_closure(dim, map, exact, project);
1725
28
1726
28
  group = setup_groups(map->ctx, map->p, map->n, &set, &n);
1727
28
  if (!group)
1728
0
    goto error;
1729
28
1730
140
  
for (i = 0; 28
i < 2 * map->n;
++i112
)
1731
112
    isl_set_free(set[i]);
1732
28
1733
28
  free(set);
1734
28
1735
28
  return floyd_warshall_with_groups(dim, map, exact, project, group, n);
1736
0
error:
1737
0
  isl_space_free(dim);
1738
0
  return NULL;
1739
114
}
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
 */
1748
struct 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
 */
1771
static isl_bool basic_map_follows(int i, int j, void *user)
1772
481
{
1773
481
  struct isl_tc_follows_data *data = user;
1774
481
  struct isl_map *map12 = NULL;
1775
481
  struct isl_map *map21 = NULL;
1776
481
  isl_bool subset;
1777
481
1778
481
  if (!isl_space_tuple_is_equal(data->list[i]->dim, isl_dim_in,
1779
481
            data->list[j]->dim, isl_dim_out))
1780
377
    return isl_bool_false;
1781
104
1782
104
  map21 = isl_map_from_basic_map(
1783
104
      isl_basic_map_apply_range(
1784
104
        isl_basic_map_copy(data->list[j]),
1785
104
        isl_basic_map_copy(data->list[i])));
1786
104
  subset = isl_map_is_empty(map21);
1787
104
  if (subset < 0)
1788
0
    goto error;
1789
104
  if (subset) {
1790
3
    isl_map_free(map21);
1791
3
    return isl_bool_false;
1792
3
  }
1793
101
1794
101
  if (!isl_space_tuple_is_equal(data->list[i]->dim, isl_dim_in,
1795
101
            data->list[i]->dim, isl_dim_out) ||
1796
101
      !isl_space_tuple_is_equal(data->list[j]->dim, isl_dim_in,
1797
101
            data->list[j]->dim, isl_dim_out)) {
1798
0
    isl_map_free(map21);
1799
0
    return isl_bool_true;
1800
0
  }
1801
101
1802
101
  map12 = isl_map_from_basic_map(
1803
101
      isl_basic_map_apply_range(
1804
101
        isl_basic_map_copy(data->list[i]),
1805
101
        isl_basic_map_copy(data->list[j])));
1806
101
1807
101
  subset = isl_map_is_subset(map21, map12);
1808
101
1809
101
  isl_map_free(map12);
1810
101
  isl_map_free(map21);
1811
101
1812
101
  if (subset)
1813
6
    data->check_closed = 1;
1814
101
1815
101
  return subset < 0 ? 
isl_bool_error0
: !subset;
1816
0
error:
1817
0
  isl_map_free(map21);
1818
0
  return isl_bool_error;
1819
481
}
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
 */
1854
static __isl_give isl_map *construct_power_components(__isl_take isl_space *dim,
1855
  __isl_keep isl_map *map, int *exact, int project)
1856
110
{
1857
110
  int i, n, c;
1858
110
  struct isl_map *path = NULL;
1859
110
  struct isl_tc_follows_data data;
1860
110
  struct isl_tarjan_graph *g = NULL;
1861
110
  int *orig_exact;
1862
110
  int local_exact;
1863
110
1864
110
  if (!map)
1865
0
    goto error;
1866
110
  if (map->n <= 1)
1867
80
    return floyd_warshall(dim, map, exact, project);
1868
30
1869
30
  data.list = map->p;
1870
30
  data.check_closed = 0;
1871
30
  g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data);
1872
30
  if (!g)
1873
0
    goto error;
1874
30
1875
30
  orig_exact = exact;
1876
30
  if (data.check_closed && 
!exact4
)
1877
0
    exact = &local_exact;
1878
30
1879
30
  c = 0;
1880
30
  i = 0;
1881
30
  n = map->n;
1882
30
  if (project)
1883
27
    path = isl_map_empty(isl_map_get_space(map));
1884
3
  else
1885
3
    path = isl_map_empty(isl_space_copy(dim));
1886
30
  path = anonymize(path);
1887
64
  while (n) {
1888
34
    struct isl_map *comp;
1889
34
    isl_map *path_comp, *path_comb;
1890
34
    comp = isl_map_alloc_space(isl_map_get_space(map), n, 0);
1891
96
    while (g->order[i] != -1) {
1892
62
      comp = isl_map_add_basic_map(comp,
1893
62
            isl_basic_map_copy(map->p[g->order[i]]));
1894
62
      --n;
1895
62
      ++i;
1896
62
    }
1897
34
    path_comp = floyd_warshall(isl_space_copy(dim),
1898
34
            comp, exact, project);
1899
34
    path_comp = anonymize(path_comp);
1900
34
    path_comb = isl_map_apply_range(isl_map_copy(path),
1901
34
            isl_map_copy(path_comp));
1902
34
    path = isl_map_union(path, path_comp);
1903
34
    path = isl_map_union(path, path_comb);
1904
34
    isl_map_free(comp);
1905
34
    ++i;
1906
34
    ++c;
1907
34
  }
1908
30
1909
30
  if (c > 1 && 
data.check_closed4
&&
!*exact4
) {
1910
0
    int closed;
1911
0
1912
0
    closed = isl_map_is_transitively_closed(path);
1913
0
    if (closed < 0)
1914
0
      goto error;
1915
0
    if (!closed) {
1916
0
      isl_tarjan_graph_free(g);
1917
0
      isl_map_free(path);
1918
0
      return floyd_warshall(dim, map, orig_exact, project);
1919
0
    }
1920
30
  }
1921
30
1922
30
  isl_tarjan_graph_free(g);
1923
30
  isl_space_free(dim);
1924
30
1925
30
  return path;
1926
0
error:
1927
0
  isl_tarjan_graph_free(g);
1928
0
  isl_space_free(dim);
1929
0
  isl_map_free(path);
1930
0
  return NULL;
1931
110
}
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
 */
1962
static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
1963
  int *exact, int project)
1964
110
{
1965
110
  struct isl_map *app = NULL;
1966
110
  isl_space *dim = NULL;
1967
110
1968
110
  if (!map)
1969
0
    return NULL;
1970
110
1971
110
  dim = isl_map_get_space(map);
1972
110
1973
110
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
1974
110
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
1975
110
1976
110
  app = construct_power_components(isl_space_copy(dim), map,
1977
110
          exact, project);
1978
110
1979
110
  isl_space_free(dim);
1980
110
1981
110
  return app;
1982
110
}
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
 */
1992
static __isl_give isl_map *map_power(__isl_take isl_map *map,
1993
  int *exact, int project)
1994
110
{
1995
110
  struct isl_map *app = NULL;
1996
110
1997
110
  if (exact)
1998
17
    *exact = 1;
1999
110
2000
110
  if (!map)
2001
0
    return NULL;
2002
110
2003
110
  isl_assert(map->ctx,
2004
110
    isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
2005
110
    goto error);
2006
110
2007
110
  app = construct_power(map, exact, project);
2008
110
2009
110
  isl_map_free(map);
2010
110
  return app;
2011
0
error:
2012
0
  isl_map_free(map);
2013
0
  isl_map_free(app);
2014
0
  return NULL;
2015
110
}
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
3
{
2028
3
  isl_space *target_dim;
2029
3
  isl_space *dim;
2030
3
  isl_map *diff;
2031
3
  unsigned d;
2032
3
  unsigned param;
2033
3
2034
3
  if (!map)
2035
0
    return NULL;
2036
3
2037
3
  d = isl_map_dim(map, isl_dim_in);
2038
3
  param = isl_map_dim(map, isl_dim_param);
2039
3
2040
3
  map = isl_map_compute_divs(map);
2041
3
  map = isl_map_coalesce(map);
2042
3
2043
3
  if (isl_map_plain_is_empty(map)) {
2044
0
    map = isl_map_from_range(isl_map_wrap(map));
2045
0
    map = isl_map_add_dims(map, isl_dim_in, 1);
2046
0
    map = isl_map_set_dim_name(map, isl_dim_in, 0, "k");
2047
0
    return map;
2048
0
  }
2049
3
2050
3
  target_dim = isl_map_get_space(map);
2051
3
  target_dim = isl_space_from_range(isl_space_wrap(target_dim));
2052
3
  target_dim = isl_space_add_dims(target_dim, isl_dim_in, 1);
2053
3
  target_dim = isl_space_set_dim_name(target_dim, isl_dim_in, 0, "k");
2054
3
2055
3
  map = map_power(map, exact, 0);
2056
3
2057
3
  map = isl_map_add_dims(map, isl_dim_param, 1);
2058
3
  dim = isl_map_get_space(map);
2059
3
  diff = equate_parameter_to_length(dim, param);
2060
3
  map = isl_map_intersect(map, diff);
2061
3
  map = isl_map_project_out(map, isl_dim_in, d, 1);
2062
3
  map = isl_map_project_out(map, isl_dim_out, d, 1);
2063
3
  map = isl_map_from_range(isl_map_wrap(map));
2064
3
  map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1);
2065
3
2066
3
  map = isl_map_reset_space(map, target_dim);
2067
3
2068
3
  return map;
2069
3
}
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
0
{
2083
0
  isl_space *dim;
2084
0
  isl_map *diff;
2085
0
  unsigned d;
2086
0
  unsigned param;
2087
0
2088
0
  if (!map)
2089
0
    return NULL;
2090
0
2091
0
  d = isl_map_dim(map, isl_dim_in);
2092
0
  param = isl_map_dim(map, isl_dim_param);
2093
0
2094
0
  map = isl_map_compute_divs(map);
2095
0
  map = isl_map_coalesce(map);
2096
0
2097
0
  if (isl_map_plain_is_empty(map)) {
2098
0
    if (exact)
2099
0
      *exact = 1;
2100
0
    map = isl_map_project_out(map, isl_dim_out, 0, d);
2101
0
    map = isl_map_add_dims(map, isl_dim_out, 1);
2102
0
    return map;
2103
0
  }
2104
0
2105
0
  map = map_power(map, exact, 0);
2106
0
2107
0
  map = isl_map_add_dims(map, isl_dim_param, 1);
2108
0
  dim = isl_map_get_space(map);
2109
0
  diff = equate_parameter_to_length(dim, param);
2110
0
  map = isl_map_intersect(map, diff);
2111
0
  map = isl_map_project_out(map, isl_dim_in, 0, d + 1);
2112
0
  map = isl_map_project_out(map, isl_dim_out, d, 1);
2113
0
  map = isl_map_reverse(map);
2114
0
  map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1);
2115
0
2116
0
  return map;
2117
0
}
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
 */
2148
static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,
2149
  __isl_take isl_set *dom, __isl_take isl_set *ran, int with_id)
2150
0
{
2151
0
  int i;
2152
0
  int k;
2153
0
  unsigned d;
2154
0
  unsigned nparam;
2155
0
  unsigned total;
2156
0
  isl_space *dim;
2157
0
  isl_set *delta;
2158
0
  isl_map *app = NULL;
2159
0
  isl_basic_set *aff = NULL;
2160
0
  isl_basic_map *bmap = NULL;
2161
0
  isl_vec *obj = NULL;
2162
0
  isl_int opt;
2163
0
2164
0
  isl_int_init(opt);
2165
0
2166
0
  delta = isl_map_deltas(isl_map_copy(map));
2167
0
2168
0
  aff = isl_set_affine_hull(isl_set_copy(delta));
2169
0
  if (!aff)
2170
0
    goto error;
2171
0
  dim = isl_map_get_space(map);
2172
0
  d = isl_space_dim(dim, isl_dim_in);
2173
0
  nparam = isl_space_dim(dim, isl_dim_param);
2174
0
  total = isl_space_dim(dim, isl_dim_all);
2175
0
  bmap = isl_basic_map_alloc_space(dim,
2176
0
          aff->n_div + 1, aff->n_div, 2 * d + 1);
2177
0
  for (i = 0; i < aff->n_div + 1; ++i) {
2178
0
    k = isl_basic_map_alloc_div(bmap);
2179
0
    if (k < 0)
2180
0
      goto error;
2181
0
    isl_int_set_si(bmap->div[k][0], 0);
2182
0
  }
2183
0
  for (i = 0; i < aff->n_eq; ++i) {
2184
0
    if (!isl_basic_set_eq_is_stride(aff, i))
2185
0
      continue;
2186
0
    k = isl_basic_map_alloc_equality(bmap);
2187
0
    if (k < 0)
2188
0
      goto error;
2189
0
    isl_seq_clr(bmap->eq[k], 1 + nparam);
2190
0
    isl_seq_cpy(bmap->eq[k] + 1 + nparam + d,
2191
0
        aff->eq[i] + 1 + nparam, d);
2192
0
    isl_seq_neg(bmap->eq[k] + 1 + nparam,
2193
0
        aff->eq[i] + 1 + nparam, d);
2194
0
    isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d,
2195
0
        aff->eq[i] + 1 + nparam + d, aff->n_div);
2196
0
    isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);
2197
0
  }
2198
0
  obj = isl_vec_alloc(map->ctx, 1 + nparam + d);
2199
0
  if (!obj)
2200
0
    goto error;
2201
0
  isl_seq_clr(obj->el, 1 + nparam + d);
2202
0
  for (i = 0; i < d; ++ i) {
2203
0
    enum isl_lp_result res;
2204
0
2205
0
    isl_int_set_si(obj->el[1 + nparam + i], 1);
2206
0
2207
0
    res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt,
2208
0
          NULL, NULL);
2209
0
    if (res == isl_lp_error)
2210
0
      goto error;
2211
0
    if (res == isl_lp_ok) {
2212
0
      k = isl_basic_map_alloc_inequality(bmap);
2213
0
      if (k < 0)
2214
0
        goto error;
2215
0
      isl_seq_clr(bmap->ineq[k],
2216
0
          1 + nparam + 2 * d + bmap->n_div);
2217
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1);
2218
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1);
2219
0
      isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2220
0
    }
2221
0
2222
0
    res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt,
2223
0
          NULL, NULL);
2224
0
    if (res == isl_lp_error)
2225
0
      goto error;
2226
0
    if (res == isl_lp_ok) {
2227
0
      k = isl_basic_map_alloc_inequality(bmap);
2228
0
      if (k < 0)
2229
0
        goto error;
2230
0
      isl_seq_clr(bmap->ineq[k],
2231
0
          1 + nparam + 2 * d + bmap->n_div);
2232
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1);
2233
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1);
2234
0
      isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2235
0
    }
2236
0
2237
0
    isl_int_set_si(obj->el[1 + nparam + i], 0);
2238
0
  }
2239
0
  k = isl_basic_map_alloc_inequality(bmap);
2240
0
  if (k < 0)
2241
0
    goto error;
2242
0
  isl_seq_clr(bmap->ineq[k],
2243
0
      1 + nparam + 2 * d + bmap->n_div);
2244
0
  if (!with_id)
2245
0
    isl_int_set_si(bmap->ineq[k][0], -1);
2246
0
  isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1);
2247
0
2248
0
  app = isl_map_from_domain_and_range(dom, ran);
2249
0
2250
0
  isl_vec_free(obj);
2251
0
  isl_basic_set_free(aff);
2252
0
  isl_map_free(map);
2253
0
  bmap = isl_basic_map_finalize(bmap);
2254
0
  isl_set_free(delta);
2255
0
  isl_int_clear(opt);
2256
0
2257
0
  map = isl_map_from_basic_map(bmap);
2258
0
  map = isl_map_intersect(map, app);
2259
0
2260
0
  return map;
2261
0
error:
2262
0
  isl_vec_free(obj);
2263
0
  isl_basic_map_free(bmap);
2264
0
  isl_basic_set_free(aff);
2265
0
  isl_set_free(dom);
2266
0
  isl_set_free(ran);
2267
0
  isl_map_free(map);
2268
0
  isl_set_free(delta);
2269
0
  isl_int_clear(opt);
2270
0
  return NULL;
2271
0
}
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
 */
2286
static __isl_give isl_map *box_closure(__isl_take isl_map *map)
2287
0
{
2288
0
  isl_set *domain;
2289
0
  isl_set *range;
2290
0
2291
0
  domain = isl_map_domain(isl_map_copy(map));
2292
0
  domain = isl_set_coalesce(domain);
2293
0
  range = isl_map_range(isl_map_copy(map));
2294
0
  range = isl_set_coalesce(range);
2295
0
2296
0
  return box_closure_on_domain(map, domain, range, 0);
2297
0
}
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
 */
2312
static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map,
2313
  __isl_take isl_set *dom)
2314
0
{
2315
0
  return box_closure_on_domain(map, dom, isl_set_copy(dom), 1);
2316
0
}
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
 */
2324
static int check_exactness_omega(__isl_keep isl_map *map,
2325
  __isl_keep isl_map *app)
2326
0
{
2327
0
  isl_set *delta;
2328
0
  int i;
2329
0
  int is_empty, is_exact;
2330
0
  unsigned d;
2331
0
  isl_map *test;
2332
0
2333
0
  delta = isl_map_deltas(isl_map_copy(app));
2334
0
  d = isl_set_dim(delta, isl_dim_set);
2335
0
  for (i = 0; i < d; ++i)
2336
0
    delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
2337
0
  is_empty = isl_set_is_empty(delta);
2338
0
  isl_set_free(delta);
2339
0
  if (is_empty < 0)
2340
0
    return -1;
2341
0
  if (!is_empty)
2342
0
    return 0;
2343
0
2344
0
  test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map));
2345
0
  test = isl_map_union(test, isl_map_copy(map));
2346
0
  is_exact = isl_map_is_subset(app, test);
2347
0
  isl_map_free(test);
2348
0
2349
0
  return is_exact;
2350
0
}
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
 */
2380
static int can_be_split_off(__isl_keep isl_map *map, int i,
2381
  __isl_give isl_map **tc, __isl_give isl_map **qc)
2382
0
{
2383
0
  isl_map *map_i, *id = NULL;
2384
0
  int j = -1;
2385
0
  isl_set *C;
2386
0
2387
0
  *tc = NULL;
2388
0
  *qc = NULL;
2389
0
2390
0
  C = isl_set_union(isl_map_domain(isl_map_copy(map)),
2391
0
        isl_map_range(isl_map_copy(map)));
2392
0
  C = isl_set_from_basic_set(isl_set_simple_hull(C));
2393
0
  if (!C)
2394
0
    goto error;
2395
0
2396
0
  map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
2397
0
  *tc = box_closure(isl_map_copy(map_i));
2398
0
  *qc = box_closure_with_identity(map_i, C);
2399
0
  id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc));
2400
0
2401
0
  if (!id || !*qc)
2402
0
    goto error;
2403
0
  if (id->n != 1 || (*qc)->n != 1)
2404
0
    goto done;
2405
0
2406
0
  for (j = 0; j < map->n; ++j) {
2407
0
    isl_map *map_j, *test;
2408
0
    int is_ok;
2409
0
2410
0
    if (i == j)
2411
0
      continue;
2412
0
    map_j = isl_map_from_basic_map(
2413
0
          isl_basic_map_copy(map->p[j]));
2414
0
    test = isl_map_apply_range(isl_map_copy(id),
2415
0
            isl_map_copy(map_j));
2416
0
    test = isl_map_apply_range(test, isl_map_copy(id));
2417
0
    is_ok = isl_map_is_equal(test, map_j);
2418
0
    isl_map_free(map_j);
2419
0
    isl_map_free(test);
2420
0
    if (is_ok < 0)
2421
0
      goto error;
2422
0
    if (!is_ok)
2423
0
      break;
2424
0
  }
2425
0
2426
0
done:
2427
0
  isl_map_free(id);
2428
0
  if (j == map->n)
2429
0
    return 1;
2430
0
2431
0
  isl_map_free(*qc);
2432
0
  isl_map_free(*tc);
2433
0
  *qc = NULL;
2434
0
  *tc = NULL;
2435
0
2436
0
  return 0;
2437
0
error:
2438
0
  isl_map_free(id);
2439
0
  isl_map_free(*qc);
2440
0
  isl_map_free(*tc);
2441
0
  *qc = NULL;
2442
0
  *tc = NULL;
2443
0
  return -1;
2444
0
}
2445
2446
static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map,
2447
  int *exact)
2448
0
{
2449
0
  isl_map *app;
2450
0
2451
0
  app = box_closure(isl_map_copy(map));
2452
0
  if (exact)
2453
0
    *exact = check_exactness_omega(map, app);
2454
0
2455
0
  isl_map_free(map);
2456
0
  return app;
2457
0
}
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
 */
2477
static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,
2478
  int *exact)
2479
0
{
2480
0
  int i, j;
2481
0
  int exact_i;
2482
0
  isl_map *app;
2483
0
2484
0
  if (!map)
2485
0
    return NULL;
2486
0
  if (map->n == 1)
2487
0
    return box_closure_with_check(map, exact);
2488
0
2489
0
  for (i = 0; i < map->n; ++i) {
2490
0
    int ok;
2491
0
    isl_map *qc, *tc;
2492
0
    ok = can_be_split_off(map, i, &tc, &qc);
2493
0
    if (ok < 0)
2494
0
      goto error;
2495
0
    if (!ok)
2496
0
      continue;
2497
0
2498
0
    app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0);
2499
0
2500
0
    for (j = 0; j < map->n; ++j) {
2501
0
      if (j == i)
2502
0
        continue;
2503
0
      app = isl_map_add_basic_map(app,
2504
0
            isl_basic_map_copy(map->p[j]));
2505
0
    }
2506
0
2507
0
    app = isl_map_apply_range(isl_map_copy(qc), app);
2508
0
    app = isl_map_apply_range(app, qc);
2509
0
2510
0
    app = isl_map_union(tc, transitive_closure_omega(app, NULL));
2511
0
    exact_i = check_exactness_omega(map, app);
2512
0
    if (exact_i == 1) {
2513
0
      if (exact)
2514
0
        *exact = exact_i;
2515
0
      isl_map_free(map);
2516
0
      return app;
2517
0
    }
2518
0
    isl_map_free(app);
2519
0
    if (exact_i < 0)
2520
0
      goto error;
2521
0
  }
2522
0
2523
0
  return box_closure_with_check(map, exact);
2524
0
error:
2525
0
  isl_map_free(map);
2526
0
  return NULL;
2527
0
}
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
115
{
2538
115
  isl_space *target_dim;
2539
115
  int closed;
2540
115
2541
115
  if (!map)
2542
0
    goto error;
2543
115
2544
115
  if (map->ctx->opt->closure == ISL_CLOSURE_BOX)
2545
115
    
return transitive_closure_omega(map, exact)0
;
2546
115
2547
115
  map = isl_map_compute_divs(map);
2548
115
  map = isl_map_coalesce(map);
2549
115
  closed = isl_map_is_transitively_closed(map);
2550
115
  if (closed < 0)
2551
0
    goto error;
2552
115
  if (closed) {
2553
8
    if (exact)
2554
6
      *exact = 1;
2555
8
    return map;
2556
8
  }
2557
107
2558
107
  target_dim = isl_map_get_space(map);
2559
107
  map = map_power(map, exact, 1);
2560
107
  map = isl_map_reset_space(map, target_dim);
2561
107
2562
107
  return map;
2563
0
error:
2564
0
  isl_map_free(map);
2565
0
  return NULL;
2566
115
}
2567
2568
static isl_stat inc_count(__isl_take isl_map *map, void *user)
2569
183
{
2570
183
  int *n = user;
2571
183
2572
183
  *n += map->n;
2573
183
2574
183
  isl_map_free(map);
2575
183
2576
183
  return isl_stat_ok;
2577
183
}
2578
2579
static isl_stat collect_basic_map(__isl_take isl_map *map, void *user)
2580
126
{
2581
126
  int i;
2582
126
  isl_basic_map ***next = user;
2583
126
2584
293
  for (i = 0; i < map->n; 
++i167
) {
2585
167
    **next = isl_basic_map_copy(map->p[i]);
2586
167
    if (!**next)
2587
0
      goto error;
2588
167
    (*next)++;
2589
167
  }
2590
126
2591
126
  isl_map_free(map);
2592
126
  return isl_stat_ok;
2593
0
error:
2594
0
  isl_map_free(map);
2595
0
  return isl_stat_error;
2596
126
}
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
 */
2606
static __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
93
{
2609
93
  int i, j, k;
2610
93
  int n_group;
2611
93
  int *group = NULL;
2612
93
  isl_set **set = NULL;
2613
93
  isl_map ***grid = NULL;
2614
93
  isl_union_map *app;
2615
93
2616
93
  group = setup_groups(ctx, list, n, &set, &n_group);
2617
93
  if (!group)
2618
0
    goto error;
2619
93
2620
93
  grid = isl_calloc_array(ctx, isl_map **, n_group);
2621
93
  if (!grid)
2622
0
    goto error;
2623
187
  
for (i = 0; 93
i < n_group;
++i94
) {
2624
94
    grid[i] = isl_calloc_array(ctx, isl_map *, n_group);
2625
94
    if (!grid[i])
2626
0
      goto error;
2627
190
    
for (j = 0; 94
j < n_group;
++j96
) {
2628
96
      isl_space *dim1, *dim2, *dim;
2629
96
      dim1 = isl_space_reverse(isl_set_get_space(set[i]));
2630
96
      dim2 = isl_set_get_space(set[j]);
2631
96
      dim = isl_space_join(dim1, dim2);
2632
96
      grid[i][j] = isl_map_empty(dim);
2633
96
    }
2634
94
  }
2635
93
2636
205
  
for (k = 0; 93
k < n;
++k112
) {
2637
112
    i = group[2 * k];
2638
112
    j = group[2 * k + 1];
2639
112
    grid[i][j] = isl_map_union(grid[i][j],
2640
112
        isl_map_from_basic_map(
2641
112
          isl_basic_map_copy(list[k])));
2642
112
  }
2643
93
  
2644
93
  floyd_warshall_iterate(grid, n_group, exact);
2645
93
2646
93
  app = isl_union_map_empty(isl_map_get_space(grid[0][0]));
2647
93
2648
187
  for (i = 0; i < n_group; 
++i94
) {
2649
190
    for (j = 0; j < n_group; 
++j96
)
2650
96
      app = isl_union_map_add_map(app, grid[i][j]);
2651
94
    free(grid[i]);
2652
94
  }
2653
93
  free(grid);
2654
93
2655
317
  for (i = 0; i < 2 * n; 
++i224
)
2656
224
    isl_set_free(set[i]);
2657
93
  free(set);
2658
93
2659
93
  free(group);
2660
93
  return app;
2661
0
error:
2662
0
  if (grid)
2663
0
    for (i = 0; i < n_group; ++i) {
2664
0
      if (!grid[i])
2665
0
        continue;
2666
0
      for (j = 0; j < n_group; ++j)
2667
0
        isl_map_free(grid[i][j]);
2668
0
      free(grid[i]);
2669
0
    }
2670
0
  free(grid);
2671
0
  if (set) {
2672
0
    for (i = 0; i < 2 * n; ++i)
2673
0
      isl_set_free(set[i]);
2674
0
    free(set);
2675
0
  }
2676
0
  free(group);
2677
0
  return NULL;
2678
93
}
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
 */
2686
static __isl_give isl_union_map *union_floyd_warshall(
2687
  __isl_take isl_union_map *umap, int *exact)
2688
93
{
2689
93
  int i, n;
2690
93
  isl_ctx *ctx;
2691
93
  isl_basic_map **list = NULL;
2692
93
  isl_basic_map **next;
2693
93
  isl_union_map *res;
2694
93
2695
93
  n = 0;
2696
93
  if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
2697
0
    goto error;
2698
93
2699
93
  ctx = isl_union_map_get_ctx(umap);
2700
93
  list = isl_calloc_array(ctx, isl_basic_map *, n);
2701
93
  if (!list)
2702
0
    goto error;
2703
93
2704
93
  next = list;
2705
93
  if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
2706
0
    goto error;
2707
93
2708
93
  res = union_floyd_warshall_on_list(ctx, list, n, exact);
2709
93
2710
93
  if (list) {
2711
205
    for (i = 0; i < n; 
++i112
)
2712
112
      isl_basic_map_free(list[i]);
2713
93
    free(list);
2714
93
  }
2715
93
2716
93
  isl_union_map_free(umap);
2717
93
  return res;
2718
0
error:
2719
0
  if (list) {
2720
0
    for (i = 0; i < n; ++i)
2721
0
      isl_basic_map_free(list[i]);
2722
0
    free(list);
2723
0
  }
2724
0
  isl_union_map_free(umap);
2725
0
  return NULL;
2726
93
}
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
 */
2733
static __isl_give isl_union_map *union_components(
2734
  __isl_take isl_union_map *umap, int *exact)
2735
70
{
2736
70
  int i;
2737
70
  int n;
2738
70
  isl_ctx *ctx;
2739
70
  isl_basic_map **list = NULL;
2740
70
  isl_basic_map **next;
2741
70
  isl_union_map *path = NULL;
2742
70
  struct isl_tc_follows_data data;
2743
70
  struct isl_tarjan_graph *g = NULL;
2744
70
  int c, l;
2745
70
  int recheck = 0;
2746
70
2747
70
  n = 0;
2748
70
  if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
2749
0
    goto error;
2750
70
2751
70
  if (n == 0)
2752
0
    return umap;
2753
70
  if (n <= 1)
2754
57
    return union_floyd_warshall(umap, exact);
2755
13
2756
13
  ctx = isl_union_map_get_ctx(umap);
2757
13
  list = isl_calloc_array(ctx, isl_basic_map *, n);
2758
13
  if (!list)
2759
0
    goto error;
2760
13
2761
13
  next = list;
2762
13
  if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
2763
0
    goto error;
2764
13
2765
13
  data.list = list;
2766
13
  data.check_closed = 0;
2767
13
  g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data);
2768
13
  if (!g)
2769
0
    goto error;
2770
13
2771
13
  c = 0;
2772
13
  i = 0;
2773
13
  l = n;
2774
13
  path = isl_union_map_empty(isl_union_map_get_space(umap));
2775
49
  while (l) {
2776
36
    isl_union_map *comp;
2777
36
    isl_union_map *path_comp, *path_comb;
2778
36
    comp = isl_union_map_empty(isl_union_map_get_space(umap));
2779
91
    while (g->order[i] != -1) {
2780
55
      comp = isl_union_map_add_map(comp,
2781
55
            isl_map_from_basic_map(
2782
55
          isl_basic_map_copy(list[g->order[i]])));
2783
55
      --l;
2784
55
      ++i;
2785
55
    }
2786
36
    path_comp = union_floyd_warshall(comp, exact);
2787
36
    path_comb = isl_union_map_apply_range(isl_union_map_copy(path),
2788
36
            isl_union_map_copy(path_comp));
2789
36
    path = isl_union_map_union(path, path_comp);
2790
36
    path = isl_union_map_union(path, path_comb);
2791
36
    ++i;
2792
36
    ++c;
2793
36
  }
2794
13
2795
13
  if (c > 1 && 
data.check_closed8
&&
!*exact0
) {
2796
0
    int closed;
2797
0
2798
0
    closed = isl_union_map_is_transitively_closed(path);
2799
0
    if (closed < 0)
2800
0
      goto error;
2801
0
    recheck = !closed;
2802
0
  }
2803
13
2804
13
  isl_tarjan_graph_free(g);
2805
13
2806
68
  for (i = 0; i < n; 
++i55
)
2807
55
    isl_basic_map_free(list[i]);
2808
13
  free(list);
2809
13
2810
13
  if (recheck) {
2811
0
    isl_union_map_free(path);
2812
0
    return union_floyd_warshall(umap, exact);
2813
0
  }
2814
13
2815
13
  isl_union_map_free(umap);
2816
13
2817
13
  return path;
2818
0
error:
2819
0
  isl_tarjan_graph_free(g);
2820
0
  if (list) {
2821
0
    for (i = 0; i < n; ++i)
2822
0
      isl_basic_map_free(list[i]);
2823
0
    free(list);
2824
0
  }
2825
0
  isl_union_map_free(umap);
2826
0
  isl_union_map_free(path);
2827
0
  return NULL;
2828
70
}
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
72
{
2836
72
  int closed;
2837
72
2838
72
  if (!umap)
2839
0
    return NULL;
2840
72
2841
72
  if (exact)
2842
0
    *exact = 1;
2843
72
2844
72
  umap = isl_union_map_compute_divs(umap);
2845
72
  umap = isl_union_map_coalesce(umap);
2846
72
  closed = isl_union_map_is_transitively_closed(umap);
2847
72
  if (closed < 0)
2848
0
    goto error;
2849
72
  if (closed)
2850
2
    return umap;
2851
70
  umap = union_components(umap, exact);
2852
70
  return umap;
2853
0
error:
2854
0
  isl_union_map_free(umap);
2855
0
  return NULL;
2856
72
}
2857
2858
struct isl_union_power {
2859
  isl_union_map *pow;
2860
  int *exact;
2861
};
2862
2863
static isl_stat power(__isl_take isl_map *map, void *user)
2864
0
{
2865
0
  struct isl_union_power *up = user;
2866
0
2867
0
  map = isl_map_power(map, up->exact);
2868
0
  up->pow = isl_union_map_from_map(map);
2869
0
2870
0
  return isl_stat_error;
2871
0
}
2872
2873
/* Construct a map [x] -> [x+1], with parameters prescribed by "dim".
2874
 */
2875
static __isl_give isl_union_map *increment(__isl_take isl_space *dim)
2876
0
{
2877
0
  int k;
2878
0
  isl_basic_map *bmap;
2879
0
2880
0
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
2881
0
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
2882
0
  bmap = isl_basic_map_alloc_space(dim, 0, 1, 0);
2883
0
  k = isl_basic_map_alloc_equality(bmap);
2884
0
  if (k < 0)
2885
0
    goto error;
2886
0
  isl_seq_clr(bmap->eq[k], isl_basic_map_total_dim(bmap));
2887
0
  isl_int_set_si(bmap->eq[k][0], 1);
2888
0
  isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);
2889
0
  isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);
2890
0
  return isl_union_map_from_map(isl_map_from_basic_map(bmap));
2891
0
error:
2892
0
  isl_basic_map_free(bmap);
2893
0
  return NULL;
2894
0
}
2895
2896
/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim".
2897
 */
2898
static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim)
2899
0
{
2900
0
  isl_basic_map *bmap;
2901
0
2902
0
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
2903
0
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
2904
0
  bmap = isl_basic_map_universe(dim);
2905
0
  bmap = isl_basic_map_deltas_map(bmap);
2906
0
2907
0
  return isl_union_map_from_map(isl_map_from_basic_map(bmap));
2908
0
}
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
0
{
2917
0
  int n;
2918
0
  isl_union_map *inc;
2919
0
  isl_union_map *dm;
2920
0
2921
0
  if (!umap)
2922
0
    return NULL;
2923
0
  n = isl_union_map_n_map(umap);
2924
0
  if (n == 0)
2925
0
    return umap;
2926
0
  if (n == 1) {
2927
0
    struct isl_union_power up = { NULL, exact };
2928
0
    isl_union_map_foreach_map(umap, &power, &up);
2929
0
    isl_union_map_free(umap);
2930
0
    return up.pow;
2931
0
  }
2932
0
  inc = increment(isl_union_map_get_space(umap));
2933
0
  umap = isl_union_map_product(inc, umap);
2934
0
  umap = isl_union_map_transitive_closure(umap, exact);
2935
0
  umap = isl_union_map_zip(umap);
2936
0
  dm = deltas_map(isl_union_map_get_space(umap));
2937
0
  umap = isl_union_map_apply_domain(umap, dm);
2938
0
  
2939
0
  return umap;
2940
0
}
2941
2942
#undef TYPE
2943
1
#define TYPE isl_map
2944
#include "isl_power_templ.c"
2945
2946
#undef TYPE
2947
0
#define TYPE isl_union_map
2948
#include "isl_power_templ.c"