Coverage Report

Created: 2017-10-03 07:32

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/isl_convex_hull.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2014      INRIA Rocquencourt
4
 *
5
 * Use of this software is governed by the MIT license
6
 *
7
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
8
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9
 * and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
10
 * B.P. 105 - 78153 Le Chesnay, France
11
 */
12
13
#include <isl_ctx_private.h>
14
#include <isl_map_private.h>
15
#include <isl_lp_private.h>
16
#include <isl/map.h>
17
#include <isl_mat_private.h>
18
#include <isl_vec_private.h>
19
#include <isl/set.h>
20
#include <isl_seq.h>
21
#include <isl_options_private.h>
22
#include "isl_equalities.h"
23
#include "isl_tab.h"
24
#include <isl_sort.h>
25
26
#include <bset_to_bmap.c>
27
#include <bset_from_bmap.c>
28
#include <set_to_map.c>
29
30
static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(
31
  __isl_take isl_set *set);
32
33
/* Remove redundant
34
 * constraints.  If the minimal value along the normal of a constraint
35
 * is the same if the constraint is removed, then the constraint is redundant.
36
 *
37
 * Since some constraints may be mutually redundant, sort the constraints
38
 * first such that constraints that involve existentially quantified
39
 * variables are considered for removal before those that do not.
40
 * The sorting is also needed for the use in map_simple_hull.
41
 *
42
 * Note that isl_tab_detect_implicit_equalities may also end up
43
 * marking some constraints as redundant.  Make sure the constraints
44
 * are preserved and undo those marking such that isl_tab_detect_redundant
45
 * can consider the constraints in the sorted order.
46
 *
47
 * Alternatively, we could have intersected the basic map with the
48
 * corresponding equality and then checked if the dimension was that
49
 * of a facet.
50
 */
51
__isl_give isl_basic_map *isl_basic_map_remove_redundancies(
52
  __isl_take isl_basic_map *bmap)
53
300k
{
54
300k
  struct isl_tab *tab;
55
300k
56
300k
  if (!bmap)
57
0
    return NULL;
58
300k
59
300k
  bmap = isl_basic_map_gauss(bmap, NULL);
60
300k
  if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
61
1.16k
    return bmap;
62
299k
  
if (299k
ISL_F_ISSET299k
(bmap, ISL_BASIC_MAP_NO_REDUNDANT))
63
93.7k
    return bmap;
64
205k
  
if (205k
bmap->n_ineq <= 1205k
)
65
95.2k
    return bmap;
66
110k
67
110k
  bmap = isl_basic_map_sort_constraints(bmap);
68
110k
  tab = isl_tab_from_basic_map(bmap, 0);
69
110k
  if (!tab)
70
0
    goto error;
71
110k
  tab->preserve = 1;
72
110k
  if (isl_tab_detect_implicit_equalities(tab) < 0)
73
0
    goto error;
74
110k
  
if (110k
isl_tab_restore_redundant(tab) < 0110k
)
75
0
    goto error;
76
110k
  tab->preserve = 0;
77
110k
  if (isl_tab_detect_redundant(tab) < 0)
78
0
    goto error;
79
110k
  bmap = isl_basic_map_update_from_tab(bmap, tab);
80
110k
  isl_tab_free(tab);
81
110k
  if (!bmap)
82
0
    return NULL;
83
110k
  
ISL_F_SET110k
(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
84
110k
  ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
85
110k
  return bmap;
86
0
error:
87
0
  isl_tab_free(tab);
88
0
  isl_basic_map_free(bmap);
89
0
  return NULL;
90
300k
}
91
92
__isl_give isl_basic_set *isl_basic_set_remove_redundancies(
93
  __isl_take isl_basic_set *bset)
94
3.86k
{
95
3.86k
  return bset_from_bmap(
96
3.86k
    isl_basic_map_remove_redundancies(bset_to_bmap(bset)));
97
3.86k
}
98
99
/* Remove redundant constraints in each of the basic maps.
100
 */
101
__isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map)
102
1.06k
{
103
1.06k
  return isl_map_inline_foreach_basic_map(map,
104
1.06k
              &isl_basic_map_remove_redundancies);
105
1.06k
}
106
107
__isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set)
108
1.04k
{
109
1.04k
  return isl_map_remove_redundancies(set);
110
1.04k
}
111
112
/* Check if the set set is bound in the direction of the affine
113
 * constraint c and if so, set the constant term such that the
114
 * resulting constraint is a bounding constraint for the set.
115
 */
116
static int uset_is_bound(__isl_keep isl_set *set, isl_int *c, unsigned len)
117
3.88k
{
118
3.88k
  int first;
119
3.88k
  int j;
120
3.88k
  isl_int opt;
121
3.88k
  isl_int opt_denom;
122
3.88k
123
3.88k
  isl_int_init(opt);
124
3.88k
  isl_int_init(opt_denom);
125
3.88k
  first = 1;
126
11.6k
  for (j = 0; 
j < set->n11.6k
;
++j7.76k
) {
127
7.76k
    enum isl_lp_result res;
128
7.76k
129
7.76k
    if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
130
0
      continue;
131
7.76k
132
7.76k
    res = isl_basic_set_solve_lp(set->p[j],
133
7.76k
        0, c, set->ctx->one, &opt, &opt_denom, NULL);
134
7.76k
    if (res == isl_lp_unbounded)
135
0
      break;
136
7.76k
    
if (7.76k
res == isl_lp_error7.76k
)
137
0
      goto error;
138
7.76k
    
if (7.76k
res == isl_lp_empty7.76k
) {
139
0
      set->p[j] = isl_basic_set_set_to_empty(set->p[j]);
140
0
      if (!set->p[j])
141
0
        goto error;
142
0
      continue;
143
0
    }
144
7.76k
    
if (7.76k
first || 7.76k
isl_int_is_neg3.88k
(opt)) {
145
4.75k
      if (
!4.75k
isl_int_is_one4.75k
(opt_denom))
146
4.54k
        isl_seq_scale(c, c, opt_denom, len);
147
4.75k
      isl_int_sub(c[0], c[0], opt);
148
4.75k
    }
149
7.76k
    first = 0;
150
7.76k
  }
151
3.88k
  
isl_int_clear3.88k
(opt);
152
3.88k
  isl_int_clear(opt_denom);
153
3.88k
  return j >= set->n;
154
0
error:
155
0
  isl_int_clear(opt);
156
0
  isl_int_clear(opt_denom);
157
0
  return -1;
158
3.88k
}
159
160
static struct isl_basic_set *isl_basic_set_add_equality(
161
  struct isl_basic_set *bset, isl_int *c)
162
42.3k
{
163
42.3k
  int i;
164
42.3k
  unsigned dim;
165
42.3k
166
42.3k
  if (!bset)
167
0
    return NULL;
168
42.3k
169
42.3k
  
if (42.3k
ISL_F_ISSET42.3k
(bset, ISL_BASIC_SET_EMPTY))
170
0
    return bset;
171
42.3k
172
42.3k
  
isl_assert42.3k
(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
173
42.3k
  
isl_assert42.3k
(bset->ctx, bset->n_div == 0, goto error);
174
42.3k
  dim = isl_basic_set_n_dim(bset);
175
42.3k
  bset = isl_basic_set_cow(bset);
176
42.3k
  bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0);
177
42.3k
  i = isl_basic_set_alloc_equality(bset);
178
42.3k
  if (i < 0)
179
0
    goto error;
180
42.3k
  isl_seq_cpy(bset->eq[i], c, 1 + dim);
181
42.3k
  return bset;
182
0
error:
183
0
  isl_basic_set_free(bset);
184
0
  return NULL;
185
42.3k
}
186
187
static __isl_give isl_set *isl_set_add_basic_set_equality(
188
  __isl_take isl_set *set, isl_int *c)
189
7.59k
{
190
7.59k
  int i;
191
7.59k
192
7.59k
  set = isl_set_cow(set);
193
7.59k
  if (!set)
194
0
    return NULL;
195
22.7k
  
for (i = 0; 7.59k
i < set->n22.7k
;
++i15.1k
) {
196
15.1k
    set->p[i] = isl_basic_set_add_equality(set->p[i], c);
197
15.1k
    if (!set->p[i])
198
0
      goto error;
199
15.1k
  }
200
7.59k
  return set;
201
0
error:
202
0
  isl_set_free(set);
203
0
  return NULL;
204
7.59k
}
205
206
/* Given a union of basic sets, construct the constraints for wrapping
207
 * a facet around one of its ridges.
208
 * In particular, if each of n the d-dimensional basic sets i in "set"
209
 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0
210
 * and is defined by the constraints
211
 *            [ 1 ]
212
 *        A_i [ x ]  >= 0
213
 *
214
 * then the resulting set is of dimension n*(1+d) and has as constraints
215
 *
216
 *            [ a_i ]
217
 *        A_i [ x_i ] >= 0
218
 *
219
 *              a_i   >= 0
220
 *
221
 *      \sum_i x_{i,1} = 1
222
 */
223
static __isl_give isl_basic_set *wrap_constraints(__isl_keep isl_set *set)
224
18.5k
{
225
18.5k
  struct isl_basic_set *lp;
226
18.5k
  unsigned n_eq;
227
18.5k
  unsigned n_ineq;
228
18.5k
  int i, j, k;
229
18.5k
  unsigned dim, lp_dim;
230
18.5k
231
18.5k
  if (!set)
232
0
    return NULL;
233
18.5k
234
18.5k
  dim = 1 + isl_set_n_dim(set);
235
18.5k
  n_eq = 1;
236
18.5k
  n_ineq = set->n;
237
50.5k
  for (i = 0; 
i < set->n50.5k
;
++i31.9k
) {
238
31.9k
    n_eq += set->p[i]->n_eq;
239
31.9k
    n_ineq += set->p[i]->n_ineq;
240
31.9k
  }
241
18.5k
  lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq);
242
18.5k
  lp = isl_basic_set_set_rational(lp);
243
18.5k
  if (!lp)
244
0
    return NULL;
245
18.5k
  lp_dim = isl_basic_set_n_dim(lp);
246
18.5k
  k = isl_basic_set_alloc_equality(lp);
247
18.5k
  isl_int_set_si(lp->eq[k][0], -1);
248
50.5k
  for (i = 0; 
i < set->n50.5k
;
++i31.9k
) {
249
31.9k
    isl_int_set_si(lp->eq[k][1+dim*i], 0);
250
31.9k
    isl_int_set_si(lp->eq[k][1+dim*i+1], 1);
251
31.9k
    isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2);
252
31.9k
  }
253
50.5k
  for (i = 0; 
i < set->n50.5k
;
++i31.9k
) {
254
31.9k
    k = isl_basic_set_alloc_inequality(lp);
255
31.9k
    isl_seq_clr(lp->ineq[k], 1+lp_dim);
256
31.9k
    isl_int_set_si(lp->ineq[k][1+dim*i], 1);
257
31.9k
258
76.0k
    for (j = 0; 
j < set->p[i]->n_eq76.0k
;
++j44.1k
) {
259
44.1k
      k = isl_basic_set_alloc_equality(lp);
260
44.1k
      isl_seq_clr(lp->eq[k], 1+dim*i);
261
44.1k
      isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim);
262
44.1k
      isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1));
263
44.1k
    }
264
31.9k
265
111k
    for (j = 0; 
j < set->p[i]->n_ineq111k
;
++j79.9k
) {
266
79.9k
      k = isl_basic_set_alloc_inequality(lp);
267
79.9k
      isl_seq_clr(lp->ineq[k], 1+dim*i);
268
79.9k
      isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim);
269
79.9k
      isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1));
270
79.9k
    }
271
31.9k
  }
272
18.5k
  return lp;
273
18.5k
}
274
275
/* Given a facet "facet" of the convex hull of "set" and a facet "ridge"
276
 * of that facet, compute the other facet of the convex hull that contains
277
 * the ridge.
278
 *
279
 * We first transform the set such that the facet constraint becomes
280
 *
281
 *      x_1 >= 0
282
 *
283
 * I.e., the facet lies in
284
 *
285
 *      x_1 = 0
286
 *
287
 * and on that facet, the constraint that defines the ridge is
288
 *
289
 *      x_2 >= 0
290
 *
291
 * (This transformation is not strictly needed, all that is needed is
292
 * that the ridge contains the origin.)
293
 *
294
 * Since the ridge contains the origin, the cone of the convex hull
295
 * will be of the form
296
 *
297
 *      x_1 >= 0
298
 *      x_2 >= a x_1
299
 *
300
 * with this second constraint defining the new facet.
301
 * The constant a is obtained by settting x_1 in the cone of the
302
 * convex hull to 1 and minimizing x_2.
303
 * Now, each element in the cone of the convex hull is the sum
304
 * of elements in the cones of the basic sets.
305
 * If a_i is the dilation factor of basic set i, then the problem
306
 * we need to solve is
307
 *
308
 *      min \sum_i x_{i,2}
309
 *      st
310
 *        \sum_i x_{i,1} = 1
311
 *            a_i   >= 0
312
 *          [ a_i ]
313
 *        A [ x_i ] >= 0
314
 *
315
 * with
316
 *            [  1  ]
317
 *        A_i [ x_i ] >= 0
318
 *
319
 * the constraints of each (transformed) basic set.
320
 * If a = n/d, then the constraint defining the new facet (in the transformed
321
 * space) is
322
 *
323
 *      -n x_1 + d x_2 >= 0
324
 *
325
 * In the original space, we need to take the same combination of the
326
 * corresponding constraints "facet" and "ridge".
327
 *
328
 * If a = -infty = "-1/0", then we just return the original facet constraint.
329
 * This means that the facet is unbounded, but has a bounded intersection
330
 * with the union of sets.
331
 */
332
isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
333
  isl_int *facet, isl_int *ridge)
334
18.5k
{
335
18.5k
  int i;
336
18.5k
  isl_ctx *ctx;
337
18.5k
  struct isl_mat *T = NULL;
338
18.5k
  struct isl_basic_set *lp = NULL;
339
18.5k
  struct isl_vec *obj;
340
18.5k
  enum isl_lp_result res;
341
18.5k
  isl_int num, den;
342
18.5k
  unsigned dim;
343
18.5k
344
18.5k
  if (!set)
345
0
    return NULL;
346
18.5k
  ctx = set->ctx;
347
18.5k
  set = isl_set_copy(set);
348
18.5k
  set = isl_set_set_rational(set);
349
18.5k
350
18.5k
  dim = 1 + isl_set_n_dim(set);
351
18.5k
  T = isl_mat_alloc(ctx, 3, dim);
352
18.5k
  if (!T)
353
0
    goto error;
354
18.5k
  
isl_int_set_si18.5k
(T->row[0][0], 1);
355
18.5k
  isl_seq_clr(T->row[0]+1, dim - 1);
356
18.5k
  isl_seq_cpy(T->row[1], facet, dim);
357
18.5k
  isl_seq_cpy(T->row[2], ridge, dim);
358
18.5k
  T = isl_mat_right_inverse(T);
359
18.5k
  set = isl_set_preimage(set, T);
360
18.5k
  T = NULL;
361
18.5k
  if (!set)
362
0
    goto error;
363
18.5k
  lp = wrap_constraints(set);
364
18.5k
  obj = isl_vec_alloc(ctx, 1 + dim*set->n);
365
18.5k
  if (!obj)
366
0
    goto error;
367
18.5k
  
isl_int_set_si18.5k
(obj->block.data[0], 0);
368
50.5k
  for (i = 0; 
i < set->n50.5k
;
++i31.9k
) {
369
31.9k
    isl_seq_clr(obj->block.data + 1 + dim*i, 2);
370
31.9k
    isl_int_set_si(obj->block.data[1 + dim*i+2], 1);
371
31.9k
    isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3);
372
31.9k
  }
373
18.5k
  isl_int_init(num);
374
18.5k
  isl_int_init(den);
375
18.5k
  res = isl_basic_set_solve_lp(lp, 0,
376
18.5k
          obj->block.data, ctx->one, &num, &den, NULL);
377
18.5k
  if (
res == isl_lp_ok18.5k
) {
378
15.8k
    isl_int_neg(num, num);
379
15.8k
    isl_seq_combine(facet, num, facet, den, ridge, dim);
380
15.8k
    isl_seq_normalize(ctx, facet, dim);
381
15.8k
  }
382
18.5k
  isl_int_clear(num);
383
18.5k
  isl_int_clear(den);
384
18.5k
  isl_vec_free(obj);
385
18.5k
  isl_basic_set_free(lp);
386
18.5k
  isl_set_free(set);
387
18.5k
  if (res == isl_lp_error)
388
0
    return NULL;
389
18.5k
  
isl_assert18.5k
(ctx, res == isl_lp_ok || res == isl_lp_unbounded,
390
18.5k
       return NULL);
391
18.5k
  return facet;
392
0
error:
393
0
  isl_basic_set_free(lp);
394
0
  isl_mat_free(T);
395
0
  isl_set_free(set);
396
0
  return NULL;
397
18.5k
}
398
399
/* Compute the constraint of a facet of "set".
400
 *
401
 * We first compute the intersection with a bounding constraint
402
 * that is orthogonal to one of the coordinate axes.
403
 * If the affine hull of this intersection has only one equality,
404
 * we have found a facet.
405
 * Otherwise, we wrap the current bounding constraint around
406
 * one of the equalities of the face (one that is not equal to
407
 * the current bounding constraint).
408
 * This process continues until we have found a facet.
409
 * The dimension of the intersection increases by at least
410
 * one on each iteration, so termination is guaranteed.
411
 */
412
static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
413
3.88k
{
414
3.88k
  struct isl_set *slice = NULL;
415
3.88k
  struct isl_basic_set *face = NULL;
416
3.88k
  int i;
417
3.88k
  unsigned dim = isl_set_n_dim(set);
418
3.88k
  int is_bound;
419
3.88k
  isl_mat *bounds = NULL;
420
3.88k
421
3.88k
  isl_assert(set->ctx, set->n > 0, goto error);
422
3.88k
  bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
423
3.88k
  if (!bounds)
424
0
    return NULL;
425
3.88k
426
3.88k
  isl_seq_clr(bounds->row[0], dim);
427
3.88k
  isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
428
3.88k
  is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
429
3.88k
  if (is_bound < 0)
430
0
    goto error;
431
3.88k
  
isl_assert3.88k
(set->ctx, is_bound, goto error);
432
3.88k
  isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
433
3.88k
  bounds->n_row = 1;
434
3.88k
435
7.59k
  for (;;) {
436
7.59k
    slice = isl_set_copy(set);
437
7.59k
    slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
438
7.59k
    face = isl_set_affine_hull(slice);
439
7.59k
    if (!face)
440
0
      goto error;
441
7.59k
    
if (7.59k
face->n_eq == 17.59k
) {
442
3.88k
      isl_basic_set_free(face);
443
3.88k
      break;
444
3.88k
    }
445
6.98k
    
for (i = 0; 3.71k
i < face->n_eq6.98k
;
++i3.27k
)
446
6.98k
      
if (6.98k
!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
447
3.71k
          !isl_seq_is_neg(bounds->row[0],
448
3.71k
            face->eq[i], 1 + dim))
449
3.71k
        break;
450
3.71k
    isl_assert(set->ctx, i < face->n_eq, goto error);
451
3.71k
    
if (3.71k
!isl_set_wrap_facet(set, bounds->row[0], face->eq[i])3.71k
)
452
0
      goto error;
453
3.71k
    isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
454
3.71k
    isl_basic_set_free(face);
455
3.71k
  }
456
3.88k
457
3.88k
  return bounds;
458
0
error:
459
0
  isl_basic_set_free(face);
460
0
  isl_mat_free(bounds);
461
0
  return NULL;
462
3.88k
}
463
464
/* Given the bounding constraint "c" of a facet of the convex hull of "set",
465
 * compute a hyperplane description of the facet, i.e., compute the facets
466
 * of the facet.
467
 *
468
 * We compute an affine transformation that transforms the constraint
469
 *
470
 *        [ 1 ]
471
 *      c [ x ] = 0
472
 *
473
 * to the constraint
474
 *
475
 *         z_1  = 0
476
 *
477
 * by computing the right inverse U of a matrix that starts with the rows
478
 *
479
 *      [ 1 0 ]
480
 *      [  c  ]
481
 *
482
 * Then
483
 *      [ 1 ]     [ 1 ]
484
 *      [ x ] = U [ z ]
485
 * and
486
 *      [ 1 ]     [ 1 ]
487
 *      [ z ] = Q [ x ]
488
 *
489
 * with Q = U^{-1}
490
 * Since z_1 is zero, we can drop this variable as well as the corresponding
491
 * column of U to obtain
492
 *
493
 *      [ 1 ]      [ 1  ]
494
 *      [ x ] = U' [ z' ]
495
 * and
496
 *      [ 1  ]      [ 1 ]
497
 *      [ z' ] = Q' [ x ]
498
 *
499
 * with Q' equal to Q, but without the corresponding row.
500
 * After computing the facets of the facet in the z' space,
501
 * we convert them back to the x space through Q.
502
 */
503
static __isl_give isl_basic_set *compute_facet(__isl_keep isl_set *set,
504
  isl_int *c)
505
13.5k
{
506
13.5k
  struct isl_mat *m, *U, *Q;
507
13.5k
  struct isl_basic_set *facet = NULL;
508
13.5k
  struct isl_ctx *ctx;
509
13.5k
  unsigned dim;
510
13.5k
511
13.5k
  ctx = set->ctx;
512
13.5k
  set = isl_set_copy(set);
513
13.5k
  dim = isl_set_n_dim(set);
514
13.5k
  m = isl_mat_alloc(set->ctx, 2, 1 + dim);
515
13.5k
  if (!m)
516
0
    goto error;
517
13.5k
  
isl_int_set_si13.5k
(m->row[0][0], 1);
518
13.5k
  isl_seq_clr(m->row[0]+1, dim);
519
13.5k
  isl_seq_cpy(m->row[1], c, 1+dim);
520
13.5k
  U = isl_mat_right_inverse(m);
521
13.5k
  Q = isl_mat_right_inverse(isl_mat_copy(U));
522
13.5k
  U = isl_mat_drop_cols(U, 1, 1);
523
13.5k
  Q = isl_mat_drop_rows(Q, 1, 1);
524
13.5k
  set = isl_set_preimage(set, U);
525
13.5k
  facet = uset_convex_hull_wrap_bounded(set);
526
13.5k
  facet = isl_basic_set_preimage(facet, Q);
527
13.5k
  if (
facet && 13.5k
facet->n_eq != 013.5k
)
528
0
    isl_die(ctx, isl_error_internal, "unexpected equality",
529
13.5k
      return isl_basic_set_free(facet));
530
13.5k
  return facet;
531
0
error:
532
0
  isl_basic_set_free(facet);
533
0
  isl_set_free(set);
534
0
  return NULL;
535
13.5k
}
536
537
/* Given an initial facet constraint, compute the remaining facets.
538
 * We do this by running through all facets found so far and computing
539
 * the adjacent facets through wrapping, adding those facets that we
540
 * hadn't already found before.
541
 *
542
 * For each facet we have found so far, we first compute its facets
543
 * in the resulting convex hull.  That is, we compute the ridges
544
 * of the resulting convex hull contained in the facet.
545
 * We also compute the corresponding facet in the current approximation
546
 * of the convex hull.  There is no need to wrap around the ridges
547
 * in this facet since that would result in a facet that is already
548
 * present in the current approximation.
549
 *
550
 * This function can still be significantly optimized by checking which of
551
 * the facets of the basic sets are also facets of the convex hull and
552
 * using all the facets so far to help in constructing the facets of the
553
 * facets
554
 * and/or
555
 * using the technique in section "3.1 Ridge Generation" of
556
 * "Extended Convex Hull" by Fukuda et al.
557
 */
558
static __isl_give isl_basic_set *extend(__isl_take isl_basic_set *hull,
559
  __isl_keep isl_set *set)
560
3.90k
{
561
3.90k
  int i, j, f;
562
3.90k
  int k;
563
3.90k
  struct isl_basic_set *facet = NULL;
564
3.90k
  struct isl_basic_set *hull_facet = NULL;
565
3.90k
  unsigned dim;
566
3.90k
567
3.90k
  if (!hull)
568
0
    return NULL;
569
3.90k
570
3.90k
  
isl_assert3.90k
(set->ctx, set->n > 0, goto error);
571
3.90k
572
3.90k
  dim = isl_set_n_dim(set);
573
3.90k
574
17.4k
  for (i = 0; 
i < hull->n_ineq17.4k
;
++i13.5k
) {
575
13.5k
    facet = compute_facet(set, hull->ineq[i]);
576
13.5k
    facet = isl_basic_set_add_equality(facet, hull->ineq[i]);
577
13.5k
    facet = isl_basic_set_gauss(facet, NULL);
578
13.5k
    facet = isl_basic_set_normalize_constraints(facet);
579
13.5k
    hull_facet = isl_basic_set_copy(hull);
580
13.5k
    hull_facet = isl_basic_set_add_equality(hull_facet, hull->ineq[i]);
581
13.5k
    hull_facet = isl_basic_set_gauss(hull_facet, NULL);
582
13.5k
    hull_facet = isl_basic_set_normalize_constraints(hull_facet);
583
13.5k
    if (
!facet || 13.5k
!hull_facet13.5k
)
584
0
      goto error;
585
13.5k
    hull = isl_basic_set_cow(hull);
586
13.5k
    hull = isl_basic_set_extend_space(hull,
587
13.5k
      isl_space_copy(hull->dim), 0, 0, facet->n_ineq);
588
13.5k
    if (!hull)
589
0
      goto error;
590
48.2k
    
for (j = 0; 13.5k
j < facet->n_ineq48.2k
;
++j34.7k
) {
591
67.1k
      for (f = 0; 
f < hull_facet->n_ineq67.1k
;
++f32.4k
)
592
57.4k
        
if (57.4k
isl_seq_eq(facet->ineq[j],
593
57.4k
            hull_facet->ineq[f], 1 + dim))
594
25.0k
          break;
595
34.7k
      if (f < hull_facet->n_ineq)
596
25.0k
        continue;
597
9.65k
      k = isl_basic_set_alloc_inequality(hull);
598
9.65k
      if (k < 0)
599
0
        goto error;
600
9.65k
      isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
601
9.65k
      if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
602
0
        goto error;
603
34.7k
    }
604
13.5k
    isl_basic_set_free(hull_facet);
605
13.5k
    isl_basic_set_free(facet);
606
13.5k
  }
607
3.90k
  hull = isl_basic_set_simplify(hull);
608
3.90k
  hull = isl_basic_set_finalize(hull);
609
3.90k
  return hull;
610
0
error:
611
0
  isl_basic_set_free(hull_facet);
612
0
  isl_basic_set_free(facet);
613
0
  isl_basic_set_free(hull);
614
0
  return NULL;
615
3.90k
}
616
617
/* Special case for computing the convex hull of a one dimensional set.
618
 * We simply collect the lower and upper bounds of each basic set
619
 * and the biggest of those.
620
 */
621
static __isl_give isl_basic_set *convex_hull_1d(__isl_take isl_set *set)
622
4.06k
{
623
4.06k
  struct isl_mat *c = NULL;
624
4.06k
  isl_int *lower = NULL;
625
4.06k
  isl_int *upper = NULL;
626
4.06k
  int i, j, k;
627
4.06k
  isl_int a, b;
628
4.06k
  struct isl_basic_set *hull;
629
4.06k
630
12.2k
  for (i = 0; 
i < set->n12.2k
;
++i8.13k
) {
631
8.13k
    set->p[i] = isl_basic_set_simplify(set->p[i]);
632
8.13k
    if (!set->p[i])
633
0
      goto error;
634
8.13k
  }
635
4.06k
  set = isl_set_remove_empty_parts(set);
636
4.06k
  if (!set)
637
0
    goto error;
638
4.06k
  
isl_assert4.06k
(set->ctx, set->n > 0, goto error);
639
4.06k
  c = isl_mat_alloc(set->ctx, 2, 2);
640
4.06k
  if (!c)
641
0
    goto error;
642
4.06k
643
4.06k
  
if (4.06k
set->p[0]->n_eq > 04.06k
) {
644
4.06k
    isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error);
645
4.06k
    lower = c->row[0];
646
4.06k
    upper = c->row[1];
647
4.06k
    if (
isl_int_is_pos4.06k
(set->p[0]->eq[0][1])) {
648
4.06k
      isl_seq_cpy(lower, set->p[0]->eq[0], 2);
649
4.06k
      isl_seq_neg(upper, set->p[0]->eq[0], 2);
650
4.06k
    } else {
651
0
      isl_seq_neg(lower, set->p[0]->eq[0], 2);
652
0
      isl_seq_cpy(upper, set->p[0]->eq[0], 2);
653
0
    }
654
4.06k
  } else {
655
16
    for (j = 0; 
j < set->p[0]->n_ineq16
;
++j10
) {
656
10
      if (
isl_int_is_pos10
(set->p[0]->ineq[j][1])) {
657
6
        lower = c->row[0];
658
6
        isl_seq_cpy(lower, set->p[0]->ineq[j], 2);
659
10
      } else {
660
4
        upper = c->row[1];
661
4
        isl_seq_cpy(upper, set->p[0]->ineq[j], 2);
662
4
      }
663
10
    }
664
6
  }
665
4.06k
666
4.06k
  
isl_int_init4.06k
(a);
667
4.06k
  isl_int_init(b);
668
12.2k
  for (i = 0; 
i < set->n12.2k
;
++i8.13k
) {
669
8.13k
    struct isl_basic_set *bset = set->p[i];
670
8.13k
    int has_lower = 0;
671
8.13k
    int has_upper = 0;
672
8.13k
673
15.7k
    for (j = 0; 
j < bset->n_eq15.7k
;
++j7.59k
) {
674
7.59k
      has_lower = 1;
675
7.59k
      has_upper = 1;
676
7.59k
      if (
lower7.59k
) {
677
7.59k
        isl_int_mul(a, lower[0], bset->eq[j][1]);
678
7.59k
        isl_int_mul(b, lower[1], bset->eq[j][0]);
679
7.59k
        if (
isl_int_lt7.59k
(a, b) && 7.59k
isl_int_is_pos1.64k
(bset->eq[j][1]))
680
1.64k
          isl_seq_cpy(lower, bset->eq[j], 2);
681
7.59k
        if (
isl_int_gt7.59k
(a, b) && 7.59k
isl_int_is_neg1.89k
(bset->eq[j][1]))
682
0
          isl_seq_neg(lower, bset->eq[j], 2);
683
7.59k
      }
684
7.59k
      if (
upper7.59k
) {
685
7.59k
        isl_int_mul(a, upper[0], bset->eq[j][1]);
686
7.59k
        isl_int_mul(b, upper[1], bset->eq[j][0]);
687
7.59k
        if (
isl_int_lt7.59k
(a, b) && 7.59k
isl_int_is_pos1.89k
(bset->eq[j][1]))
688
1.89k
          isl_seq_neg(upper, bset->eq[j], 2);
689
7.59k
        if (
isl_int_gt7.59k
(a, b) && 7.59k
isl_int_is_neg1.64k
(bset->eq[j][1]))
690
0
          isl_seq_cpy(upper, bset->eq[j], 2);
691
7.59k
      }
692
7.59k
    }
693
9.20k
    for (j = 0; 
j < bset->n_ineq9.20k
;
++j1.07k
) {
694
1.07k
      if (isl_int_is_pos(bset->ineq[j][1]))
695
537
        has_lower = 1;
696
1.07k
      if (isl_int_is_neg(bset->ineq[j][1]))
697
535
        has_upper = 1;
698
1.07k
      if (
lower && 1.07k
isl_int_is_pos1.07k
(bset->ineq[j][1])) {
699
537
        isl_int_mul(a, lower[0], bset->ineq[j][1]);
700
537
        isl_int_mul(b, lower[1], bset->ineq[j][0]);
701
537
        if (isl_int_lt(a, b))
702
286
          isl_seq_cpy(lower, bset->ineq[j], 2);
703
537
      }
704
1.07k
      if (
upper && 1.07k
isl_int_is_neg1.07k
(bset->ineq[j][1])) {
705
535
        isl_int_mul(a, upper[0], bset->ineq[j][1]);
706
535
        isl_int_mul(b, upper[1], bset->ineq[j][0]);
707
535
        if (isl_int_gt(a, b))
708
245
          isl_seq_cpy(upper, bset->ineq[j], 2);
709
535
      }
710
1.07k
    }
711
8.13k
    if (!has_lower)
712
0
      lower = NULL;
713
8.13k
    if (!has_upper)
714
2
      upper = NULL;
715
8.13k
  }
716
4.06k
  isl_int_clear(a);
717
4.06k
  isl_int_clear(b);
718
4.06k
719
4.06k
  hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2);
720
4.06k
  hull = isl_basic_set_set_rational(hull);
721
4.06k
  if (!hull)
722
0
    goto error;
723
4.06k
  
if (4.06k
lower4.06k
) {
724
4.06k
    k = isl_basic_set_alloc_inequality(hull);
725
4.06k
    isl_seq_cpy(hull->ineq[k], lower, 2);
726
4.06k
  }
727
4.06k
  if (
upper4.06k
) {
728
4.06k
    k = isl_basic_set_alloc_inequality(hull);
729
4.06k
    isl_seq_cpy(hull->ineq[k], upper, 2);
730
4.06k
  }
731
4.06k
  hull = isl_basic_set_finalize(hull);
732
4.06k
  isl_set_free(set);
733
4.06k
  isl_mat_free(c);
734
4.06k
  return hull;
735
0
error:
736
0
  isl_set_free(set);
737
0
  isl_mat_free(c);
738
0
  return NULL;
739
4.06k
}
740
741
static __isl_give isl_basic_set *convex_hull_0d(__isl_take isl_set *set)
742
0
{
743
0
  struct isl_basic_set *convex_hull;
744
0
745
0
  if (!set)
746
0
    return NULL;
747
0
748
0
  
if (0
isl_set_is_empty(set)0
)
749
0
    convex_hull = isl_basic_set_empty(isl_space_copy(set->dim));
750
0
  else
751
0
    convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
752
0
  isl_set_free(set);
753
0
  return convex_hull;
754
0
}
755
756
/* Compute the convex hull of a pair of basic sets without any parameters or
757
 * integer divisions using Fourier-Motzkin elimination.
758
 * The convex hull is the set of all points that can be written as
759
 * the sum of points from both basic sets (in homogeneous coordinates).
760
 * We set up the constraints in a space with dimensions for each of
761
 * the three sets and then project out the dimensions corresponding
762
 * to the two original basic sets, retaining only those corresponding
763
 * to the convex hull.
764
 */
765
static __isl_give isl_basic_set *convex_hull_pair_elim(
766
  __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
767
14
{
768
14
  int i, j, k;
769
14
  struct isl_basic_set *bset[2];
770
14
  struct isl_basic_set *hull = NULL;
771
14
  unsigned dim;
772
14
773
14
  if (
!bset1 || 14
!bset214
)
774
0
    goto error;
775
14
776
14
  dim = isl_basic_set_n_dim(bset1);
777
14
  hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0,
778
14
        1 + dim + bset1->n_eq + bset2->n_eq,
779
14
        2 + bset1->n_ineq + bset2->n_ineq);
780
14
  bset[0] = bset1;
781
14
  bset[1] = bset2;
782
42
  for (i = 0; 
i < 242
;
++i28
) {
783
41
    for (j = 0; 
j < bset[i]->n_eq41
;
++j13
) {
784
13
      k = isl_basic_set_alloc_equality(hull);
785
13
      if (k < 0)
786
0
        goto error;
787
13
      isl_seq_clr(hull->eq[k], (i+1) * (1+dim));
788
13
      isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
789
13
      isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j],
790
13
          1+dim);
791
13
    }
792
126
    
for (j = 0; 28
j < bset[i]->n_ineq126
;
++j98
) {
793
98
      k = isl_basic_set_alloc_inequality(hull);
794
98
      if (k < 0)
795
0
        goto error;
796
98
      isl_seq_clr(hull->ineq[k], (i+1) * (1+dim));
797
98
      isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim));
798
98
      isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim),
799
98
          bset[i]->ineq[j], 1+dim);
800
98
    }
801
28
    k = isl_basic_set_alloc_inequality(hull);
802
28
    if (k < 0)
803
0
      goto error;
804
28
    isl_seq_clr(hull->ineq[k], 1+2+3*dim);
805
28
    isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1);
806
28
  }
807
65
  
for (j = 0; 14
j < 1+dim65
;
++j51
) {
808
51
    k = isl_basic_set_alloc_equality(hull);
809
51
    if (k < 0)
810
0
      goto error;
811
51
    isl_seq_clr(hull->eq[k], 1+2+3*dim);
812
51
    isl_int_set_si(hull->eq[k][j], -1);
813
51
    isl_int_set_si(hull->eq[k][1+dim+j], 1);
814
51
    isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1);
815
51
  }
816
14
  hull = isl_basic_set_set_rational(hull);
817
14
  hull = isl_basic_set_remove_dims(hull, isl_dim_set, dim, 2*(1+dim));
818
14
  hull = isl_basic_set_remove_redundancies(hull);
819
14
  isl_basic_set_free(bset1);
820
14
  isl_basic_set_free(bset2);
821
14
  return hull;
822
0
error:
823
0
  isl_basic_set_free(bset1);
824
0
  isl_basic_set_free(bset2);
825
0
  isl_basic_set_free(hull);
826
0
  return NULL;
827
14
}
828
829
/* Is the set bounded for each value of the parameters?
830
 */
831
isl_bool isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset)
832
2.83k
{
833
2.83k
  struct isl_tab *tab;
834
2.83k
  isl_bool bounded;
835
2.83k
836
2.83k
  if (!bset)
837
0
    return isl_bool_error;
838
2.83k
  
if (2.83k
isl_basic_set_plain_is_empty(bset)2.83k
)
839
0
    return isl_bool_true;
840
2.83k
841
2.83k
  tab = isl_tab_from_recession_cone(bset, 1);
842
2.83k
  bounded = isl_tab_cone_is_bounded(tab);
843
2.83k
  isl_tab_free(tab);
844
2.83k
  return bounded;
845
2.83k
}
846
847
/* Is the image bounded for each value of the parameters and
848
 * the domain variables?
849
 */
850
isl_bool isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap)
851
0
{
852
0
  unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param);
853
0
  unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in);
854
0
  isl_bool bounded;
855
0
856
0
  bmap = isl_basic_map_copy(bmap);
857
0
  bmap = isl_basic_map_cow(bmap);
858
0
  bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam,
859
0
          isl_dim_in, 0, n_in);
860
0
  bounded = isl_basic_set_is_bounded(bset_from_bmap(bmap));
861
0
  isl_basic_map_free(bmap);
862
0
863
0
  return bounded;
864
0
}
865
866
/* Is the set bounded for each value of the parameters?
867
 */
868
isl_bool isl_set_is_bounded(__isl_keep isl_set *set)
869
234
{
870
234
  int i;
871
234
872
234
  if (!set)
873
0
    return isl_bool_error;
874
234
875
712
  
for (i = 0; 234
i < set->n712
;
++i478
) {
876
510
    isl_bool bounded = isl_basic_set_is_bounded(set->p[i]);
877
510
    if (
!bounded || 510
bounded < 0478
)
878
32
      return bounded;
879
510
  }
880
202
  return isl_bool_true;
881
234
}
882
883
/* Compute the lineality space of the convex hull of bset1 and bset2.
884
 *
885
 * We first compute the intersection of the recession cone of bset1
886
 * with the negative of the recession cone of bset2 and then compute
887
 * the linear hull of the resulting cone.
888
 */
889
static __isl_give isl_basic_set *induced_lineality_space(
890
  __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
891
4
{
892
4
  int i, k;
893
4
  struct isl_basic_set *lin = NULL;
894
4
  unsigned dim;
895
4
896
4
  if (
!bset1 || 4
!bset24
)
897
0
    goto error;
898
4
899
4
  dim = isl_basic_set_total_dim(bset1);
900
4
  lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset1), 0,
901
4
          bset1->n_eq + bset2->n_eq,
902
4
          bset1->n_ineq + bset2->n_ineq);
903
4
  lin = isl_basic_set_set_rational(lin);
904
4
  if (!lin)
905
0
    goto error;
906
4
  
for (i = 0; 4
i < bset1->n_eq4
;
++i0
) {
907
0
    k = isl_basic_set_alloc_equality(lin);
908
0
    if (k < 0)
909
0
      goto error;
910
0
    
isl_int_set_si0
(lin->eq[k][0], 0);
911
0
    isl_seq_cpy(lin->eq[k] + 1, bset1->eq[i] + 1, dim);
912
0
  }
913
25
  
for (i = 0; 4
i < bset1->n_ineq25
;
++i21
) {
914
21
    k = isl_basic_set_alloc_inequality(lin);
915
21
    if (k < 0)
916
0
      goto error;
917
21
    
isl_int_set_si21
(lin->ineq[k][0], 0);
918
21
    isl_seq_cpy(lin->ineq[k] + 1, bset1->ineq[i] + 1, dim);
919
21
  }
920
4
  
for (i = 0; 4
i < bset2->n_eq4
;
++i0
) {
921
0
    k = isl_basic_set_alloc_equality(lin);
922
0
    if (k < 0)
923
0
      goto error;
924
0
    
isl_int_set_si0
(lin->eq[k][0], 0);
925
0
    isl_seq_neg(lin->eq[k] + 1, bset2->eq[i] + 1, dim);
926
0
  }
927
29
  
for (i = 0; 4
i < bset2->n_ineq29
;
++i25
) {
928
25
    k = isl_basic_set_alloc_inequality(lin);
929
25
    if (k < 0)
930
0
      goto error;
931
25
    
isl_int_set_si25
(lin->ineq[k][0], 0);
932
25
    isl_seq_neg(lin->ineq[k] + 1, bset2->ineq[i] + 1, dim);
933
25
  }
934
4
935
4
  isl_basic_set_free(bset1);
936
4
  isl_basic_set_free(bset2);
937
4
  return isl_basic_set_affine_hull(lin);
938
0
error:
939
0
  isl_basic_set_free(lin);
940
0
  isl_basic_set_free(bset1);
941
0
  isl_basic_set_free(bset2);
942
0
  return NULL;
943
4
}
944
945
static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set);
946
947
/* Given a set and a linear space "lin" of dimension n > 0,
948
 * project the linear space from the set, compute the convex hull
949
 * and then map the set back to the original space.
950
 *
951
 * Let
952
 *
953
 *  M x = 0
954
 *
955
 * describe the linear space.  We first compute the Hermite normal
956
 * form H = M U of M = H Q, to obtain
957
 *
958
 *  H Q x = 0
959
 *
960
 * The last n rows of H will be zero, so the last n variables of x' = Q x
961
 * are the one we want to project out.  We do this by transforming each
962
 * basic set A x >= b to A U x' >= b and then removing the last n dimensions.
963
 * After computing the convex hull in x'_1, i.e., A' x'_1 >= b',
964
 * we transform the hull back to the original space as A' Q_1 x >= b',
965
 * with Q_1 all but the last n rows of Q.
966
 */
967
static __isl_give isl_basic_set *modulo_lineality(__isl_take isl_set *set,
968
  __isl_take isl_basic_set *lin)
969
12
{
970
12
  unsigned total = isl_basic_set_total_dim(lin);
971
12
  unsigned lin_dim;
972
12
  struct isl_basic_set *hull;
973
12
  struct isl_mat *M, *U, *Q;
974
12
975
12
  if (
!set || 12
!lin12
)
976
0
    goto error;
977
12
  lin_dim = total - lin->n_eq;
978
12
  M = isl_mat_sub_alloc6(set->ctx, lin->eq, 0, lin->n_eq, 1, total);
979
12
  M = isl_mat_left_hermite(M, 0, &U, &Q);
980
12
  if (!M)
981
0
    goto error;
982
12
  isl_mat_free(M);
983
12
  isl_basic_set_free(lin);
984
12
985
12
  Q = isl_mat_drop_rows(Q, Q->n_row - lin_dim, lin_dim);
986
12
987
12
  U = isl_mat_lin_to_aff(U);
988
12
  Q = isl_mat_lin_to_aff(Q);
989
12
990
12
  set = isl_set_preimage(set, U);
991
12
  set = isl_set_remove_dims(set, isl_dim_set, total - lin_dim, lin_dim);
992
12
  hull = uset_convex_hull(set);
993
12
  hull = isl_basic_set_preimage(hull, Q);
994
12
995
12
  return hull;
996
0
error:
997
0
  isl_basic_set_free(lin);
998
0
  isl_set_free(set);
999
0
  return NULL;
1000
12
}
1001
1002
/* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space,
1003
 * set up an LP for solving
1004
 *
1005
 *  \sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j}
1006
 *
1007
 * \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0
1008
 * The next \alpha{ij} correspond to the equalities and come in pairs.
1009
 * The final \alpha{ij} correspond to the inequalities.
1010
 */
1011
static __isl_give isl_basic_set *valid_direction_lp(
1012
  __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1013
8
{
1014
8
  isl_space *dim;
1015
8
  struct isl_basic_set *lp;
1016
8
  unsigned d;
1017
8
  int n;
1018
8
  int i, j, k;
1019
8
1020
8
  if (
!bset1 || 8
!bset28
)
1021
0
    goto error;
1022
8
  d = 1 + isl_basic_set_total_dim(bset1);
1023
8
  n = 2 +
1024
8
      2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq;
1025
8
  dim = isl_space_set_alloc(bset1->ctx, 0, n);
1026
8
  lp = isl_basic_set_alloc_space(dim, 0, d, n);
1027
8
  if (!lp)
1028
0
    goto error;
1029
105
  
for (i = 0; 8
i < n105
;
++i97
) {
1030
97
    k = isl_basic_set_alloc_inequality(lp);
1031
97
    if (k < 0)
1032
0
      goto error;
1033
97
    isl_seq_clr(lp->ineq[k] + 1, n);
1034
97
    isl_int_set_si(lp->ineq[k][0], -1);
1035
97
    isl_int_set_si(lp->ineq[k][1 + i], 1);
1036
97
  }
1037
38
  
for (i = 0; 8
i < d38
;
++i30
) {
1038
30
    k = isl_basic_set_alloc_equality(lp);
1039
30
    if (k < 0)
1040
0
      goto error;
1041
30
    n = 0;
1042
30
    isl_int_set_si(lp->eq[k][n], 0); n++;
1043
30
    /* positivity constraint 1 >= 0 */
1044
30
    isl_int_set_si(lp->eq[k][n], i == 0); n++;
1045
44
    for (j = 0; 
j < bset1->n_eq44
;
++j14
) {
1046
14
      isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++;
1047
14
      isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++;
1048
14
    }
1049
180
    for (j = 0; 
j < bset1->n_ineq180
;
++j150
) {
1050
150
      isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++;
1051
150
    }
1052
30
    /* positivity constraint 1 >= 0 */
1053
30
    isl_int_set_si(lp->eq[k][n], -(i == 0)); n++;
1054
44
    for (j = 0; 
j < bset2->n_eq44
;
++j14
) {
1055
14
      isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++;
1056
14
      isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++;
1057
14
    }
1058
204
    for (j = 0; 
j < bset2->n_ineq204
;
++j174
) {
1059
174
      isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++;
1060
174
    }
1061
30
  }
1062
8
  lp = isl_basic_set_gauss(lp, NULL);
1063
8
  isl_basic_set_free(bset1);
1064
8
  isl_basic_set_free(bset2);
1065
8
  return lp;
1066
0
error:
1067
0
  isl_basic_set_free(bset1);
1068
0
  isl_basic_set_free(bset2);
1069
0
  return NULL;
1070
8
}
1071
1072
/* Compute a vector s in the homogeneous space such that <s, r> > 0
1073
 * for all rays in the homogeneous space of the two cones that correspond
1074
 * to the input polyhedra bset1 and bset2.
1075
 *
1076
 * We compute s as a vector that satisfies
1077
 *
1078
 *  s = \sum_j \alpha_{ij} h_{ij} for i = 1,2     (*)
1079
 *
1080
 * with h_{ij} the normals of the facets of polyhedron i
1081
 * (including the "positivity constraint" 1 >= 0) and \alpha_{ij}
1082
 * strictly positive numbers.  For simplicity we impose \alpha_{ij} >= 1.
1083
 * We first set up an LP with as variables the \alpha{ij}.
1084
 * In this formulation, for each polyhedron i,
1085
 * the first constraint is the positivity constraint, followed by pairs
1086
 * of variables for the equalities, followed by variables for the inequalities.
1087
 * We then simply pick a feasible solution and compute s using (*).
1088
 *
1089
 * Note that we simply pick any valid direction and make no attempt
1090
 * to pick a "good" or even the "best" valid direction.
1091
 */
1092
static __isl_give isl_vec *valid_direction(
1093
  __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1094
8
{
1095
8
  struct isl_basic_set *lp;
1096
8
  struct isl_tab *tab;
1097
8
  struct isl_vec *sample = NULL;
1098
8
  struct isl_vec *dir;
1099
8
  unsigned d;
1100
8
  int i;
1101
8
  int n;
1102
8
1103
8
  if (
!bset1 || 8
!bset28
)
1104
0
    goto error;
1105
8
  lp = valid_direction_lp(isl_basic_set_copy(bset1),
1106
8
        isl_basic_set_copy(bset2));
1107
8
  tab = isl_tab_from_basic_set(lp, 0);
1108
8
  sample = isl_tab_get_sample_value(tab);
1109
8
  isl_tab_free(tab);
1110
8
  isl_basic_set_free(lp);
1111
8
  if (!sample)
1112
0
    goto error;
1113
8
  d = isl_basic_set_total_dim(bset1);
1114
8
  dir = isl_vec_alloc(bset1->ctx, 1 + d);
1115
8
  if (!dir)
1116
0
    goto error;
1117
8
  isl_seq_clr(dir->block.data + 1, dir->size - 1);
1118
8
  n = 1;
1119
8
  /* positivity constraint 1 >= 0 */
1120
8
  isl_int_set(dir->block.data[0], sample->block.data[n]); n++;
1121
12
  for (i = 0; 
i < bset1->n_eq12
;
++i4
) {
1122
4
    isl_int_sub(sample->block.data[n],
1123
4
          sample->block.data[n], sample->block.data[n+1]);
1124
4
    isl_seq_combine(dir->block.data,
1125
4
        bset1->ctx->one, dir->block.data,
1126
4
        sample->block.data[n], bset1->eq[i], 1 + d);
1127
4
1128
4
    n += 2;
1129
4
  }
1130
39
  for (i = 0; 
i < bset1->n_ineq39
;
++i31
)
1131
31
    isl_seq_combine(dir->block.data,
1132
31
        bset1->ctx->one, dir->block.data,
1133
31
        sample->block.data[n++], bset1->ineq[i], 1 + d);
1134
8
  isl_vec_free(sample);
1135
8
  isl_seq_normalize(bset1->ctx, dir->el, dir->size);
1136
8
  isl_basic_set_free(bset1);
1137
8
  isl_basic_set_free(bset2);
1138
8
  return dir;
1139
0
error:
1140
0
  isl_vec_free(sample);
1141
0
  isl_basic_set_free(bset1);
1142
0
  isl_basic_set_free(bset2);
1143
0
  return NULL;
1144
8
}
1145
1146
/* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1},
1147
 * compute b_i' + A_i' x' >= 0, with
1148
 *
1149
 *  [ b_i A_i ]        [ y' ]                 [ y' ]
1150
 *  [  1   0  ] S^{-1} [ x' ] >= 0  or  [ b_i' A_i' ] [ x' ] >= 0
1151
 *
1152
 * In particular, add the "positivity constraint" and then perform
1153
 * the mapping.
1154
 */
1155
static __isl_give isl_basic_set *homogeneous_map(__isl_take isl_basic_set *bset,
1156
  __isl_take isl_mat *T)
1157
16
{
1158
16
  int k;
1159
16
1160
16
  if (!bset)
1161
0
    goto error;
1162
16
  bset = isl_basic_set_extend_constraints(bset, 0, 1);
1163
16
  k = isl_basic_set_alloc_inequality(bset);
1164
16
  if (k < 0)
1165
0
    goto error;
1166
16
  isl_seq_clr(bset->ineq[k] + 1, isl_basic_set_total_dim(bset));
1167
16
  isl_int_set_si(bset->ineq[k][0], 1);
1168
16
  bset = isl_basic_set_preimage(bset, T);
1169
16
  return bset;
1170
0
error:
1171
0
  isl_mat_free(T);
1172
0
  isl_basic_set_free(bset);
1173
0
  return NULL;
1174
16
}
1175
1176
/* Compute the convex hull of a pair of basic sets without any parameters or
1177
 * integer divisions, where the convex hull is known to be pointed,
1178
 * but the basic sets may be unbounded.
1179
 *
1180
 * We turn this problem into the computation of a convex hull of a pair
1181
 * _bounded_ polyhedra by "changing the direction of the homogeneous
1182
 * dimension".  This idea is due to Matthias Koeppe.
1183
 *
1184
 * Consider the cones in homogeneous space that correspond to the
1185
 * input polyhedra.  The rays of these cones are also rays of the
1186
 * polyhedra if the coordinate that corresponds to the homogeneous
1187
 * dimension is zero.  That is, if the inner product of the rays
1188
 * with the homogeneous direction is zero.
1189
 * The cones in the homogeneous space can also be considered to
1190
 * correspond to other pairs of polyhedra by chosing a different
1191
 * homogeneous direction.  To ensure that both of these polyhedra
1192
 * are bounded, we need to make sure that all rays of the cones
1193
 * correspond to vertices and not to rays.
1194
 * Let s be a direction such that <s, r> > 0 for all rays r of both cones.
1195
 * Then using s as a homogeneous direction, we obtain a pair of polytopes.
1196
 * The vector s is computed in valid_direction.
1197
 *
1198
 * Note that we need to consider _all_ rays of the cones and not just
1199
 * the rays that correspond to rays in the polyhedra.  If we were to
1200
 * only consider those rays and turn them into vertices, then we
1201
 * may inadvertently turn some vertices into rays.
1202
 *
1203
 * The standard homogeneous direction is the unit vector in the 0th coordinate.
1204
 * We therefore transform the two polyhedra such that the selected
1205
 * direction is mapped onto this standard direction and then proceed
1206
 * with the normal computation.
1207
 * Let S be a non-singular square matrix with s as its first row,
1208
 * then we want to map the polyhedra to the space
1209
 *
1210
 *  [ y' ]     [ y ]    [ y ]          [ y' ]
1211
 *  [ x' ] = S [ x ]  i.e., [ x ] = S^{-1} [ x' ]
1212
 *
1213
 * We take S to be the unimodular completion of s to limit the growth
1214
 * of the coefficients in the following computations.
1215
 *
1216
 * Let b_i + A_i x >= 0 be the constraints of polyhedron i.
1217
 * We first move to the homogeneous dimension
1218
 *
1219
 *  b_i y + A_i x >= 0    [ b_i A_i ] [ y ]    [ 0 ]
1220
 *      y         >= 0  or  [  1   0  ] [ x ] >= [ 0 ]
1221
 *
1222
 * Then we change directoin
1223
 *
1224
 *  [ b_i A_i ]        [ y' ]                 [ y' ]
1225
 *  [  1   0  ] S^{-1} [ x' ] >= 0  or  [ b_i' A_i' ] [ x' ] >= 0
1226
 *
1227
 * Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 0
1228
 * resulting in b' + A' x' >= 0, which we then convert back
1229
 *
1230
 *              [ y ]           [ y ]
1231
 *  [ b' A' ] S [ x ] >= 0  or  [ b A ] [ x ] >= 0
1232
 *
1233
 * The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra.
1234
 */
1235
static __isl_give isl_basic_set *convex_hull_pair_pointed(
1236
  __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1237
8
{
1238
8
  struct isl_ctx *ctx = NULL;
1239
8
  struct isl_vec *dir = NULL;
1240
8
  struct isl_mat *T = NULL;
1241
8
  struct isl_mat *T2 = NULL;
1242
8
  struct isl_basic_set *hull;
1243
8
  struct isl_set *set;
1244
8
1245
8
  if (
!bset1 || 8
!bset28
)
1246
0
    goto error;
1247
8
  ctx = isl_basic_set_get_ctx(bset1);
1248
8
  dir = valid_direction(isl_basic_set_copy(bset1),
1249
8
        isl_basic_set_copy(bset2));
1250
8
  if (!dir)
1251
0
    goto error;
1252
8
  T = isl_mat_alloc(ctx, dir->size, dir->size);
1253
8
  if (!T)
1254
0
    goto error;
1255
8
  isl_seq_cpy(T->row[0], dir->block.data, dir->size);
1256
8
  T = isl_mat_unimodular_complete(T, 1);
1257
8
  T2 = isl_mat_right_inverse(isl_mat_copy(T));
1258
8
1259
8
  bset1 = homogeneous_map(bset1, isl_mat_copy(T2));
1260
8
  bset2 = homogeneous_map(bset2, T2);
1261
8
  set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1262
8
  set = isl_set_add_basic_set(set, bset1);
1263
8
  set = isl_set_add_basic_set(set, bset2);
1264
8
  hull = uset_convex_hull(set);
1265
8
  hull = isl_basic_set_preimage(hull, T);
1266
8
   
1267
8
  isl_vec_free(dir);
1268
8
1269
8
  return hull;
1270
0
error:
1271
0
  isl_vec_free(dir);
1272
0
  isl_basic_set_free(bset1);
1273
0
  isl_basic_set_free(bset2);
1274
0
  return NULL;
1275
8
}
1276
1277
static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set);
1278
static __isl_give isl_basic_set *modulo_affine_hull(
1279
  __isl_take isl_set *set, __isl_take isl_basic_set *affine_hull);
1280
1281
/* Compute the convex hull of a pair of basic sets without any parameters or
1282
 * integer divisions.
1283
 *
1284
 * This function is called from uset_convex_hull_unbounded, which
1285
 * means that the complete convex hull is unbounded.  Some pairs
1286
 * of basic sets may still be bounded, though.
1287
 * They may even lie inside a lower dimensional space, in which
1288
 * case they need to be handled inside their affine hull since
1289
 * the main algorithm assumes that the result is full-dimensional.
1290
 *
1291
 * If the convex hull of the two basic sets would have a non-trivial
1292
 * lineality space, we first project out this lineality space.
1293
 */
1294
static __isl_give isl_basic_set *convex_hull_pair(
1295
  __isl_take isl_basic_set *bset1, __isl_take isl_basic_set *bset2)
1296
26
{
1297
26
  isl_basic_set *lin, *aff;
1298
26
  int bounded1, bounded2;
1299
26
1300
26
  if (
bset1->ctx->opt->convex == 26
ISL_CONVEX_HULL_FM26
)
1301
14
    return convex_hull_pair_elim(bset1, bset2);
1302
12
1303
12
  aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1),
1304
12
                isl_basic_set_copy(bset2)));
1305
12
  if (!aff)
1306
0
    goto error;
1307
12
  
if (12
aff->n_eq != 012
)
1308
1
    return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff);
1309
11
  isl_basic_set_free(aff);
1310
11
1311
11
  bounded1 = isl_basic_set_is_bounded(bset1);
1312
11
  bounded2 = isl_basic_set_is_bounded(bset2);
1313
11
1314
11
  if (
bounded1 < 0 || 11
bounded2 < 011
)
1315
0
    goto error;
1316
11
1317
11
  
if (11
bounded1 && 11
bounded24
)
1318
1
    return uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2));
1319
10
1320
10
  
if (10
bounded1 || 10
bounded27
)
1321
6
    return convex_hull_pair_pointed(bset1, bset2);
1322
4
1323
4
  lin = induced_lineality_space(isl_basic_set_copy(bset1),
1324
4
              isl_basic_set_copy(bset2));
1325
4
  if (!lin)
1326
0
    goto error;
1327
4
  
if (4
isl_basic_set_plain_is_universe(lin)4
) {
1328
0
    isl_basic_set_free(bset1);
1329
0
    isl_basic_set_free(bset2);
1330
0
    return lin;
1331
0
  }
1332
4
  
if (4
lin->n_eq < isl_basic_set_total_dim(lin)4
) {
1333
2
    struct isl_set *set;
1334
2
    set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0);
1335
2
    set = isl_set_add_basic_set(set, bset1);
1336
2
    set = isl_set_add_basic_set(set, bset2);
1337
2
    return modulo_lineality(set, lin);
1338
2
  }
1339
2
  isl_basic_set_free(lin);
1340
2
1341
2
  return convex_hull_pair_pointed(bset1, bset2);
1342
0
error:
1343
0
  isl_basic_set_free(bset1);
1344
0
  isl_basic_set_free(bset2);
1345
0
  return NULL;
1346
26
}
1347
1348
/* Compute the lineality space of a basic set.
1349
 * We basically just drop the constants and turn every inequality
1350
 * into an equality.
1351
 * Any explicit representations of local variables are removed
1352
 * because they may no longer be valid representations
1353
 * in the lineality space.
1354
 */
1355
__isl_give isl_basic_set *isl_basic_set_lineality_space(
1356
  __isl_take isl_basic_set *bset)
1357
120
{
1358
120
  int i, k;
1359
120
  struct isl_basic_set *lin = NULL;
1360
120
  unsigned n_div, dim;
1361
120
1362
120
  if (!bset)
1363
0
    goto error;
1364
120
  n_div = isl_basic_set_dim(bset, isl_dim_div);
1365
120
  dim = isl_basic_set_total_dim(bset);
1366
120
1367
120
  lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset),
1368
120
          n_div, dim, 0);
1369
120
  for (i = 0; 
i < n_div120
;
++i0
)
1370
0
    
if (0
isl_basic_set_alloc_div(lin) < 00
)
1371
0
      goto error;
1372
120
  
if (120
!lin120
)
1373
0
    goto error;
1374
238
  
for (i = 0; 120
i < bset->n_eq238
;
++i118
) {
1375
118
    k = isl_basic_set_alloc_equality(lin);
1376
118
    if (k < 0)
1377
0
      goto error;
1378
118
    
isl_int_set_si118
(lin->eq[k][0], 0);
1379
118
    isl_seq_cpy(lin->eq[k] + 1, bset->eq[i] + 1, dim);
1380
118
  }
1381
120
  lin = isl_basic_set_gauss(lin, NULL);
1382
120
  if (!lin)
1383
0
    goto error;
1384
306
  
for (i = 0; 120
i < bset->n_ineq && 306
lin->n_eq < dim223
;
++i186
) {
1385
186
    k = isl_basic_set_alloc_equality(lin);
1386
186
    if (k < 0)
1387
0
      goto error;
1388
186
    
isl_int_set_si186
(lin->eq[k][0], 0);
1389
186
    isl_seq_cpy(lin->eq[k] + 1, bset->ineq[i] + 1, dim);
1390
186
    lin = isl_basic_set_gauss(lin, NULL);
1391
186
    if (!lin)
1392
0
      goto error;
1393
186
  }
1394
120
  isl_basic_set_free(bset);
1395
120
  return lin;
1396
0
error:
1397
0
  isl_basic_set_free(lin);
1398
0
  isl_basic_set_free(bset);
1399
0
  return NULL;
1400
120
}
1401
1402
/* Compute the (linear) hull of the lineality spaces of the basic sets in the
1403
 * set "set".
1404
 */
1405
__isl_give isl_basic_set *isl_set_combined_lineality_space(
1406
  __isl_take isl_set *set)
1407
56
{
1408
56
  int i;
1409
56
  struct isl_set *lin = NULL;
1410
56
1411
56
  if (!set)
1412
0
    return NULL;
1413
56
  
if (56
set->n == 056
) {
1414
0
    isl_space *space = isl_set_get_space(set);
1415
0
    isl_set_free(set);
1416
0
    return isl_basic_set_empty(space);
1417
0
  }
1418
56
1419
56
  lin = isl_set_alloc_space(isl_set_get_space(set), set->n, 0);
1420
172
  for (i = 0; 
i < set->n172
;
++i116
)
1421
116
    lin = isl_set_add_basic_set(lin,
1422
116
        isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i])));
1423
56
  isl_set_free(set);
1424
56
  return isl_set_affine_hull(lin);
1425
56
}
1426
1427
/* Compute the convex hull of a set without any parameters or
1428
 * integer divisions.
1429
 * In each step, we combined two basic sets until only one
1430
 * basic set is left.
1431
 * The input basic sets are assumed not to have a non-trivial
1432
 * lineality space.  If any of the intermediate results has
1433
 * a non-trivial lineality space, it is projected out.
1434
 */
1435
static __isl_give isl_basic_set *uset_convex_hull_unbounded(
1436
  __isl_take isl_set *set)
1437
24
{
1438
24
  isl_basic_set_list *list;
1439
24
1440
24
  list = isl_set_get_basic_set_list(set);
1441
24
  isl_set_free(set);
1442
24
1443
26
  while (
list26
) {
1444
26
    int n;
1445
26
    struct isl_basic_set *t;
1446
26
    isl_basic_set *bset1, *bset2;
1447
26
1448
26
    n = isl_basic_set_list_n_basic_set(list);
1449
26
    if (n < 2)
1450
0
      isl_die(isl_basic_set_list_get_ctx(list),
1451
26
        isl_error_internal,
1452
26
        "expecting at least two elements", goto error);
1453
26
    bset1 = isl_basic_set_list_get_basic_set(list, n - 1);
1454
26
    bset2 = isl_basic_set_list_get_basic_set(list, n - 2);
1455
26
    bset1 = convex_hull_pair(bset1, bset2);
1456
26
    if (
n == 226
) {
1457
22
      isl_basic_set_list_free(list);
1458
22
      return bset1;
1459
22
    }
1460
4
    bset1 = isl_basic_set_underlying_set(bset1);
1461
4
    list = isl_basic_set_list_drop(list, n - 2, 2);
1462
4
    list = isl_basic_set_list_add(list, bset1);
1463
4
1464
4
    t = isl_basic_set_list_get_basic_set(list, n - 2);
1465
4
    t = isl_basic_set_lineality_space(t);
1466
4
    if (!t)
1467
0
      goto error;
1468
4
    
if (4
isl_basic_set_plain_is_universe(t)4
) {
1469
0
      isl_basic_set_list_free(list);
1470
0
      return t;
1471
0
    }
1472
4
    
if (4
t->n_eq < isl_basic_set_total_dim(t)4
) {
1473
2
      set = isl_basic_set_list_union(list);
1474
2
      return modulo_lineality(set, t);
1475
2
    }
1476
2
    isl_basic_set_free(t);
1477
2
  }
1478
24
1479
0
  return NULL;
1480
0
error:
1481
0
  isl_basic_set_list_free(list);
1482
0
  return NULL;
1483
24
}
1484
1485
/* Compute an initial hull for wrapping containing a single initial
1486
 * facet.
1487
 * This function assumes that the given set is bounded.
1488
 */
1489
static __isl_give isl_basic_set *initial_hull(__isl_take isl_basic_set *hull,
1490
  __isl_keep isl_set *set)
1491
3.88k
{
1492
3.88k
  struct isl_mat *bounds = NULL;
1493
3.88k
  unsigned dim;
1494
3.88k
  int k;
1495
3.88k
1496
3.88k
  if (!hull)
1497
0
    goto error;
1498
3.88k
  bounds = initial_facet_constraint(set);
1499
3.88k
  if (!bounds)
1500
0
    goto error;
1501
3.88k
  k = isl_basic_set_alloc_inequality(hull);
1502
3.88k
  if (k < 0)
1503
0
    goto error;
1504
3.88k
  dim = isl_set_n_dim(set);
1505
3.88k
  isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error);
1506
3.88k
  isl_seq_cpy(hull->ineq[k], bounds->row[0], bounds->n_col);
1507
3.88k
  isl_mat_free(bounds);
1508
3.88k
1509
3.88k
  return hull;
1510
0
error:
1511
0
  isl_basic_set_free(hull);
1512
0
  isl_mat_free(bounds);
1513
0
  return NULL;
1514
3.88k
}
1515
1516
struct max_constraint {
1517
  struct isl_mat *c;
1518
  int   count;
1519
  int   ineq;
1520
};
1521
1522
static int max_constraint_equal(const void *entry, const void *val)
1523
37
{
1524
37
  struct max_constraint *a = (struct max_constraint *)entry;
1525
37
  isl_int *b = (isl_int *)val;
1526
37
1527
37
  return isl_seq_eq(a->c->row[0] + 1, b, a->c->n_col - 1);
1528
37
}
1529
1530
static void update_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1531
  isl_int *con, unsigned len, int n, int ineq)
1532
3.08k
{
1533
3.08k
  struct isl_hash_table_entry *entry;
1534
3.08k
  struct max_constraint *c;
1535
3.08k
  uint32_t c_hash;
1536
3.08k
1537
3.08k
  c_hash = isl_seq_get_hash(con + 1, len);
1538
3.08k
  entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1539
3.08k
      con + 1, 0);
1540
3.08k
  if (!entry)
1541
3.05k
    return;
1542
37
  c = entry->data;
1543
37
  if (
c->count < n37
) {
1544
0
    isl_hash_table_remove(ctx, table, entry);
1545
0
    return;
1546
0
  }
1547
37
  c->count++;
1548
37
  if (isl_int_gt(c->c->row[0][0], con[0]))
1549
3
    return;
1550
34
  
if (34
isl_int_eq34
(c->c->row[0][0], con[0])) {
1551
34
    if (ineq)
1552
16
      c->ineq = ineq;
1553
34
    return;
1554
34
  }
1555
0
  c->c = isl_mat_cow(c->c);
1556
0
  isl_int_set(c->c->row[0][0], con[0]);
1557
3.08k
  c->ineq = ineq;
1558
3.08k
}
1559
1560
/* Check whether the constraint hash table "table" constains the constraint
1561
 * "con".
1562
 */
1563
static int has_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
1564
  isl_int *con, unsigned len, int n)
1565
0
{
1566
0
  struct isl_hash_table_entry *entry;
1567
0
  struct max_constraint *c;
1568
0
  uint32_t c_hash;
1569
0
1570
0
  c_hash = isl_seq_get_hash(con + 1, len);
1571
0
  entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
1572
0
      con + 1, 0);
1573
0
  if (!entry)
1574
0
    return 0;
1575
0
  c = entry->data;
1576
0
  if (c->count < n)
1577
0
    return 0;
1578
0
  
return 0
isl_int_eq0
(c->c->row[0][0], con[0]);
1579
0
}
1580
1581
/* Check for inequality constraints of a basic set without equalities
1582
 * such that the same or more stringent copies of the constraint appear
1583
 * in all of the basic sets.  Such constraints are necessarily facet
1584
 * constraints of the convex hull.
1585
 *
1586
 * If the resulting basic set is by chance identical to one of
1587
 * the basic sets in "set", then we know that this basic set contains
1588
 * all other basic sets and is therefore the convex hull of set.
1589
 * In this case we set *is_hull to 1.
1590
 */
1591
static __isl_give isl_basic_set *common_constraints(
1592
  __isl_take isl_basic_set *hull, __isl_keep isl_set *set, int *is_hull)
1593
3.90k
{
1594
3.90k
  int i, j, s, n;
1595
3.90k
  int min_constraints;
1596
3.90k
  int best;
1597
3.90k
  struct max_constraint *constraints = NULL;
1598
3.90k
  struct isl_hash_table *table = NULL;
1599
3.90k
  unsigned total;
1600
3.90k
1601
3.90k
  *is_hull = 0;
1602
3.90k
1603
11.0k
  for (i = 0; 
i < set->n11.0k
;
++i7.10k
)
1604
7.77k
    
if (7.77k
set->p[i]->n_eq == 07.77k
)
1605
674
      break;
1606
3.90k
  if (i >= set->n)
1607
3.23k
    return hull;
1608
674
  min_constraints = set->p[i]->n_ineq;
1609
674
  best = i;
1610
710
  for (i = best + 1; 
i < set->n710
;
++i36
) {
1611
36
    if (set->p[i]->n_eq != 0)
1612
8
      continue;
1613
28
    
if (28
set->p[i]->n_ineq >= min_constraints28
)
1614
27
      continue;
1615
1
    min_constraints = set->p[i]->n_ineq;
1616
1
    best = i;
1617
1
  }
1618
674
  constraints = isl_calloc_array(hull->ctx, struct max_constraint,
1619
674
          min_constraints);
1620
674
  if (!constraints)
1621
0
    return hull;
1622
674
  
table = 674
isl_alloc_type674
(hull->ctx, struct isl_hash_table);
1623
674
  if (isl_hash_table_init(hull->ctx, table, min_constraints))
1624
0
    goto error;
1625
674
1626
674
  total = isl_space_dim(set->dim, isl_dim_all);
1627
3.15k
  for (i = 0; 
i < set->p[best]->n_ineq3.15k
;
++i2.48k
) {
1628
2.48k
    constraints[i].c = isl_mat_sub_alloc6(hull->ctx,
1629
2.48k
      set->p[best]->ineq + i, 0, 1, 0, 1 + total);
1630
2.48k
    if (!constraints[i].c)
1631
0
      goto error;
1632
2.48k
    constraints[i].ineq = 1;
1633
2.48k
  }
1634
3.15k
  
for (i = 0; 674
i < min_constraints3.15k
;
++i2.48k
) {
1635
2.48k
    struct isl_hash_table_entry *entry;
1636
2.48k
    uint32_t c_hash;
1637
2.48k
    c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total);
1638
2.48k
    entry = isl_hash_table_find(hull->ctx, table, c_hash,
1639
2.48k
      max_constraint_equal, constraints[i].c->row[0] + 1, 1);
1640
2.48k
    if (!entry)
1641
0
      goto error;
1642
2.48k
    
isl_assert2.48k
(hull->ctx, !entry->data, goto error);
1643
2.48k
    entry->data = &constraints[i];
1644
2.48k
  }
1645
674
1646
674
  n = 0;
1647
2.02k
  for (s = 0; 
s < set->n2.02k
;
++s1.34k
) {
1648
1.34k
    if (s == best)
1649
674
      continue;
1650
674
1651
1.37k
    
for (i = 0; 674
i < set->p[s]->n_eq1.37k
;
++i696
) {
1652
696
      isl_int *eq = set->p[s]->eq[i];
1653
2.08k
      for (j = 0; 
j < 22.08k
;
++j1.39k
) {
1654
1.39k
        isl_seq_neg(eq, eq, 1 + total);
1655
1.39k
        update_constraint(hull->ctx, table,
1656
1.39k
                  eq, total, n, 0);
1657
1.39k
      }
1658
696
    }
1659
2.37k
    for (i = 0; 
i < set->p[s]->n_ineq2.37k
;
++i1.69k
) {
1660
1.69k
      isl_int *ineq = set->p[s]->ineq[i];
1661
1.69k
      update_constraint(hull->ctx, table, ineq, total, n,
1662
1.69k
        set->p[s]->n_eq == 0);
1663
1.69k
    }
1664
1.34k
    ++n;
1665
1.34k
  }
1666
674
1667
3.15k
  for (i = 0; 
i < min_constraints3.15k
;
++i2.48k
) {
1668
2.48k
    if (constraints[i].count < n)
1669
2.44k
      continue;
1670
37
    
if (37
!constraints[i].ineq37
)
1671
0
      continue;
1672
37
    j = isl_basic_set_alloc_inequality(hull);
1673
37
    if (j < 0)
1674
0
      goto error;
1675
37
    isl_seq_cpy(hull->ineq[j], constraints[i].c->row[0], 1 + total);
1676
37
  }
1677
674
1678
2.02k
  
for (s = 0; 674
s < set->n2.02k
;
++s1.34k
) {
1679
1.34k
    if (set->p[s]->n_eq)
1680
646
      continue;
1681
702
    
if (702
set->p[s]->n_ineq != hull->n_ineq702
)
1682
702
      continue;
1683
0
    
for (i = 0; 0
i < set->p[s]->n_ineq0
;
++i0
) {
1684
0
      isl_int *ineq = set->p[s]->ineq[i];
1685
0
      if (!has_constraint(hull->ctx, table, ineq, total, n))
1686
0
        break;
1687
0
    }
1688
0
    if (i == set->p[s]->n_ineq)
1689
0
      *is_hull = 1;
1690
1.34k
  }
1691
674
1692
674
  isl_hash_table_clear(table);
1693
3.15k
  for (i = 0; 
i < min_constraints3.15k
;
++i2.48k
)
1694
2.48k
    isl_mat_free(constraints[i].c);
1695
674
  free(constraints);
1696
674
  free(table);
1697
674
  return hull;
1698
0
error:
1699
0
  isl_hash_table_clear(table);
1700
0
  free(table);
1701
0
  if (constraints)
1702
0
    
for (i = 0; 0
i < min_constraints0
;
++i0
)
1703
0
      isl_mat_free(constraints[i].c);
1704
0
  free(constraints);
1705
0
  return hull;
1706
3.90k
}
1707
1708
/* Create a template for the convex hull of "set" and fill it up
1709
 * obvious facet constraints, if any.  If the result happens to
1710
 * be the convex hull of "set" then *is_hull is set to 1.
1711
 */
1712
static __isl_give isl_basic_set *proto_hull(__isl_keep isl_set *set,
1713
  int *is_hull)
1714
3.90k
{
1715
3.90k
  struct isl_basic_set *hull;
1716
3.90k
  unsigned n_ineq;
1717
3.90k
  int i;
1718
3.90k
1719
3.90k
  n_ineq = 1;
1720
11.7k
  for (i = 0; 
i < set->n11.7k
;
++i7.81k
) {
1721
7.81k
    n_ineq += set->p[i]->n_eq;
1722
7.81k
    n_ineq += set->p[i]->n_ineq;
1723
7.81k
  }
1724
3.90k
  hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
1725
3.90k
  hull = isl_basic_set_set_rational(hull);
1726
3.90k
  if (!hull)
1727
0
    return NULL;
1728
3.90k
  return common_constraints(hull, set, is_hull);
1729
3.90k
}
1730
1731
static __isl_give isl_basic_set *uset_convex_hull_wrap(__isl_take isl_set *set)
1732
3.90k
{
1733
3.90k
  struct isl_basic_set *hull;
1734
3.90k
  int is_hull;
1735
3.90k
1736
3.90k
  hull = proto_hull(set, &is_hull);
1737
3.90k
  if (
hull && 3.90k
!is_hull3.90k
) {
1738
3.90k
    if (hull->n_ineq == 0)
1739
3.88k
      hull = initial_hull(hull, set);
1740
3.90k
    hull = extend(hull, set);
1741
3.90k
  }
1742
3.90k
  isl_set_free(set);
1743
3.90k
1744
3.90k
  return hull;
1745
3.90k
}
1746
1747
/* Compute the convex hull of a set without any parameters or
1748
 * integer divisions.  Depending on whether the set is bounded,
1749
 * we pass control to the wrapping based convex hull or
1750
 * the Fourier-Motzkin elimination based convex hull.
1751
 * We also handle a few special cases before checking the boundedness.
1752
 */
1753
static __isl_give isl_basic_set *uset_convex_hull(__isl_take isl_set *set)
1754
59
{
1755
59
  isl_bool bounded;
1756
59
  struct isl_basic_set *convex_hull = NULL;
1757
59
  struct isl_basic_set *lin;
1758
59
1759
59
  if (isl_set_n_dim(set) == 0)
1760
0
    return convex_hull_0d(set);
1761
59
1762
59
  set = isl_set_coalesce(set);
1763
59
  set = isl_set_set_rational(set);
1764
59
1765
59
  if (!set)
1766
0
    return NULL;
1767
59
  
if (59
set->n == 159
) {
1768
14
    convex_hull = isl_basic_set_copy(set->p[0]);
1769
14
    isl_set_free(set);
1770
14
    return convex_hull;
1771
14
  }
1772
45
  
if (45
isl_set_n_dim(set) == 145
)
1773
2
    return convex_hull_1d(set);
1774
43
1775
43
  bounded = isl_set_is_bounded(set);
1776
43
  if (bounded < 0)
1777
0
    goto error;
1778
43
  
if (43
bounded && 43
set->ctx->opt->convex == 14
ISL_CONVEX_HULL_WRAP14
)
1779
11
    return uset_convex_hull_wrap(set);
1780
32
1781
32
  lin = isl_set_combined_lineality_space(isl_set_copy(set));
1782
32
  if (!lin)
1783
0
    goto error;
1784
32
  
if (32
isl_basic_set_plain_is_universe(lin)32
) {
1785
0
    isl_set_free(set);
1786
0
    return lin;
1787
0
  }
1788
32
  
if (32
lin->n_eq < isl_basic_set_total_dim(lin)32
)
1789
8
    return modulo_lineality(set, lin);
1790
24
  isl_basic_set_free(lin);
1791
24
1792
24
  return uset_convex_hull_unbounded(set);
1793
0
error:
1794
0
  isl_set_free(set);
1795
0
  isl_basic_set_free(convex_hull);
1796
0
  return NULL;
1797
59
}
1798
1799
/* This is the core procedure, where "set" is a "pure" set, i.e.,
1800
 * without parameters or divs and where the convex hull of set is
1801
 * known to be full-dimensional.
1802
 */
1803
static __isl_give isl_basic_set *uset_convex_hull_wrap_bounded(
1804
  __isl_take isl_set *set)
1805
13.5k
{
1806
13.5k
  struct isl_basic_set *convex_hull = NULL;
1807
13.5k
1808
13.5k
  if (!set)
1809
0
    goto error;
1810
13.5k
1811
13.5k
  
if (13.5k
isl_set_n_dim(set) == 013.5k
) {
1812
0
    convex_hull = isl_basic_set_universe(isl_space_copy(set->dim));
1813
0
    isl_set_free(set);
1814
0
    convex_hull = isl_basic_set_set_rational(convex_hull);
1815
0
    return convex_hull;
1816
0
  }
1817
13.5k
1818
13.5k
  set = isl_set_set_rational(set);
1819
13.5k
  set = isl_set_coalesce(set);
1820
13.5k
  if (!set)
1821
0
    goto error;
1822
13.5k
  
if (13.5k
set->n == 113.5k
) {
1823
5.61k
    convex_hull = isl_basic_set_copy(set->p[0]);
1824
5.61k
    isl_set_free(set);
1825
5.61k
    convex_hull = isl_basic_map_remove_redundancies(convex_hull);
1826
5.61k
    return convex_hull;
1827
5.61k
  }
1828
7.96k
  
if (7.96k
isl_set_n_dim(set) == 17.96k
)
1829
4.06k
    return convex_hull_1d(set);
1830
3.89k
1831
3.89k
  return uset_convex_hull_wrap(set);
1832
0
error:
1833
0
  isl_set_free(set);
1834
0
  return NULL;
1835
13.5k
}
1836
1837
/* Compute the convex hull of set "set" with affine hull "affine_hull",
1838
 * We first remove the equalities (transforming the set), compute the
1839
 * convex hull of the transformed set and then add the equalities back
1840
 * (after performing the inverse transformation.
1841
 */
1842
static __isl_give isl_basic_set *modulo_affine_hull(
1843
  __isl_take isl_set *set, __isl_take isl_basic_set *affine_hull)
1844
3
{
1845
3
  struct isl_mat *T;
1846
3
  struct isl_mat *T2;
1847
3
  struct isl_basic_set *dummy;
1848
3
  struct isl_basic_set *convex_hull;
1849
3
1850
3
  dummy = isl_basic_set_remove_equalities(
1851
3
      isl_basic_set_copy(affine_hull), &T, &T2);
1852
3
  if (!dummy)
1853
0
    goto error;
1854
3
  isl_basic_set_free(dummy);
1855
3
  set = isl_set_preimage(set, T);
1856
3
  convex_hull = uset_convex_hull(set);
1857
3
  convex_hull = isl_basic_set_preimage(convex_hull, T2);
1858
3
  convex_hull = isl_basic_set_intersect(convex_hull, affine_hull);
1859
3
  return convex_hull;
1860
0
error:
1861
0
  isl_mat_free(T);
1862
0
  isl_mat_free(T2);
1863
0
  isl_basic_set_free(affine_hull);
1864
0
  isl_set_free(set);
1865
0
  return NULL;
1866
3
}
1867
1868
/* Return an empty basic map living in the same space as "map".
1869
 */
1870
static __isl_give isl_basic_map *replace_map_by_empty_basic_map(
1871
  __isl_take isl_map *map)
1872
796
{
1873
796
  isl_space *space;
1874
796
1875
796
  space = isl_map_get_space(map);
1876
796
  isl_map_free(map);
1877
796
  return isl_basic_map_empty(space);
1878
796
}
1879
1880
/* Compute the convex hull of a map.
1881
 *
1882
 * The implementation was inspired by "Extended Convex Hull" by Fukuda et al.,
1883
 * specifically, the wrapping of facets to obtain new facets.
1884
 */
1885
__isl_give isl_basic_map *isl_map_convex_hull(__isl_take isl_map *map)
1886
40
{
1887
40
  struct isl_basic_set *bset;
1888
40
  struct isl_basic_map *model = NULL;
1889
40
  struct isl_basic_set *affine_hull = NULL;
1890
40
  struct isl_basic_map *convex_hull = NULL;
1891
40
  struct isl_set *set = NULL;
1892
40
1893
40
  map = isl_map_detect_equalities(map);
1894
40
  map = isl_map_align_divs_internal(map);
1895
40
  if (!map)
1896
0
    goto error;
1897
40
1898
40
  
if (40
map->n == 040
)
1899
2
    return replace_map_by_empty_basic_map(map);
1900
38
1901
38
  model = isl_basic_map_copy(map->p[0]);
1902
38
  set = isl_map_underlying_set(map);
1903
38
  if (!set)
1904
0
    goto error;
1905
38
1906
38
  affine_hull = isl_set_affine_hull(isl_set_copy(set));
1907
38
  if (!affine_hull)
1908
0
    goto error;
1909
38
  
if (38
affine_hull->n_eq != 038
)
1910
2
    bset = modulo_affine_hull(set, affine_hull);
1911
36
  else {
1912
36
    isl_basic_set_free(affine_hull);
1913
36
    bset = uset_convex_hull(set);
1914
36
  }
1915
38
1916
38
  convex_hull = isl_basic_map_overlying_set(bset, model);
1917
38
  if (!convex_hull)
1918
0
    return NULL;
1919
38
1920
38
  
ISL_F_SET38
(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT);
1921
38
  ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES);
1922
38
  ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL);
1923
38
  return convex_hull;
1924
0
error:
1925
0
  isl_set_free(set);
1926
0
  isl_basic_map_free(model);
1927
0
  return NULL;
1928
40
}
1929
1930
struct isl_basic_set *isl_set_convex_hull(struct isl_set *set)
1931
40
{
1932
40
  return bset_from_bmap(isl_map_convex_hull(set_to_map(set)));
1933
40
}
1934
1935
__isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map)
1936
0
{
1937
0
  isl_basic_map *hull;
1938
0
1939
0
  hull = isl_map_convex_hull(map);
1940
0
  return isl_basic_map_remove_divs(hull);
1941
0
}
1942
1943
__isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set)
1944
0
{
1945
0
  return bset_from_bmap(isl_map_polyhedral_hull(set_to_map(set)));
1946
0
}
1947
1948
struct sh_data_entry {
1949
  struct isl_hash_table *table;
1950
  struct isl_tab    *tab;
1951
};
1952
1953
/* Holds the data needed during the simple hull computation.
1954
 * In particular,
1955
 *  n   the number of basic sets in the original set
1956
 *  hull_table  a hash table of already computed constraints
1957
 *      in the simple hull
1958
 *  p   for each basic set,
1959
 *    table   a hash table of the constraints
1960
 *    tab   the tableau corresponding to the basic set
1961
 */
1962
struct sh_data {
1963
  struct isl_ctx    *ctx;
1964
  unsigned    n;
1965
  struct isl_hash_table *hull_table;
1966
  struct sh_data_entry  p[1];
1967
};
1968
1969
static void sh_data_free(struct sh_data *data)
1970
5.02k
{
1971
5.02k
  int i;
1972
5.02k
1973
5.02k
  if (!data)
1974
0
    return;
1975
5.02k
  isl_hash_table_free(data->ctx, data->hull_table);
1976
15.9k
  for (i = 0; 
i < data->n15.9k
;
++i10.9k
) {
1977
10.9k
    isl_hash_table_free(data->ctx, data->p[i].table);
1978
10.9k
    isl_tab_free(data->p[i].tab);
1979
10.9k
  }
1980
5.02k
  free(data);
1981
5.02k
}
1982
1983
struct ineq_cmp_data {
1984
  unsigned  len;
1985
  isl_int   *p;
1986
};
1987
1988
static int has_ineq(const void *entry, const void *val)
1989
156k
{
1990
156k
  isl_int *row = (isl_int *)entry;
1991
156k
  struct ineq_cmp_data *v = (struct ineq_cmp_data *)val;
1992
156k
1993
156k
  return isl_seq_eq(row + 1, v->p + 1, v->len) ||
1994
5.50k
         isl_seq_is_neg(row + 1, v->p + 1, v->len);
1995
156k
}
1996
1997
static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,
1998
      isl_int *ineq, unsigned len)
1999
121k
{
2000
121k
  uint32_t c_hash;
2001
121k
  struct ineq_cmp_data v;
2002
121k
  struct isl_hash_table_entry *entry;
2003
121k
2004
121k
  v.len = len;
2005
121k
  v.p = ineq;
2006
121k
  c_hash = isl_seq_get_hash(ineq + 1, len);
2007
121k
  entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1);
2008
121k
  if (!entry)
2009
0
    return - 1;
2010
121k
  entry->data = ineq;
2011
121k
  return 0;
2012
121k
}
2013
2014
/* Fill hash table "table" with the constraints of "bset".
2015
 * Equalities are added as two inequalities.
2016
 * The value in the hash table is a pointer to the (in)equality of "bset".
2017
 */
2018
static int hash_basic_set(struct isl_hash_table *table,
2019
  __isl_keep isl_basic_set *bset)
2020
10.9k
{
2021
10.9k
  int i, j;
2022
10.9k
  unsigned dim = isl_basic_set_total_dim(bset);
2023
10.9k
2024
16.6k
  for (i = 0; 
i < bset->n_eq16.6k
;
++i5.63k
) {
2025
16.8k
    for (j = 0; 
j < 216.8k
;
++j11.2k
) {
2026
11.2k
      isl_seq_neg(bset->eq[i], bset->eq[i], 1 + dim);
2027
11.2k
      if (hash_ineq(bset->ctx, table, bset->eq[i], dim) < 0)
2028
0
        return -1;
2029
11.2k
    }
2030
5.63k
  }
2031
120k
  
for (i = 0; 10.9k
i < bset->n_ineq120k
;
++i109k
) {
2032
109k
    if (hash_ineq(bset->ctx, table, bset->ineq[i], dim) < 0)
2033
0
      return -1;
2034
109k
  }
2035
10.9k
  return 0;
2036
10.9k
}
2037
2038
static struct sh_data *sh_data_alloc(__isl_keep isl_set *set, unsigned n_ineq)
2039
5.02k
{
2040
5.02k
  struct sh_data *data;
2041
5.02k
  int i;
2042
5.02k
2043
5.02k
  data = isl_calloc(set->ctx, struct sh_data,
2044
5.02k
    sizeof(struct sh_data) +
2045
5.02k
    (set->n - 1) * sizeof(struct sh_data_entry));
2046
5.02k
  if (!data)
2047
0
    return NULL;
2048
5.02k
  data->ctx = set->ctx;
2049
5.02k
  data->n = set->n;
2050
5.02k
  data->hull_table = isl_hash_table_alloc(set->ctx, n_ineq);
2051
5.02k
  if (!data->hull_table)
2052
0
    goto error;
2053
15.9k
  
for (i = 0; 5.02k
i < set->n15.9k
;
++i10.9k
) {
2054
10.9k
    data->p[i].table = isl_hash_table_alloc(set->ctx,
2055
10.9k
            2 * set->p[i]->n_eq + set->p[i]->n_ineq);
2056
10.9k
    if (!data->p[i].table)
2057
0
      goto error;
2058
10.9k
    
if (10.9k
hash_basic_set(data->p[i].table, set->p[i]) < 010.9k
)
2059
0
      goto error;
2060
10.9k
  }
2061
5.02k
  return data;
2062
0
error:
2063
0
  sh_data_free(data);
2064
0
  return NULL;
2065
5.02k
}
2066
2067
/* Check if inequality "ineq" is a bound for basic set "j" or if
2068
 * it can be relaxed (by increasing the constant term) to become
2069
 * a bound for that basic set.  In the latter case, the constant
2070
 * term is updated.
2071
 * Relaxation of the constant term is only allowed if "shift" is set.
2072
 *
2073
 * Return 1 if "ineq" is a bound
2074
 *    0 if "ineq" may attain arbitrarily small values on basic set "j"
2075
 *   -1 if some error occurred
2076
 */
2077
static int is_bound(struct sh_data *data, __isl_keep isl_set *set, int j,
2078
  isl_int *ineq, int shift)
2079
8.17k
{
2080
8.17k
  enum isl_lp_result res;
2081
8.17k
  isl_int opt;
2082
8.17k
2083
8.17k
  if (
!data->p[j].tab8.17k
) {
2084
6.00k
    data->p[j].tab = isl_tab_from_basic_set(set->p[j], 0);
2085
6.00k
    if (!data->p[j].tab)
2086
0
      return -1;
2087
8.17k
  }
2088
8.17k
2089
8.17k
  
isl_int_init8.17k
(opt);
2090
8.17k
2091
8.17k
  res = isl_tab_min(data->p[j].tab, ineq, data->ctx->one,
2092
8.17k
        &opt, NULL, 0);
2093
8.17k
  if (
res == isl_lp_ok && 8.17k
isl_int_is_neg1.54k
(opt)) {
2094
667
    if (shift)
2095
41
      isl_int_sub(ineq[0], ineq[0], opt);
2096
667
    else
2097
626
      res = isl_lp_unbounded;
2098
667
  }
2099
8.17k
2100
8.17k
  isl_int_clear(opt);
2101
8.17k
2102
8.17k
  return (res == isl_lp_ok || 
res == isl_lp_empty7.25k
) ?
1921
:
2103
7.25k
         
res == isl_lp_unbounded ? 7.25k
07.25k
:
-10
;
2104
8.17k
}
2105
2106
/* Set the constant term of "ineq" to the maximum of those of the constraints
2107
 * in the basic sets of "set" following "i" that are parallel to "ineq".
2108
 * That is, if any of the basic sets of "set" following "i" have a more
2109
 * relaxed copy of "ineq", then replace "ineq" by the most relaxed copy.
2110
 * "c_hash" is the hash value of the linear part of "ineq".
2111
 * "v" has been set up for use by has_ineq.
2112
 *
2113
 * Note that the two inequality constraints corresponding to an equality are
2114
 * represented by the same inequality constraint in data->p[j].table
2115
 * (but with different hash values).  This means the constraint (or at
2116
 * least its constant term) may need to be temporarily negated to get
2117
 * the actually hashed constraint.
2118
 */
2119
static void set_max_constant_term(struct sh_data *data, __isl_keep isl_set *set,
2120
  int i, isl_int *ineq, uint32_t c_hash, struct ineq_cmp_data *v)
2121
42.7k
{
2122
42.7k
  int j;
2123
42.7k
  isl_ctx *ctx;
2124
42.7k
  struct isl_hash_table_entry *entry;
2125
42.7k
2126
42.7k
  ctx = isl_set_get_ctx(set);
2127
83.7k
  for (j = i + 1; 
j < set->n83.7k
;
++j41.0k
) {
2128
41.0k
    int neg;
2129
41.0k
    isl_int *ineq_j;
2130
41.0k
2131
41.0k
    entry = isl_hash_table_find(ctx, data->p[j].table,
2132
41.0k
            c_hash, &has_ineq, v, 0);
2133
41.0k
    if (!entry)
2134
2.71k
      continue;
2135
38.2k
2136
38.2k
    ineq_j = entry->data;
2137
38.2k
    neg = isl_seq_is_neg(ineq_j + 1, ineq + 1, v->len);
2138
38.2k
    if (neg)
2139
2.66k
      isl_int_neg(ineq_j[0], ineq_j[0]);
2140
38.2k
    if (isl_int_gt(ineq_j[0], ineq[0]))
2141
355
      isl_int_set(ineq[0], ineq_j[0]);
2142
38.2k
    if (neg)
2143
2.66k
      isl_int_neg(ineq_j[0], ineq_j[0]);
2144
41.0k
  }
2145
42.7k
}
2146
2147
/* Check if inequality "ineq" from basic set "i" is or can be relaxed to
2148
 * become a bound on the whole set.  If so, add the (relaxed) inequality
2149
 * to "hull".  Relaxation is only allowed if "shift" is set.
2150
 *
2151
 * We first check if "hull" already contains a translate of the inequality.
2152
 * If so, we are done.
2153
 * Then, we check if any of the previous basic sets contains a translate
2154
 * of the inequality.  If so, then we have already considered this
2155
 * inequality and we are done.
2156
 * Otherwise, for each basic set other than "i", we check if the inequality
2157
 * is a bound on the basic set, but first replace the constant term
2158
 * by the maximal value of any translate of the inequality in any
2159
 * of the following basic sets.
2160
 * For previous basic sets, we know that they do not contain a translate
2161
 * of the inequality, so we directly call is_bound.
2162
 * For following basic sets, we first check if a translate of the
2163
 * inequality appears in its description.  If so, the constant term
2164
 * of the inequality has already been updated with respect to this
2165
 * translate and the inequality is therefore known to be a bound
2166
 * of this basic set.
2167
 */
2168
static __isl_give isl_basic_set *add_bound(__isl_take isl_basic_set *hull,
2169
  struct sh_data *data, __isl_keep isl_set *set, int i, isl_int *ineq,
2170
  int shift)
2171
81.0k
{
2172
81.0k
  uint32_t c_hash;
2173
81.0k
  struct ineq_cmp_data v;
2174
81.0k
  struct isl_hash_table_entry *entry;
2175
81.0k
  int j, k;
2176
81.0k
2177
81.0k
  if (!hull)
2178
0
    return NULL;
2179
81.0k
2180
81.0k
  v.len = isl_basic_set_total_dim(hull);
2181
81.0k
  v.p = ineq;
2182
81.0k
  c_hash = isl_seq_get_hash(ineq + 1, v.len);
2183
81.0k
2184
81.0k
  entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2185
81.0k
          has_ineq, &v, 0);
2186
81.0k
  if (entry)
2187
38.2k
    return hull;
2188
42.7k
2189
45.5k
  
for (j = 0; 42.7k
j < i45.5k
;
++j2.75k
) {
2190
2.84k
    entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2191
2.84k
            c_hash, has_ineq, &v, 0);
2192
2.84k
    if (entry)
2193
91
      break;
2194
2.84k
  }
2195
42.7k
  if (j < i)
2196
91
    return hull;
2197
42.7k
2198
42.7k
  k = isl_basic_set_alloc_inequality(hull);
2199
42.7k
  if (k < 0)
2200
0
    goto error;
2201
42.7k
  isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2202
42.7k
2203
42.7k
  set_max_constant_term(data, set, i, hull->ineq[k], c_hash, &v);
2204
42.8k
  for (j = 0; 
j < i42.8k
;
++j124
) {
2205
2.58k
    int bound;
2206
2.58k
    bound = is_bound(data, set, j, hull->ineq[k], shift);
2207
2.58k
    if (bound < 0)
2208
0
      goto error;
2209
2.58k
    
if (2.58k
!bound2.58k
)
2210
2.45k
      break;
2211
2.58k
  }
2212
42.7k
  
if (42.7k
j < i42.7k
) {
2213
2.45k
    isl_basic_set_free_inequality(hull, 1);
2214
2.45k
    return hull;
2215
2.45k
  }
2216
40.2k
2217
78.6k
  
for (j = i + 1; 40.2k
j < set->n78.6k
;
++j38.3k
) {
2218
40.7k
    int bound;
2219
40.7k
    entry = isl_hash_table_find(hull->ctx, data->p[j].table,
2220
40.7k
            c_hash, has_ineq, &v, 0);
2221
40.7k
    if (entry)
2222
38.2k
      continue;
2223
2.57k
    bound = is_bound(data, set, j, hull->ineq[k], shift);
2224
2.57k
    if (bound < 0)
2225
0
      goto error;
2226
2.57k
    
if (2.57k
!bound2.57k
)
2227
2.41k
      break;
2228
40.7k
  }
2229
40.2k
  
if (40.2k
j < set->n40.2k
) {
2230
2.41k
    isl_basic_set_free_inequality(hull, 1);
2231
2.41k
    return hull;
2232
2.41k
  }
2233
37.8k
2234
37.8k
  entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
2235
37.8k
          has_ineq, &v, 1);
2236
37.8k
  if (!entry)
2237
0
    goto error;
2238
37.8k
  entry->data = hull->ineq[k];
2239
37.8k
2240
37.8k
  return hull;
2241
0
error:
2242
0
  isl_basic_set_free(hull);
2243
0
  return NULL;
2244
81.0k
}
2245
2246
/* Check if any inequality from basic set "i" is or can be relaxed to
2247
 * become a bound on the whole set.  If so, add the (relaxed) inequality
2248
 * to "hull".  Relaxation is only allowed if "shift" is set.
2249
 */
2250
static __isl_give isl_basic_set *add_bounds(__isl_take isl_basic_set *bset,
2251
  struct sh_data *data, __isl_keep isl_set *set, int i, int shift)
2252
8.69k
{
2253
8.69k
  int j, k;
2254
8.69k
  unsigned dim = isl_basic_set_total_dim(bset);
2255
8.69k
2256
13.9k
  for (j = 0; 
j < set->p[i]->n_eq13.9k
;
++j5.29k
) {
2257
15.8k
    for (k = 0; 
k < 215.8k
;
++k10.5k
) {
2258
10.5k
      isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim);
2259
10.5k
      bset = add_bound(bset, data, set, i, set->p[i]->eq[j],
2260
10.5k
              shift);
2261
10.5k
    }
2262
5.29k
  }
2263
79.1k
  for (j = 0; 
j < set->p[i]->n_ineq79.1k
;
++j70.4k
)
2264
70.4k
    bset = add_bound(bset, data, set, i, set->p[i]->ineq[j], shift);
2265
8.69k
  return bset;
2266
8.69k
}
2267
2268
/* Compute a superset of the convex hull of set that is described
2269
 * by only (translates of) the constraints in the constituents of set.
2270
 * Translation is only allowed if "shift" is set.
2271
 */
2272
static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set,
2273
  int shift)
2274
4.28k
{
2275
4.28k
  struct sh_data *data = NULL;
2276
4.28k
  struct isl_basic_set *hull = NULL;
2277
4.28k
  unsigned n_ineq;
2278
4.28k
  int i;
2279
4.28k
2280
4.28k
  if (!set)
2281
0
    return NULL;
2282
4.28k
2283
4.28k
  n_ineq = 0;
2284
12.9k
  for (i = 0; 
i < set->n12.9k
;
++i8.69k
) {
2285
8.69k
    if (!set->p[i])
2286
0
      goto error;
2287
8.69k
    n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq;
2288
8.69k
  }
2289
4.28k
2290
4.28k
  hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq);
2291
4.28k
  if (!hull)
2292
0
    goto error;
2293
4.28k
2294
4.28k
  data = sh_data_alloc(set, n_ineq);
2295
4.28k
  if (!data)
2296
0
    goto error;
2297
4.28k
2298
12.9k
  
for (i = 0; 4.28k
i < set->n12.9k
;
++i8.69k
)
2299
8.69k
    hull = add_bounds(hull, data, set, i, shift);
2300
4.28k
2301
4.28k
  sh_data_free(data);
2302
4.28k
  isl_set_free(set);
2303
4.28k
2304
4.28k
  return hull;
2305
0
error:
2306
0
  sh_data_free(data);
2307
0
  isl_basic_set_free(hull);
2308
0
  isl_set_free(set);
2309
0
  return NULL;
2310
4.28k
}
2311
2312
/* Compute a superset of the convex hull of map that is described
2313
 * by only (translates of) the constraints in the constituents of map.
2314
 * Handle trivial cases where map is NULL or contains at most one disjunct.
2315
 */
2316
static __isl_give isl_basic_map *map_simple_hull_trivial(
2317
  __isl_take isl_map *map)
2318
43.3k
{
2319
43.3k
  isl_basic_map *hull;
2320
43.3k
2321
43.3k
  if (!map)
2322
0
    return NULL;
2323
43.3k
  
if (43.3k
map->n == 043.3k
)
2324
794
    return replace_map_by_empty_basic_map(map);
2325
42.5k
2326
42.5k
  hull = isl_basic_map_copy(map->p[0]);
2327
42.5k
  isl_map_free(map);
2328
42.5k
  return hull;
2329
42.5k
}
2330
2331
/* Return a copy of the simple hull cached inside "map".
2332
 * "shift" determines whether to return the cached unshifted or shifted
2333
 * simple hull.
2334
 */
2335
static __isl_give isl_basic_map *cached_simple_hull(__isl_take isl_map *map,
2336
  int shift)
2337
39
{
2338
39
  isl_basic_map *hull;
2339
39
2340
39
  hull = isl_basic_map_copy(map->cached_simple_hull[shift]);
2341
39
  isl_map_free(map);
2342
39
2343
39
  return hull;
2344
39
}
2345
2346
/* Compute a superset of the convex hull of map that is described
2347
 * by only (translates of) the constraints in the constituents of map.
2348
 * Translation is only allowed if "shift" is set.
2349
 *
2350
 * The constraints are sorted while removing redundant constraints
2351
 * in order to indicate a preference of which constraints should
2352
 * be preserved.  In particular, pairs of constraints that are
2353
 * sorted together are preferred to either both be preserved
2354
 * or both be removed.  The sorting is performed inside
2355
 * isl_basic_map_remove_redundancies.
2356
 *
2357
 * The result of the computation is stored in map->cached_simple_hull[shift]
2358
 * such that it can be reused in subsequent calls.  The cache is cleared
2359
 * whenever the map is modified (in isl_map_cow).
2360
 * Note that the results need to be stored in the input map for there
2361
 * to be any chance that they may get reused.  In particular, they
2362
 * are stored in a copy of the input map that is saved before
2363
 * the integer division alignment.
2364
 */
2365
static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map,
2366
  int shift)
2367
47.6k
{
2368
47.6k
  struct isl_set *set = NULL;
2369
47.6k
  struct isl_basic_map *model = NULL;
2370
47.6k
  struct isl_basic_map *hull;
2371
47.6k
  struct isl_basic_map *affine_hull;
2372
47.6k
  struct isl_basic_set *bset = NULL;
2373
47.6k
  isl_map *input;
2374
47.6k
2375
47.6k
  if (
!map || 47.6k
map->n <= 147.6k
)
2376
42.5k
    return map_simple_hull_trivial(map);
2377
5.11k
2378
5.11k
  
if (5.11k
map->cached_simple_hull[shift]5.11k
)
2379
39
    return cached_simple_hull(map, shift);
2380
5.07k
2381
5.07k
  map = isl_map_detect_equalities(map);
2382
5.07k
  if (
!map || 5.07k
map->n <= 15.07k
)
2383
792
    return map_simple_hull_trivial(map);
2384
4.28k
  affine_hull = isl_map_affine_hull(isl_map_copy(map));
2385
4.28k
  input = isl_map_copy(map);
2386
4.28k
  map = isl_map_align_divs_internal(map);
2387
4.28k
  model = map ? isl_basic_map_copy(map->p[0]) : NULL;
2388
4.28k
2389
4.28k
  set = isl_map_underlying_set(map);
2390
4.28k
2391
4.28k
  bset = uset_simple_hull(set, shift);
2392
4.28k
2393
4.28k
  hull = isl_basic_map_overlying_set(bset, model);
2394
4.28k
2395
4.28k
  hull = isl_basic_map_intersect(hull, affine_hull);
2396
4.28k
  hull = isl_basic_map_remove_redundancies(hull);
2397
4.28k
2398
4.28k
  if (
hull4.28k
) {
2399
4.28k
    ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT);
2400
4.28k
    ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES);
2401
4.28k
  }
2402
4.28k
2403
4.28k
  hull = isl_basic_map_finalize(hull);
2404
4.28k
  if (input)
2405
4.28k
    input->cached_simple_hull[shift] = isl_basic_map_copy(hull);
2406
47.6k
  isl_map_free(input);
2407
47.6k
2408
47.6k
  return hull;
2409
47.6k
}
2410
2411
/* Compute a superset of the convex hull of map that is described
2412
 * by only translates of the constraints in the constituents of map.
2413
 */
2414
__isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map)
2415
31.0k
{
2416
31.0k
  return map_simple_hull(map, 1);
2417
31.0k
}
2418
2419
struct isl_basic_set *isl_set_simple_hull(struct isl_set *set)
2420
14.4k
{
2421
14.4k
  return bset_from_bmap(isl_map_simple_hull(set_to_map(set)));
2422
14.4k
}
2423
2424
/* Compute a superset of the convex hull of map that is described
2425
 * by only the constraints in the constituents of map.
2426
 */
2427
__isl_give isl_basic_map *isl_map_unshifted_simple_hull(
2428
  __isl_take isl_map *map)
2429
16.6k
{
2430
16.6k
  return map_simple_hull(map, 0);
2431
16.6k
}
2432
2433
__isl_give isl_basic_set *isl_set_unshifted_simple_hull(
2434
  __isl_take isl_set *set)
2435
10.0k
{
2436
10.0k
  return isl_map_unshifted_simple_hull(set);
2437
10.0k
}
2438
2439
/* Drop all inequalities from "bmap1" that do not also appear in "bmap2".
2440
 * A constraint that appears with different constant terms
2441
 * in "bmap1" and "bmap2" is also kept, with the least restrictive
2442
 * (i.e., greatest) constant term.
2443
 * "bmap1" and "bmap2" are assumed to have the same (known)
2444
 * integer divisions.
2445
 * The constraints of both "bmap1" and "bmap2" are assumed
2446
 * to have been sorted using isl_basic_map_sort_constraints.
2447
 *
2448
 * Run through the inequality constraints of "bmap1" and "bmap2"
2449
 * in sorted order.
2450
 * Each constraint of "bmap1" without a matching constraint in "bmap2"
2451
 * is removed.
2452
 * If a match is found, the constraint is kept.  If needed, the constant
2453
 * term of the constraint is adjusted.
2454
 */
2455
static __isl_give isl_basic_map *select_shared_inequalities(
2456
  __isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)
2457
652
{
2458
652
  int i1, i2;
2459
652
2460
652
  bmap1 = isl_basic_map_cow(bmap1);
2461
652
  if (
!bmap1 || 652
!bmap2652
)
2462
0
    return isl_basic_map_free(bmap1);
2463
652
2464
652
  i1 = bmap1->n_ineq - 1;
2465
652
  i2 = bmap2->n_ineq - 1;
2466
2.05k
  while (
bmap1 && 2.05k
i1 >= 02.05k
&&
i2 >= 01.56k
) {
2467
1.40k
    int cmp;
2468
1.40k
2469
1.40k
    cmp = isl_basic_map_constraint_cmp(bmap1, bmap1->ineq[i1],
2470
1.40k
              bmap2->ineq[i2]);
2471
1.40k
    if (
cmp < 01.40k
) {
2472
362
      --i2;
2473
362
      continue;
2474
362
    }
2475
1.04k
    
if (1.04k
cmp > 01.04k
) {
2476
450
      if (isl_basic_map_drop_inequality(bmap1, i1) < 0)
2477
0
        bmap1 = isl_basic_map_free(bmap1);
2478
450
      --i1;
2479
450
      continue;
2480
450
    }
2481
593
    
if (593
isl_int_lt593
(bmap1->ineq[i1][0], bmap2->ineq[i2][0]))
2482
29
      isl_int_set(bmap1->ineq[i1][0], bmap2->ineq[i2][0]);
2483
1.40k
    --i1;
2484
1.40k
    --i2;
2485
1.40k
  }
2486
829
  for (; 
i1 >= 0829
;
--i1177
)
2487
177
    
if (177
isl_basic_map_drop_inequality(bmap1, i1) < 0177
)
2488
0
      bmap1 = isl_basic_map_free(bmap1);
2489
652
2490
652
  return bmap1;
2491
652
}
2492
2493
/* Drop all equalities from "bmap1" that do not also appear in "bmap2".
2494
 * "bmap1" and "bmap2" are assumed to have the same (known)
2495
 * integer divisions.
2496
 *
2497
 * Run through the equality constraints of "bmap1" and "bmap2".
2498
 * Each constraint of "bmap1" without a matching constraint in "bmap2"
2499
 * is removed.
2500
 */
2501
static __isl_give isl_basic_map *select_shared_equalities(
2502
  __isl_take isl_basic_map *bmap1, __isl_keep isl_basic_map *bmap2)
2503
652
{
2504
652
  int i1, i2;
2505
652
  unsigned total;
2506
652
2507
652
  bmap1 = isl_basic_map_cow(bmap1);
2508
652
  if (
!bmap1 || 652
!bmap2652
)
2509
0
    return isl_basic_map_free(bmap1);
2510
652
2511
652
  total = isl_basic_map_total_dim(bmap1);
2512
652
2513
652
  i1 = bmap1->n_eq - 1;
2514
652
  i2 = bmap2->n_eq - 1;
2515
683
  while (
bmap1 && 683
i1 >= 0683
&&
i2 >= 038
) {
2516
31
    int last1, last2;
2517
31
2518
31
    last1 = isl_seq_last_non_zero(bmap1->eq[i1] + 1, total);
2519
31
    last2 = isl_seq_last_non_zero(bmap2->eq[i2] + 1, total);
2520
31
    if (
last1 > last231
) {
2521
3
      --i2;
2522
3
      continue;
2523
3
    }
2524
28
    
if (28
last1 < last228
) {
2525
3
      if (isl_basic_map_drop_equality(bmap1, i1) < 0)
2526
0
        bmap1 = isl_basic_map_free(bmap1);
2527
3
      --i1;
2528
3
      continue;
2529
3
    }
2530
25
    
if (25
!isl_seq_eq(bmap1->eq[i1], bmap2->eq[i2], 1 + total)25
) {
2531
3
      if (isl_basic_map_drop_equality(bmap1, i1) < 0)
2532
0
        bmap1 = isl_basic_map_free(bmap1);
2533
3
    }
2534
31
    --i1;
2535
31
    --i2;
2536
31
  }
2537
659
  for (; 
i1 >= 0659
;
--i17
)
2538
7
    
if (7
isl_basic_map_drop_equality(bmap1, i1) < 07
)
2539
0
      bmap1 = isl_basic_map_free(bmap1);
2540
652
2541
652
  return bmap1;
2542
652
}
2543
2544
/* Compute a superset of "bmap1" and "bmap2" that is described
2545
 * by only the constraints that appear in both "bmap1" and "bmap2".
2546
 *
2547
 * First drop constraints that involve unknown integer divisions
2548
 * since it is not trivial to check whether two such integer divisions
2549
 * in different basic maps are the same.
2550
 * Then align the remaining (known) divs and sort the constraints.
2551
 * Finally drop all inequalities and equalities from "bmap1" that
2552
 * do not also appear in "bmap2".
2553
 */
2554
__isl_give isl_basic_map *isl_basic_map_plain_unshifted_simple_hull(
2555
  __isl_take isl_basic_map *bmap1, __isl_take isl_basic_map *bmap2)
2556
652
{
2557
652
  bmap1 = isl_basic_map_drop_constraint_involving_unknown_divs(bmap1);
2558
652
  bmap2 = isl_basic_map_drop_constraint_involving_unknown_divs(bmap2);
2559
652
  bmap2 = isl_basic_map_align_divs(bmap2, bmap1);
2560
652
  bmap1 = isl_basic_map_align_divs(bmap1, bmap2);
2561
652
  bmap1 = isl_basic_map_gauss(bmap1, NULL);
2562
652
  bmap2 = isl_basic_map_gauss(bmap2, NULL);
2563
652
  bmap1 = isl_basic_map_sort_constraints(bmap1);
2564
652
  bmap2 = isl_basic_map_sort_constraints(bmap2);
2565
652
2566
652
  bmap1 = select_shared_inequalities(bmap1, bmap2);
2567
652
  bmap1 = select_shared_equalities(bmap1, bmap2);
2568
652
2569
652
  isl_basic_map_free(bmap2);
2570
652
  bmap1 = isl_basic_map_finalize(bmap1);
2571
652
  return bmap1;
2572
652
}
2573
2574
/* Compute a superset of the convex hull of "map" that is described
2575
 * by only the constraints in the constituents of "map".
2576
 * In particular, the result is composed of constraints that appear
2577
 * in each of the basic maps of "map"
2578
 *
2579
 * Constraints that involve unknown integer divisions are dropped
2580
 * since it is not trivial to check whether two such integer divisions
2581
 * in different basic maps are the same.
2582
 *
2583
 * The hull is initialized from the first basic map and then
2584
 * updated with respect to the other basic maps in turn.
2585
 */
2586
__isl_give isl_basic_map *isl_map_plain_unshifted_simple_hull(
2587
  __isl_take isl_map *map)
2588
496
{
2589
496
  int i;
2590
496
  isl_basic_map *hull;
2591
496
2592
496
  if (!map)
2593
0
    return NULL;
2594
496
  
if (496
map->n <= 1496
)
2595
0
    return map_simple_hull_trivial(map);
2596
496
  map = isl_map_drop_constraint_involving_unknown_divs(map);
2597
496
  hull = isl_basic_map_copy(map->p[0]);
2598
1.14k
  for (i = 1; 
i < map->n1.14k
;
++i652
) {
2599
652
    isl_basic_map *bmap_i;
2600
652
2601
652
    bmap_i = isl_basic_map_copy(map->p[i]);
2602
652
    hull = isl_basic_map_plain_unshifted_simple_hull(hull, bmap_i);
2603
652
  }
2604
496
2605
496
  isl_map_free(map);
2606
496
  return hull;
2607
496
}
2608
2609
/* Compute a superset of the convex hull of "set" that is described
2610
 * by only the constraints in the constituents of "set".
2611
 * In particular, the result is composed of constraints that appear
2612
 * in each of the basic sets of "set"
2613
 */
2614
__isl_give isl_basic_set *isl_set_plain_unshifted_simple_hull(
2615
  __isl_take isl_set *set)
2616
43
{
2617
43
  return isl_map_plain_unshifted_simple_hull(set);
2618
43
}
2619
2620
/* Check if "ineq" is a bound on "set" and, if so, add it to "hull".
2621
 *
2622
 * For each basic set in "set", we first check if the basic set
2623
 * contains a translate of "ineq".  If this translate is more relaxed,
2624
 * then we assume that "ineq" is not a bound on this basic set.
2625
 * Otherwise, we know that it is a bound.
2626
 * If the basic set does not contain a translate of "ineq", then
2627
 * we call is_bound to perform the test.
2628
 */
2629
static __isl_give isl_basic_set *add_bound_from_constraint(
2630
  __isl_take isl_basic_set *hull, struct sh_data *data,
2631
  __isl_keep isl_set *set, isl_int *ineq)
2632
12.6k
{
2633
12.6k
  int i, k;
2634
12.6k
  isl_ctx *ctx;
2635
12.6k
  uint32_t c_hash;
2636
12.6k
  struct ineq_cmp_data v;
2637
12.6k
2638
12.6k
  if (
!hull || 12.6k
!set12.6k
)
2639
0
    return isl_basic_set_free(hull);
2640
12.6k
2641
12.6k
  v.len = isl_basic_set_total_dim(hull);
2642
12.6k
  v.p = ineq;
2643
12.6k
  c_hash = isl_seq_get_hash(ineq + 1, v.len);
2644
12.6k
2645
12.6k
  ctx = isl_basic_set_get_ctx(hull);
2646
52.8k
  for (i = 0; 
i < set->n52.8k
;
++i40.1k
) {
2647
44.5k
    int bound;
2648
44.5k
    struct isl_hash_table_entry *entry;
2649
44.5k
2650
44.5k
    entry = isl_hash_table_find(ctx, data->p[i].table,
2651
44.5k
            c_hash, &has_ineq, &v, 0);
2652
44.5k
    if (
entry44.5k
) {
2653
41.5k
      isl_int *ineq_i = entry->data;
2654
41.5k
      int neg, more_relaxed;
2655
41.5k
2656
41.5k
      neg = isl_seq_is_neg(ineq_i + 1, ineq + 1, v.len);
2657
41.5k
      if (neg)
2658
190
        isl_int_neg(ineq_i[0], ineq_i[0]);
2659
41.5k
      more_relaxed = isl_int_gt(ineq_i[0], ineq[0]);
2660
41.5k
      if (neg)
2661
190
        isl_int_neg(ineq_i[0], ineq_i[0]);
2662
41.5k
      if (more_relaxed)
2663
2.07k
        break;
2664
41.5k
      else
2665
39.5k
        continue;
2666
3.01k
    }
2667
3.01k
    bound = is_bound(data, set, i, ineq, 0);
2668
3.01k
    if (bound < 0)
2669
0
      return isl_basic_set_free(hull);
2670
3.01k
    
if (3.01k
!bound3.01k
)
2671
2.38k
      break;
2672
44.5k
  }
2673
12.6k
  
if (12.6k
i < set->n12.6k
)
2674
4.46k
    return hull;
2675
8.21k
2676
8.21k
  k = isl_basic_set_alloc_inequality(hull);
2677
8.21k
  if (k < 0)
2678
0
    return isl_basic_set_free(hull);
2679
8.21k
  isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len);
2680
8.21k
2681
8.21k
  return hull;
2682
8.21k
}
2683
2684
/* Compute a superset of the convex hull of "set" that is described
2685
 * by only some of the "n_ineq" constraints in the list "ineq", where "set"
2686
 * has no parameters or integer divisions.
2687
 *
2688
 * The inequalities in "ineq" are assumed to have been sorted such
2689
 * that constraints with the same linear part appear together and
2690
 * that among constraints with the same linear part, those with
2691
 * smaller constant term appear first.
2692
 *
2693
 * We reuse the same data structure that is used by uset_simple_hull,
2694
 * but we do not need the hull table since we will not consider the
2695
 * same constraint more than once.  We therefore allocate it with zero size.
2696
 *
2697
 * We run through the constraints and try to add them one by one,
2698
 * skipping identical constraints.  If we have added a constraint and
2699
 * the next constraint is a more relaxed translate, then we skip this
2700
 * next constraint as well.
2701
 */
2702
static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_constraints(
2703
  __isl_take isl_set *set, int n_ineq, isl_int **ineq)
2704
738
{
2705
738
  int i;
2706
738
  int last_added = 0;
2707
738
  struct sh_data *data = NULL;
2708
738
  isl_basic_set *hull = NULL;
2709
738
  unsigned dim;
2710
738
2711
738
  hull = isl_basic_set_alloc_space(isl_set_get_space(set), 0, 0, n_ineq);
2712
738
  if (!hull)
2713
0
    goto error;
2714
738
2715
738
  data = sh_data_alloc(set, 0);
2716
738
  if (!data)
2717
0
    goto error;
2718
738
2719
738
  dim = isl_set_dim(set, isl_dim_set);
2720
44.4k
  for (i = 0; 
i < n_ineq44.4k
;
++i43.6k
) {
2721
43.6k
    int hull_n_ineq = hull->n_ineq;
2722
43.6k
    int parallel;
2723
43.6k
2724
42.9k
    parallel = i > 0 && isl_seq_eq(ineq[i - 1] + 1, ineq[i] + 1,
2725
42.9k
            dim);
2726
43.6k
    if (parallel &&
2727
33.1k
        
(last_added || 33.1k
isl_int_eq5.33k
(ineq[i - 1][0], ineq[i][0])))
2728
31.0k
      continue;
2729
12.6k
    hull = add_bound_from_constraint(hull, data, set, ineq[i]);
2730
12.6k
    if (!hull)
2731
0
      goto error;
2732
12.6k
    last_added = hull->n_ineq > hull_n_ineq;
2733
12.6k
  }
2734
738
2735
738
  sh_data_free(data);
2736
738
  isl_set_free(set);
2737
738
  return hull;
2738
0
error:
2739
0
  sh_data_free(data);
2740
0
  isl_set_free(set);
2741
0
  isl_basic_set_free(hull);
2742
0
  return NULL;
2743
738
}
2744
2745
/* Collect pointers to all the inequalities in the elements of "list"
2746
 * in "ineq".  For equalities, store both a pointer to the equality and
2747
 * a pointer to its opposite, which is first copied to "mat".
2748
 * "ineq" and "mat" are assumed to have been preallocated to the right size
2749
 * (the number of inequalities + 2 times the number of equalites and
2750
 * the number of equalities, respectively).
2751
 */
2752
static __isl_give isl_mat *collect_inequalities(__isl_take isl_mat *mat,
2753
  __isl_keep isl_basic_set_list *list, isl_int **ineq)
2754
738
{
2755
738
  int i, j, n, n_eq, n_ineq;
2756
738
2757
738
  if (!mat)
2758
0
    return NULL;
2759
738
2760
738
  n_eq = 0;
2761
738
  n_ineq = 0;
2762
738
  n = isl_basic_set_list_n_basic_set(list);
2763
3.99k
  for (i = 0; 
i < n3.99k
;
++i3.25k
) {
2764
3.25k
    isl_basic_set *bset;
2765
3.25k
    bset = isl_basic_set_list_get_basic_set(list, i);
2766
3.25k
    if (!bset)
2767
0
      return isl_mat_free(mat);
2768
4.23k
    
for (j = 0; 3.25k
j < bset->n_eq4.23k
;
++j978
) {
2769
978
      ineq[n_ineq++] = mat->row[n_eq];
2770
978
      ineq[n_ineq++] = bset->eq[j];
2771
978
      isl_seq_neg(mat->row[n_eq++], bset->eq[j], mat->n_col);
2772
978
    }
2773
44.9k
    for (j = 0; 
j < bset->n_ineq44.9k
;
++j41.7k
)
2774
41.7k
      ineq[n_ineq++] = bset->ineq[j];
2775
3.25k
    isl_basic_set_free(bset);
2776
3.25k
  }
2777
738
2778
738
  return mat;
2779
738
}
2780
2781
/* Comparison routine for use as an isl_sort callback.
2782
 *
2783
 * Constraints with the same linear part are sorted together and
2784
 * among constraints with the same linear part, those with smaller
2785
 * constant term are sorted first.
2786
 */
2787
static int cmp_ineq(const void *a, const void *b, void *arg)
2788
251k
{
2789
251k
  unsigned dim = *(unsigned *) arg;
2790
251k
  isl_int * const *ineq1 = a;
2791
251k
  isl_int * const *ineq2 = b;
2792
251k
  int cmp;
2793
251k
2794
251k
  cmp = isl_seq_cmp((*ineq1) + 1, (*ineq2) + 1, dim);
2795
251k
  if (cmp != 0)
2796
196k
    return cmp;
2797
55.1k
  
return 55.1k
isl_int_cmp55.1k
((*ineq1)[0], (*ineq2)[0]);
2798
251k
}
2799
2800
/* Compute a superset of the convex hull of "set" that is described
2801
 * by only constraints in the elements of "list", where "set" has
2802
 * no parameters or integer divisions.
2803
 *
2804
 * We collect all the constraints in those elements and then
2805
 * sort the constraints such that constraints with the same linear part
2806
 * are sorted together and that those with smaller constant term are
2807
 * sorted first.
2808
 */
2809
static __isl_give isl_basic_set *uset_unshifted_simple_hull_from_basic_set_list(
2810
  __isl_take isl_set *set, __isl_take isl_basic_set_list *list)
2811
738
{
2812
738
  int i, n, n_eq, n_ineq;
2813
738
  unsigned dim;
2814
738
  isl_ctx *ctx;
2815
738
  isl_mat *mat = NULL;
2816
738
  isl_int **ineq = NULL;
2817
738
  isl_basic_set *hull;
2818
738
2819
738
  if (!set)
2820
0
    goto error;
2821
738
  ctx = isl_set_get_ctx(set);
2822
738
2823
738
  n_eq = 0;
2824
738
  n_ineq = 0;
2825
738
  n = isl_basic_set_list_n_basic_set(list);
2826
3.99k
  for (i = 0; 
i < n3.99k
;
++i3.25k
) {
2827
3.25k
    isl_basic_set *bset;
2828
3.25k
    bset = isl_basic_set_list_get_basic_set(list, i);
2829
3.25k
    if (!bset)
2830
0
      goto error;
2831
3.25k
    n_eq += bset->n_eq;
2832
3.25k
    n_ineq += 2 * bset->n_eq + bset->n_ineq;
2833
3.25k
    isl_basic_set_free(bset);
2834
3.25k
  }
2835
738
2836
738
  
ineq = 738
isl_alloc_array738
(ctx, isl_int *, n_ineq);
2837
738
  if (
n_ineq > 0 && 738
!ineq738
)
2838
0
    goto error;
2839
738
2840
738
  dim = isl_set_dim(set, isl_dim_set);
2841
738
  mat = isl_mat_alloc(ctx, n_eq, 1 + dim);
2842
738
  mat = collect_inequalities(mat, list, ineq);
2843
738
  if (!mat)
2844
0
    goto error;
2845
738
2846
738
  
if (738
isl_sort(ineq, n_ineq, sizeof(ineq[0]), &cmp_ineq, &dim) < 0738
)
2847
0
    goto error;
2848
738
2849
738
  hull = uset_unshifted_simple_hull_from_constraints(set, n_ineq, ineq);
2850
738
2851
738
  isl_mat_free(mat);
2852
738
  free(ineq);
2853
738
  isl_basic_set_list_free(list);
2854
738
  return hull;
2855
0
error:
2856
0
  isl_mat_free(mat);
2857
0
  free(ineq);
2858
0
  isl_set_free(set);
2859
0
  isl_basic_set_list_free(list);
2860
0
  return NULL;
2861
738
}
2862
2863
/* Compute a superset of the convex hull of "map" that is described
2864
 * by only constraints in the elements of "list".
2865
 *
2866
 * If the list is empty, then we can only describe the universe set.
2867
 * If the input map is empty, then all constraints are valid, so
2868
 * we return the intersection of the elements in "list".
2869
 *
2870
 * Otherwise, we align all divs and temporarily treat them
2871
 * as regular variables, computing the unshifted simple hull in
2872
 * uset_unshifted_simple_hull_from_basic_set_list.
2873
 */
2874
static __isl_give isl_basic_map *map_unshifted_simple_hull_from_basic_map_list(
2875
  __isl_take isl_map *map, __isl_take isl_basic_map_list *list)
2876
738
{
2877
738
  isl_basic_map *model;
2878
738
  isl_basic_map *hull;
2879
738
  isl_set *set;
2880
738
  isl_basic_set_list *bset_list;
2881
738
2882
738
  if (
!map || 738
!list738
)
2883
0
    goto error;
2884
738
2885
738
  
if (738
isl_basic_map_list_n_basic_map(list) == 0738
) {
2886
0
    isl_space *space;
2887
0
2888
0
    space = isl_map_get_space(map);
2889
0
    isl_map_free(map);
2890
0
    isl_basic_map_list_free(list);
2891
0
    return isl_basic_map_universe(space);
2892
0
  }
2893
738
  
if (738
isl_map_plain_is_empty(map)738
) {
2894
0
    isl_map_free(map);
2895
0
    return isl_basic_map_list_intersect(list);
2896
0
  }
2897
738
2898
738
  map = isl_map_align_divs_to_basic_map_list(map, list);
2899
738
  if (!map)
2900
0
    goto error;
2901
738
  list = isl_basic_map_list_align_divs_to_basic_map(list, map->p[0]);
2902
738
2903
738
  model = isl_basic_map_list_get_basic_map(list, 0);
2904
738
2905
738
  set = isl_map_underlying_set(map);
2906
738
  bset_list = isl_basic_map_list_underlying_set(list);
2907
738
2908
738
  hull = uset_unshifted_simple_hull_from_basic_set_list(set, bset_list);
2909
738
  hull = isl_basic_map_overlying_set(hull, model);
2910
738
2911
738
  return hull;
2912
0
error:
2913
0
  isl_map_free(map);
2914
0
  isl_basic_map_list_free(list);
2915
0
  return NULL;
2916
738
}
2917
2918
/* Return a sequence of the basic maps that make up the maps in "list".
2919
 */
2920
static __isl_give isl_basic_map_list *collect_basic_maps(
2921
  __isl_take isl_map_list *list)
2922
738
{
2923
738
  int i, n;
2924
738
  isl_ctx *ctx;
2925
738
  isl_basic_map_list *bmap_list;
2926
738
2927
738
  if (!list)
2928
0
    return NULL;
2929
738
  n = isl_map_list_n_map(list);
2930
738
  ctx = isl_map_list_get_ctx(list);
2931
738
  bmap_list = isl_basic_map_list_alloc(ctx, 0);
2932
738
2933
2.40k
  for (i = 0; 
i < n2.40k
;
++i1.67k
) {
2934
1.67k
    isl_map *map;
2935
1.67k
    isl_basic_map_list *list_i;
2936
1.67k
2937
1.67k
    map = isl_map_list_get_map(list, i);
2938
1.67k
    map = isl_map_compute_divs(map);
2939
1.67k
    list_i = isl_map_get_basic_map_list(map);
2940
1.67k
    isl_map_free(map);
2941
1.67k
    bmap_list = isl_basic_map_list_concat(bmap_list, list_i);
2942
1.67k
  }
2943
738
2944
738
  isl_map_list_free(list);
2945
738
  return bmap_list;
2946
738
}
2947
2948
/* Compute a superset of the convex hull of "map" that is described
2949
 * by only constraints in the elements of "list".
2950
 *
2951
 * If "map" is the universe, then the convex hull (and therefore
2952
 * any superset of the convexhull) is the universe as well.
2953
 *
2954
 * Otherwise, we collect all the basic maps in the map list and
2955
 * continue with map_unshifted_simple_hull_from_basic_map_list.
2956
 */
2957
__isl_give isl_basic_map *isl_map_unshifted_simple_hull_from_map_list(
2958
  __isl_take isl_map *map, __isl_take isl_map_list *list)
2959
749
{
2960
749
  isl_basic_map_list *bmap_list;
2961
749
  int is_universe;
2962
749
2963
749
  is_universe = isl_map_plain_is_universe(map);
2964
749
  if (is_universe < 0)
2965
0
    map = isl_map_free(map);
2966
749
  if (
is_universe < 0 || 749
is_universe749
) {
2967
11
    isl_map_list_free(list);
2968
11
    return isl_map_unshifted_simple_hull(map);
2969
11
  }
2970
738
2971
738
  bmap_list = collect_basic_maps(list);
2972
738
  return map_unshifted_simple_hull_from_basic_map_list(map, bmap_list);
2973
738
}
2974
2975
/* Compute a superset of the convex hull of "set" that is described
2976
 * by only constraints in the elements of "list".
2977
 */
2978
__isl_give isl_basic_set *isl_set_unshifted_simple_hull_from_set_list(
2979
  __isl_take isl_set *set, __isl_take isl_set_list *list)
2980
63
{
2981
63
  return isl_map_unshifted_simple_hull_from_map_list(set, list);
2982
63
}
2983
2984
/* Given a set "set", return parametric bounds on the dimension "dim".
2985
 */
2986
static struct isl_basic_set *set_bounds(struct isl_set *set, int dim)
2987
0
{
2988
0
  unsigned set_dim = isl_set_dim(set, isl_dim_set);
2989
0
  set = isl_set_copy(set);
2990
0
  set = isl_set_eliminate_dims(set, dim + 1, set_dim - (dim + 1));
2991
0
  set = isl_set_eliminate_dims(set, 0, dim);
2992
0
  return isl_set_convex_hull(set);
2993
0
}
2994
2995
/* Computes a "simple hull" and then check if each dimension in the
2996
 * resulting hull is bounded by a symbolic constant.  If not, the
2997
 * hull is intersected with the corresponding bounds on the whole set.
2998
 */
2999
__isl_give isl_basic_set *isl_set_bounded_simple_hull(__isl_take isl_set *set)
3000
0
{
3001
0
  int i, j;
3002
0
  struct isl_basic_set *hull;
3003
0
  unsigned nparam, left;
3004
0
  int removed_divs = 0;
3005
0
3006
0
  hull = isl_set_simple_hull(isl_set_copy(set));
3007
0
  if (!hull)
3008
0
    goto error;
3009
0
3010
0
  nparam = isl_basic_set_dim(hull, isl_dim_param);
3011
0
  for (i = 0; 
i < isl_basic_set_dim(hull, isl_dim_set)0
;
++i0
) {
3012
0
    int lower = 0, upper = 0;
3013
0
    struct isl_basic_set *bounds;
3014
0
3015
0
    left = isl_basic_set_total_dim(hull) - nparam - i - 1;
3016
0
    for (j = 0; 
j < hull->n_eq0
;
++j0
) {
3017
0
      if (isl_int_is_zero(hull->eq[j][1 + nparam + i]))
3018
0
        continue;
3019
0
      
if (0
isl_seq_first_non_zero(hull->eq[j]+1+nparam+i+1,
3020
0
                left) == -1)
3021
0
        break;
3022
0
    }
3023
0
    if (j < hull->n_eq)
3024
0
      continue;
3025
0
3026
0
    
for (j = 0; 0
j < hull->n_ineq0
;
++j0
) {
3027
0
      if (isl_int_is_zero(hull->ineq[j][1 + nparam + i]))
3028
0
        continue;
3029
0
      
if (0
isl_seq_first_non_zero(hull->ineq[j]+1+nparam+i+1,
3030
0
                left) != -1 ||
3031
0
          isl_seq_first_non_zero(hull->ineq[j]+1+nparam,
3032
0
                i) != -1)
3033
0
        continue;
3034
0
      
if (0
isl_int_is_pos0
(hull->ineq[j][1 + nparam + i]))
3035
0
        lower = 1;
3036
0
      else
3037
0
        upper = 1;
3038
0
      if (
lower && 0
upper0
)
3039
0
        break;
3040
0
    }
3041
0
3042
0
    if (
lower && 0
upper0
)
3043
0
      continue;
3044
0
3045
0
    
if (0
!removed_divs0
) {
3046
0
      set = isl_set_remove_divs(set);
3047
0
      if (!set)
3048
0
        goto error;
3049
0
      removed_divs = 1;
3050
0
    }
3051
0
    bounds = set_bounds(set, i);
3052
0
    hull = isl_basic_set_intersect(hull, bounds);
3053
0
    if (!hull)
3054
0
      goto error;
3055
0
  }
3056
0
3057
0
  isl_set_free(set);
3058
0
  return hull;
3059
0
error:
3060
0
  isl_set_free(set);
3061
0
  isl_basic_set_free(hull);
3062
0
  return NULL;
3063
0
}