Coverage Report

Created: 2017-06-23 12:40

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/isl_sample.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 *
4
 * Use of this software is governed by the MIT license
5
 *
6
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8
 */
9
10
#include <isl_ctx_private.h>
11
#include <isl_map_private.h>
12
#include "isl_sample.h"
13
#include <isl/vec.h>
14
#include <isl/mat.h>
15
#include <isl_seq.h>
16
#include "isl_equalities.h"
17
#include "isl_tab.h"
18
#include "isl_basis_reduction.h"
19
#include <isl_factorization.h>
20
#include <isl_point_private.h>
21
#include <isl_options_private.h>
22
#include <isl_vec_private.h>
23
24
#include <bset_from_bmap.c>
25
#include <set_to_map.c>
26
27
static __isl_give isl_vec *empty_sample(__isl_take isl_basic_set *bset)
28
344
{
29
344
  struct isl_vec *vec;
30
344
31
344
  vec = isl_vec_alloc(bset->ctx, 0);
32
344
  isl_basic_set_free(bset);
33
344
  return vec;
34
344
}
35
36
/* Construct a zero sample of the same dimension as bset.
37
 * As a special case, if bset is zero-dimensional, this
38
 * function creates a zero-dimensional sample point.
39
 */
40
static __isl_give isl_vec *zero_sample(__isl_take isl_basic_set *bset)
41
52.6k
{
42
52.6k
  unsigned dim;
43
52.6k
  struct isl_vec *sample;
44
52.6k
45
52.6k
  dim = isl_basic_set_total_dim(bset);
46
52.6k
  sample = isl_vec_alloc(bset->ctx, 1 + dim);
47
52.6k
  if (
sample52.6k
)
{52.6k
48
52.6k
    isl_int_set_si(sample->el[0], 1);
49
52.6k
    isl_seq_clr(sample->el + 1, dim);
50
52.6k
  }
51
52.6k
  isl_basic_set_free(bset);
52
52.6k
  return sample;
53
52.6k
}
54
55
static __isl_give isl_vec *interval_sample(__isl_take isl_basic_set *bset)
56
111k
{
57
111k
  int i;
58
111k
  isl_int t;
59
111k
  struct isl_vec *sample;
60
111k
61
111k
  bset = isl_basic_set_simplify(bset);
62
111k
  if (!bset)
63
0
    return NULL;
64
111k
  
if (111k
isl_basic_set_plain_is_empty(bset)111k
)
65
0
    return empty_sample(bset);
66
111k
  
if (111k
bset->n_eq == 0 && 111k
bset->n_ineq == 0111k
)
67
6.14k
    return zero_sample(bset);
68
111k
69
105k
  sample = isl_vec_alloc(bset->ctx, 2);
70
105k
  if (!sample)
71
0
    goto error;
72
105k
  
if (105k
!bset105k
)
73
0
    return NULL;
74
105k
  
isl_int_set_si105k
(sample->block.data[0], 1);105k
75
105k
76
105k
  if (
bset->n_eq > 0105k
)
{439
77
439
    isl_assert(bset->ctx, bset->n_eq == 1, goto error);
78
439
    
isl_assert439
(bset->ctx, bset->n_ineq == 0, goto error);439
79
439
    
if (439
isl_int_is_one439
(bset->eq[0][1]))
80
439
      isl_int_neg(sample->el[1], bset->eq[0][0]);
81
0
    else {
82
0
      isl_assert(bset->ctx, isl_int_is_negone(bset->eq[0][1]),
83
0
           goto error);
84
0
      
isl_int_set0
(sample->el[1], bset->eq[0][0]);0
85
0
    }
86
439
    isl_basic_set_free(bset);
87
439
    return sample;
88
439
  }
89
105k
90
105k
  
isl_int_init105k
(t);105k
91
105k
  if (isl_int_is_one(bset->ineq[0][1]))
92
91.4k
    isl_int_neg(sample->block.data[1], bset->ineq[0][0]);
93
105k
  else
94
13.6k
    isl_int_set(sample->block.data[1], bset->ineq[0][0]);
95
199k
  for (i = 1; 
i < bset->n_ineq199k
;
++i94.9k
)
{94.9k
96
94.9k
    isl_seq_inner_product(sample->block.data,
97
94.9k
          bset->ineq[i], 2, &t);
98
94.9k
    if (isl_int_is_neg(t))
99
0
      break;
100
94.9k
  }
101
105k
  isl_int_clear(t);
102
105k
  if (
i < bset->n_ineq105k
)
{0
103
0
    isl_vec_free(sample);
104
0
    return empty_sample(bset);
105
0
  }
106
105k
107
105k
  isl_basic_set_free(bset);
108
105k
  return sample;
109
0
error:
110
0
  isl_basic_set_free(bset);
111
0
  isl_vec_free(sample);
112
0
  return NULL;
113
105k
}
114
115
/* Find a sample integer point, if any, in bset, which is known
116
 * to have equalities.  If bset contains no integer points, then
117
 * return a zero-length vector.
118
 * We simply remove the known equalities, compute a sample
119
 * in the resulting bset, using the specified recurse function,
120
 * and then transform the sample back to the original space.
121
 */
122
static __isl_give isl_vec *sample_eq(__isl_take isl_basic_set *bset,
123
  __isl_give isl_vec *(*recurse)(__isl_take isl_basic_set *))
124
60.5k
{
125
60.5k
  struct isl_mat *T;
126
60.5k
  struct isl_vec *sample;
127
60.5k
128
60.5k
  if (!bset)
129
0
    return NULL;
130
60.5k
131
60.5k
  bset = isl_basic_set_remove_equalities(bset, &T, NULL);
132
60.5k
  sample = recurse(bset);
133
60.5k
  if (
!sample || 60.5k
sample->size == 060.5k
)
134
1.38k
    isl_mat_free(T);
135
60.5k
  else
136
59.1k
    sample = isl_mat_vec_product(T, sample);
137
60.5k
  return sample;
138
60.5k
}
139
140
/* Return a matrix containing the equalities of the tableau
141
 * in constraint form.  The tableau is assumed to have
142
 * an associated bset that has been kept up-to-date.
143
 */
144
static struct isl_mat *tab_equalities(struct isl_tab *tab)
145
1.38k
{
146
1.38k
  int i, j;
147
1.38k
  int n_eq;
148
1.38k
  struct isl_mat *eq;
149
1.38k
  struct isl_basic_set *bset;
150
1.38k
151
1.38k
  if (!tab)
152
0
    return NULL;
153
1.38k
154
1.38k
  bset = isl_tab_peek_bset(tab);
155
1.38k
  isl_assert(tab->mat->ctx, bset, return NULL);
156
1.38k
157
1.38k
  n_eq = tab->n_var - tab->n_col + tab->n_dead;
158
1.38k
  if (
tab->empty || 1.38k
n_eq == 01.38k
)
159
231
    return isl_mat_alloc(tab->mat->ctx, 0, tab->n_var);
160
1.15k
  
if (1.15k
n_eq == tab->n_var1.15k
)
161
0
    return isl_mat_identity(tab->mat->ctx, tab->n_var);
162
1.15k
163
1.15k
  eq = isl_mat_alloc(tab->mat->ctx, n_eq, tab->n_var);
164
1.15k
  if (!eq)
165
0
    return NULL;
166
16.7k
  
for (i = 0, j = 0; 1.15k
i < tab->n_con16.7k
;
++i15.5k
)
{15.5k
167
15.5k
    if (tab->con[i].is_row)
168
9.53k
      continue;
169
6.02k
    
if (6.02k
tab->con[i].index >= 0 && 6.02k
tab->con[i].index >= tab->n_dead4.68k
)
170
2.75k
      continue;
171
3.26k
    
if (3.26k
i < bset->n_eq3.26k
)
172
228
      isl_seq_cpy(eq->row[j], bset->eq[i] + 1, tab->n_var);
173
3.26k
    else
174
3.04k
      isl_seq_cpy(eq->row[j],
175
3.04k
            bset->ineq[i - bset->n_eq] + 1, tab->n_var);
176
3.26k
    ++j;
177
3.26k
  }
178
1.15k
  isl_assert(bset->ctx, j == n_eq, goto error);
179
1.15k
  return eq;
180
0
error:
181
0
  isl_mat_free(eq);
182
0
  return NULL;
183
1.15k
}
184
185
/* Compute and return an initial basis for the bounded tableau "tab".
186
 *
187
 * If the tableau is either full-dimensional or zero-dimensional,
188
 * the we simply return an identity matrix.
189
 * Otherwise, we construct a basis whose first directions correspond
190
 * to equalities.
191
 */
192
static struct isl_mat *initial_basis(struct isl_tab *tab)
193
64.9k
{
194
64.9k
  int n_eq;
195
64.9k
  struct isl_mat *eq;
196
64.9k
  struct isl_mat *Q;
197
64.9k
198
64.9k
  tab->n_unbounded = 0;
199
64.9k
  tab->n_zero = n_eq = tab->n_var - tab->n_col + tab->n_dead;
200
64.9k
  if (
tab->empty || 64.9k
n_eq == 064.9k
||
n_eq == tab->n_var871
)
201
64.4k
    return isl_mat_identity(tab->mat->ctx, 1 + tab->n_var);
202
64.9k
203
526
  eq = tab_equalities(tab);
204
526
  eq = isl_mat_left_hermite(eq, 0, NULL, &Q);
205
526
  if (!eq)
206
0
    return NULL;
207
526
  isl_mat_free(eq);
208
526
209
526
  Q = isl_mat_lin_to_aff(Q);
210
526
  return Q;
211
526
}
212
213
/* Compute the minimum of the current ("level") basis row over "tab"
214
 * and store the result in position "level" of "min".
215
 *
216
 * This function assumes that at least one more row and at least
217
 * one more element in the constraint array are available in the tableau.
218
 */
219
static enum isl_lp_result compute_min(isl_ctx *ctx, struct isl_tab *tab,
220
  __isl_keep isl_vec *min, int level)
221
100k
{
222
100k
  return isl_tab_min(tab, tab->basis->row[1 + level],
223
100k
          ctx->one, &min->el[level], NULL, 0);
224
100k
}
225
226
/* Compute the maximum of the current ("level") basis row over "tab"
227
 * and store the result in position "level" of "max".
228
 *
229
 * This function assumes that at least one more row and at least
230
 * one more element in the constraint array are available in the tableau.
231
 */
232
static enum isl_lp_result compute_max(isl_ctx *ctx, struct isl_tab *tab,
233
  __isl_keep isl_vec *max, int level)
234
39.2k
{
235
39.2k
  enum isl_lp_result res;
236
39.2k
  unsigned dim = tab->n_var;
237
39.2k
238
39.2k
  isl_seq_neg(tab->basis->row[1 + level] + 1,
239
39.2k
        tab->basis->row[1 + level] + 1, dim);
240
39.2k
  res = isl_tab_min(tab, tab->basis->row[1 + level],
241
39.2k
        ctx->one, &max->el[level], NULL, 0);
242
39.2k
  isl_seq_neg(tab->basis->row[1 + level] + 1,
243
39.2k
        tab->basis->row[1 + level] + 1, dim);
244
39.2k
  isl_int_neg(max->el[level], max->el[level]);
245
39.2k
246
39.2k
  return res;
247
39.2k
}
248
249
/* Perform a greedy search for an integer point in the set represented
250
 * by "tab", given that the minimal rational value (rounded up to the
251
 * nearest integer) at "level" is smaller than the maximal rational
252
 * value (rounded down to the nearest integer).
253
 *
254
 * Return 1 if we have found an integer point (if tab->n_unbounded > 0
255
 * then we may have only found integer values for the bounded dimensions
256
 * and it is the responsibility of the caller to extend this solution
257
 * to the unbounded dimensions).
258
 * Return 0 if greedy search did not result in a solution.
259
 * Return -1 if some error occurred.
260
 *
261
 * We assign a value half-way between the minimum and the maximum
262
 * to the current dimension and check if the minimal value of the
263
 * next dimension is still smaller than (or equal) to the maximal value.
264
 * We continue this process until either
265
 * - the minimal value (rounded up) is greater than the maximal value
266
 *  (rounded down).  In this case, greedy search has failed.
267
 * - we have exhausted all bounded dimensions, meaning that we have
268
 *  found a solution.
269
 * - the sample value of the tableau is integral.
270
 * - some error has occurred.
271
 */
272
static int greedy_search(isl_ctx *ctx, struct isl_tab *tab,
273
  __isl_keep isl_vec *min, __isl_keep isl_vec *max, int level)
274
8.98k
{
275
8.98k
  struct isl_tab_undo *snap;
276
8.98k
  enum isl_lp_result res;
277
8.98k
278
8.98k
  snap = isl_tab_snap(tab);
279
8.98k
280
25.5k
  do {
281
25.5k
    isl_int_add(tab->basis->row[1 + level][0],
282
25.5k
          min->el[level], max->el[level]);
283
25.5k
    isl_int_fdiv_q_ui(tab->basis->row[1 + level][0],
284
25.5k
          tab->basis->row[1 + level][0], 2);
285
25.5k
    isl_int_neg(tab->basis->row[1 + level][0],
286
25.5k
          tab->basis->row[1 + level][0]);
287
25.5k
    if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
288
0
      return -1;
289
25.5k
    
isl_int_set_si25.5k
(tab->basis->row[1 + level][0], 0);25.5k
290
25.5k
291
25.5k
    if (++level >= tab->n_var - tab->n_unbounded)
292
2.46k
      return 1;
293
23.1k
    
if (23.1k
isl_tab_sample_is_integer(tab)23.1k
)
294
2.56k
      return 1;
295
23.1k
296
20.5k
    res = compute_min(ctx, tab, min, level);
297
20.5k
    if (res == isl_lp_error)
298
0
      return -1;
299
20.5k
    
if (20.5k
res != isl_lp_ok20.5k
)
300
0
      isl_die(ctx, isl_error_internal,
301
20.5k
        "expecting bounded rational solution",
302
20.5k
        return -1);
303
20.5k
    res = compute_max(ctx, tab, max, level);
304
20.5k
    if (res == isl_lp_error)
305
0
      return -1;
306
20.5k
    
if (20.5k
res != isl_lp_ok20.5k
)
307
0
      isl_die(ctx, isl_error_internal,
308
20.5k
        "expecting bounded rational solution",
309
20.5k
        return -1);
310
20.5k
  } while (isl_int_le(min->el[level], max->el[level]));
311
8.98k
312
3.94k
  
if (3.94k
isl_tab_rollback(tab, snap) < 03.94k
)
313
0
    return -1;
314
3.94k
315
3.94k
  return 0;
316
3.94k
}
317
318
/* Given a tableau representing a set, find and return
319
 * an integer point in the set, if there is any.
320
 *
321
 * We perform a depth first search
322
 * for an integer point, by scanning all possible values in the range
323
 * attained by a basis vector, where an initial basis may have been set
324
 * by the calling function.  Otherwise an initial basis that exploits
325
 * the equalities in the tableau is created.
326
 * tab->n_zero is currently ignored and is clobbered by this function.
327
 *
328
 * The tableau is allowed to have unbounded direction, but then
329
 * the calling function needs to set an initial basis, with the
330
 * unbounded directions last and with tab->n_unbounded set
331
 * to the number of unbounded directions.
332
 * Furthermore, the calling functions needs to add shifted copies
333
 * of all constraints involving unbounded directions to ensure
334
 * that any feasible rational value in these directions can be rounded
335
 * up to yield a feasible integer value.
336
 * In particular, let B define the given basis x' = B x
337
 * and let T be the inverse of B, i.e., X = T x'.
338
 * Let a x + c >= 0 be a constraint of the set represented by the tableau,
339
 * or a T x' + c >= 0 in terms of the given basis.  Assume that
340
 * the bounded directions have an integer value, then we can safely
341
 * round up the values for the unbounded directions if we make sure
342
 * that x' not only satisfies the original constraint, but also
343
 * the constraint "a T x' + c + s >= 0" with s the sum of all
344
 * negative values in the last n_unbounded entries of "a T".
345
 * The calling function therefore needs to add the constraint
346
 * a x + c + s >= 0.  The current function then scans the first
347
 * directions for an integer value and once those have been found,
348
 * it can compute "T ceil(B x)" to yield an integer point in the set.
349
 * Note that during the search, the first rows of B may be changed
350
 * by a basis reduction, but the last n_unbounded rows of B remain
351
 * unaltered and are also not mixed into the first rows.
352
 *
353
 * The search is implemented iteratively.  "level" identifies the current
354
 * basis vector.  "init" is true if we want the first value at the current
355
 * level and false if we want the next value.
356
 *
357
 * At the start of each level, we first check if we can find a solution
358
 * using greedy search.  If not, we continue with the exhaustive search.
359
 *
360
 * The initial basis is the identity matrix.  If the range in some direction
361
 * contains more than one integer value, we perform basis reduction based
362
 * on the value of ctx->opt->gbr
363
 *  - ISL_GBR_NEVER:  never perform basis reduction
364
 *  - ISL_GBR_ONCE:   only perform basis reduction the first
365
 *        time such a range is encountered
366
 *  - ISL_GBR_ALWAYS: always perform basis reduction when
367
 *        such a range is encountered
368
 *
369
 * When ctx->opt->gbr is set to ISL_GBR_ALWAYS, then we allow the basis
370
 * reduction computation to return early.  That is, as soon as it
371
 * finds a reasonable first direction.
372
 */ 
373
struct isl_vec *isl_tab_sample(struct isl_tab *tab)
374
77.4k
{
375
77.4k
  unsigned dim;
376
77.4k
  unsigned gbr;
377
77.4k
  struct isl_ctx *ctx;
378
77.4k
  struct isl_vec *sample;
379
77.4k
  struct isl_vec *min;
380
77.4k
  struct isl_vec *max;
381
77.4k
  enum isl_lp_result res;
382
77.4k
  int level;
383
77.4k
  int init;
384
77.4k
  int reduced;
385
77.4k
  struct isl_tab_undo **snap;
386
77.4k
387
77.4k
  if (!tab)
388
0
    return NULL;
389
77.4k
  
if (77.4k
tab->empty77.4k
)
390
7.00k
    return isl_vec_alloc(tab->mat->ctx, 0);
391
77.4k
392
70.4k
  
if (70.4k
!tab->basis70.4k
)
393
64.8k
    tab->basis = initial_basis(tab);
394
70.4k
  if (!tab->basis)
395
0
    return NULL;
396
70.4k
  
isl_assert70.4k
(tab->mat->ctx, tab->basis->n_row == tab->n_var + 1,70.4k
397
70.4k
        return NULL);
398
70.4k
  
isl_assert70.4k
(tab->mat->ctx, tab->basis->n_col == tab->n_var + 1,70.4k
399
70.4k
        return NULL);
400
70.4k
401
70.4k
  ctx = tab->mat->ctx;
402
70.4k
  dim = tab->n_var;
403
70.4k
  gbr = ctx->opt->gbr;
404
70.4k
405
70.4k
  if (
tab->n_unbounded == tab->n_var70.4k
)
{0
406
0
    sample = isl_tab_get_sample_value(tab);
407
0
    sample = isl_mat_vec_product(isl_mat_copy(tab->basis), sample);
408
0
    sample = isl_vec_ceil(sample);
409
0
    sample = isl_mat_vec_inverse_product(isl_mat_copy(tab->basis),
410
0
              sample);
411
0
    return sample;
412
0
  }
413
70.4k
414
70.4k
  
if (70.4k
isl_tab_extend_cons(tab, dim + 1) < 070.4k
)
415
0
    return NULL;
416
70.4k
417
70.4k
  min = isl_vec_alloc(ctx, dim);
418
70.4k
  max = isl_vec_alloc(ctx, dim);
419
70.4k
  snap = isl_alloc_array(ctx, struct isl_tab_undo *, dim);
420
70.4k
421
70.4k
  if (
!min || 70.4k
!max70.4k
||
!snap70.4k
)
422
0
    goto error;
423
70.4k
424
70.4k
  level = 0;
425
70.4k
  init = 1;
426
70.4k
  reduced = 0;
427
70.4k
428
86.2k
  while (
level >= 086.2k
)
{84.5k
429
84.5k
    if (
init84.5k
)
{80.3k
430
80.3k
      int choice;
431
80.3k
432
80.3k
      res = compute_min(ctx, tab, min, level);
433
80.3k
      if (res == isl_lp_error)
434
0
        goto error;
435
80.3k
      
if (80.3k
res != isl_lp_ok80.3k
)
436
0
        isl_die(ctx, isl_error_internal,
437
80.3k
          "expecting bounded rational solution",
438
80.3k
          goto error);
439
80.3k
      
if (80.3k
isl_tab_sample_is_integer(tab)80.3k
)
440
61.6k
        break;
441
18.6k
      res = compute_max(ctx, tab, max, level);
442
18.6k
      if (res == isl_lp_error)
443
0
        goto error;
444
18.6k
      
if (18.6k
res != isl_lp_ok18.6k
)
445
0
        isl_die(ctx, isl_error_internal,
446
18.6k
          "expecting bounded rational solution",
447
18.6k
          goto error);
448
18.6k
      
if (18.6k
isl_tab_sample_is_integer(tab)18.6k
)
449
1.89k
        break;
450
16.7k
      
choice = 16.7k
isl_int_lt16.7k
(min->el[level], max->el[level]);
451
16.7k
      if (
choice16.7k
)
{8.98k
452
8.98k
        int g;
453
8.98k
        g = greedy_search(ctx, tab, min, max, level);
454
8.98k
        if (g < 0)
455
0
          goto error;
456
8.98k
        
if (8.98k
g8.98k
)
457
5.03k
          break;
458
8.98k
      }
459
11.7k
      
if (11.7k
!reduced && 11.7k
choice8.98k
&&
460
3.30k
          
ctx->opt->gbr != 3.30k
ISL_GBR_NEVER3.30k
)
{3.30k
461
3.30k
        unsigned gbr_only_first;
462
3.30k
        if (
ctx->opt->gbr == 3.30k
ISL_GBR_ONCE3.30k
)
463
0
          
ctx->opt->gbr = 0
ISL_GBR_NEVER0
;
464
3.30k
        tab->n_zero = level;
465
3.30k
        gbr_only_first = ctx->opt->gbr_only_first;
466
3.30k
        ctx->opt->gbr_only_first =
467
3.30k
          ctx->opt->gbr == ISL_GBR_ALWAYS;
468
3.30k
        tab = isl_tab_compute_reduced_basis(tab);
469
3.30k
        ctx->opt->gbr_only_first = gbr_only_first;
470
3.30k
        if (
!tab || 3.30k
!tab->basis3.30k
)
471
0
          goto error;
472
3.30k
        reduced = 1;
473
3.30k
        continue;
474
3.30k
      }
475
8.41k
      reduced = 0;
476
8.41k
      snap[level] = isl_tab_snap(tab);
477
8.41k
    } else
478
4.17k
      isl_int_add_ui(min->el[level], min->el[level], 1);
479
84.5k
480
12.5k
    
if (12.5k
isl_int_gt12.5k
(min->el[level], max->el[level]))
{5.90k
481
5.90k
      level--;
482
5.90k
      init = 0;
483
5.90k
      if (level >= 0)
484
4.17k
        
if (4.17k
isl_tab_rollback(tab, snap[level]) < 04.17k
)
485
0
          goto error;
486
5.90k
      continue;
487
5.90k
    }
488
6.68k
    
isl_int_neg6.68k
(tab->basis->row[1 + level][0], min->el[level]);6.68k
489
6.68k
    if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
490
0
      goto error;
491
6.68k
    
isl_int_set_si6.68k
(tab->basis->row[1 + level][0], 0);6.68k
492
6.68k
    if (
level + tab->n_unbounded < dim - 16.68k
)
{6.61k
493
6.61k
      ++level;
494
6.61k
      init = 1;
495
6.61k
      continue;
496
6.61k
    }
497
71
    break;
498
6.68k
  }
499
70.4k
500
70.4k
  
if (70.4k
level >= 070.4k
)
{68.7k
501
68.7k
    sample = isl_tab_get_sample_value(tab);
502
68.7k
    if (!sample)
503
0
      goto error;
504
68.7k
    
if (68.7k
tab->n_unbounded && 68.7k
!596
isl_int_is_one596
(sample->el[0]))
{258
505
258
      sample = isl_mat_vec_product(isl_mat_copy(tab->basis),
506
258
                 sample);
507
258
      sample = isl_vec_ceil(sample);
508
258
      sample = isl_mat_vec_inverse_product(
509
258
          isl_mat_copy(tab->basis), sample);
510
258
    }
511
68.7k
  } else
512
1.72k
    sample = isl_vec_alloc(ctx, 0);
513
70.4k
514
70.4k
  ctx->opt->gbr = gbr;
515
70.4k
  isl_vec_free(min);
516
70.4k
  isl_vec_free(max);
517
70.4k
  free(snap);
518
70.4k
  return sample;
519
0
error:
520
0
  ctx->opt->gbr = gbr;
521
0
  isl_vec_free(min);
522
0
  isl_vec_free(max);
523
0
  free(snap);
524
0
  return NULL;
525
70.4k
}
526
527
static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset);
528
529
/* Compute a sample point of the given basic set, based on the given,
530
 * non-trivial factorization.
531
 */
532
static __isl_give isl_vec *factored_sample(__isl_take isl_basic_set *bset,
533
  __isl_take isl_factorizer *f)
534
25.6k
{
535
25.6k
  int i, n;
536
25.6k
  isl_vec *sample = NULL;
537
25.6k
  isl_ctx *ctx;
538
25.6k
  unsigned nparam;
539
25.6k
  unsigned nvar;
540
25.6k
541
25.6k
  ctx = isl_basic_set_get_ctx(bset);
542
25.6k
  if (!ctx)
543
0
    goto error;
544
25.6k
545
25.6k
  nparam = isl_basic_set_dim(bset, isl_dim_param);
546
25.6k
  nvar = isl_basic_set_dim(bset, isl_dim_set);
547
25.6k
548
25.6k
  sample = isl_vec_alloc(ctx, 1 + isl_basic_set_total_dim(bset));
549
25.6k
  if (!sample)
550
0
    goto error;
551
25.6k
  
isl_int_set_si25.6k
(sample->el[0], 1);25.6k
552
25.6k
553
25.6k
  bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
554
25.6k
555
118k
  for (i = 0, n = 0; 
i < f->n_group118k
;
++i92.8k
)
{93.1k
556
93.1k
    isl_basic_set *bset_i;
557
93.1k
    isl_vec *sample_i;
558
93.1k
559
93.1k
    bset_i = isl_basic_set_copy(bset);
560
93.1k
    bset_i = isl_basic_set_drop_constraints_involving(bset_i,
561
93.1k
          nparam + n + f->len[i], nvar - n - f->len[i]);
562
93.1k
    bset_i = isl_basic_set_drop_constraints_involving(bset_i,
563
93.1k
          nparam, n);
564
93.1k
    bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
565
93.1k
          n + f->len[i], nvar - n - f->len[i]);
566
93.1k
    bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
567
93.1k
568
93.1k
    sample_i = sample_bounded(bset_i);
569
93.1k
    if (!sample_i)
570
0
      goto error;
571
93.1k
    
if (93.1k
sample_i->size == 093.1k
)
{383
572
383
      isl_basic_set_free(bset);
573
383
      isl_factorizer_free(f);
574
383
      isl_vec_free(sample);
575
383
      return sample_i;
576
383
    }
577
92.8k
    isl_seq_cpy(sample->el + 1 + nparam + n,
578
92.8k
          sample_i->el + 1, f->len[i]);
579
92.8k
    isl_vec_free(sample_i);
580
92.8k
581
92.8k
    n += f->len[i];
582
92.8k
  }
583
25.6k
584
25.2k
  f->morph = isl_morph_inverse(f->morph);
585
25.2k
  sample = isl_morph_vec(isl_morph_copy(f->morph), sample);
586
25.2k
587
25.2k
  isl_basic_set_free(bset);
588
25.2k
  isl_factorizer_free(f);
589
25.2k
  return sample;
590
0
error:
591
0
  isl_basic_set_free(bset);
592
0
  isl_factorizer_free(f);
593
0
  isl_vec_free(sample);
594
0
  return NULL;
595
25.6k
}
596
597
/* Given a basic set that is known to be bounded, find and return
598
 * an integer point in the basic set, if there is any.
599
 *
600
 * After handling some trivial cases, we construct a tableau
601
 * and then use isl_tab_sample to find a sample, passing it
602
 * the identity matrix as initial basis.
603
 */ 
604
static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset)
605
189k
{
606
189k
  unsigned dim;
607
189k
  struct isl_vec *sample;
608
189k
  struct isl_tab *tab = NULL;
609
189k
  isl_factorizer *f;
610
189k
611
189k
  if (!bset)
612
0
    return NULL;
613
189k
614
189k
  
if (189k
isl_basic_set_plain_is_empty(bset)189k
)
615
21
    return empty_sample(bset);
616
189k
617
189k
  dim = isl_basic_set_total_dim(bset);
618
189k
  if (dim == 0)
619
42.8k
    return zero_sample(bset);
620
147k
  
if (147k
dim == 1147k
)
621
86.2k
    return interval_sample(bset);
622
60.8k
  
if (60.8k
bset->n_eq > 060.8k
)
623
679
    return sample_eq(bset, sample_bounded);
624
60.8k
625
60.1k
  f = isl_basic_set_factorizer(bset);
626
60.1k
  if (!f)
627
0
    goto error;
628
60.1k
  
if (60.1k
f->n_group != 060.1k
)
629
25.6k
    return factored_sample(bset, f);
630
34.5k
  isl_factorizer_free(f);
631
34.5k
632
34.5k
  tab = isl_tab_from_basic_set(bset, 1);
633
34.5k
  if (
tab && 34.5k
tab->empty34.5k
)
{4.10k
634
4.10k
    isl_tab_free(tab);
635
4.10k
    ISL_F_SET(bset, ISL_BASIC_SET_EMPTY);
636
4.10k
    sample = isl_vec_alloc(isl_basic_set_get_ctx(bset), 0);
637
4.10k
    isl_basic_set_free(bset);
638
4.10k
    return sample;
639
4.10k
  }
640
34.5k
641
30.4k
  
if (30.4k
!30.4k
ISL_F_ISSET30.4k
(bset, ISL_BASIC_SET_NO_IMPLICIT))
642
30.3k
    
if (30.3k
isl_tab_detect_implicit_equalities(tab) < 030.3k
)
643
0
      goto error;
644
30.4k
645
30.4k
  sample = isl_tab_sample(tab);
646
30.4k
  if (!sample)
647
0
    goto error;
648
30.4k
649
30.4k
  
if (30.4k
sample->size > 030.4k
)
{29.3k
650
29.3k
    isl_vec_free(bset->sample);
651
29.3k
    bset->sample = isl_vec_copy(sample);
652
29.3k
  }
653
30.4k
654
30.4k
  isl_basic_set_free(bset);
655
30.4k
  isl_tab_free(tab);
656
30.4k
  return sample;
657
0
error:
658
0
  isl_basic_set_free(bset);
659
0
  isl_tab_free(tab);
660
0
  return NULL;
661
30.4k
}
662
663
/* Given a basic set "bset" and a value "sample" for the first coordinates
664
 * of bset, plug in these values and drop the corresponding coordinates.
665
 *
666
 * We do this by computing the preimage of the transformation
667
 *
668
 *       [ 1 0 ]
669
 *  x =  [ s 0 ] x'
670
 *       [ 0 I ]
671
 *
672
 * where [1 s] is the sample value and I is the identity matrix of the
673
 * appropriate dimension.
674
 */
675
static __isl_give isl_basic_set *plug_in(__isl_take isl_basic_set *bset,
676
  __isl_take isl_vec *sample)
677
77.8k
{
678
77.8k
  int i;
679
77.8k
  unsigned total;
680
77.8k
  struct isl_mat *T;
681
77.8k
682
77.8k
  if (
!bset || 77.8k
!sample77.8k
)
683
0
    goto error;
684
77.8k
685
77.8k
  total = isl_basic_set_total_dim(bset);
686
77.8k
  T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1));
687
77.8k
  if (!T)
688
0
    goto error;
689
77.8k
690
271k
  
for (i = 0; 77.8k
i < sample->size271k
;
++i193k
)
{193k
691
193k
    isl_int_set(T->row[i][0], sample->el[i]);
692
193k
    isl_seq_clr(T->row[i] + 1, T->n_col - 1);
693
193k
  }
694
366k
  for (i = 0; 
i < T->n_col - 1366k
;
++i288k
)
{288k
695
288k
    isl_seq_clr(T->row[sample->size + i], T->n_col);
696
288k
    isl_int_set_si(T->row[sample->size + i][1 + i], 1);
697
288k
  }
698
77.8k
  isl_vec_free(sample);
699
77.8k
700
77.8k
  bset = isl_basic_set_preimage(bset, T);
701
77.8k
  return bset;
702
0
error:
703
0
  isl_basic_set_free(bset);
704
0
  isl_vec_free(sample);
705
0
  return NULL;
706
77.8k
}
707
708
/* Given a basic set "bset", return any (possibly non-integer) point
709
 * in the basic set.
710
 */
711
static __isl_give isl_vec *rational_sample(__isl_take isl_basic_set *bset)
712
81.3k
{
713
81.3k
  struct isl_tab *tab;
714
81.3k
  struct isl_vec *sample;
715
81.3k
716
81.3k
  if (!bset)
717
0
    return NULL;
718
81.3k
719
81.3k
  tab = isl_tab_from_basic_set(bset, 0);
720
81.3k
  sample = isl_tab_get_sample_value(tab);
721
81.3k
  isl_tab_free(tab);
722
81.3k
723
81.3k
  isl_basic_set_free(bset);
724
81.3k
725
81.3k
  return sample;
726
81.3k
}
727
728
/* Given a linear cone "cone" and a rational point "vec",
729
 * construct a polyhedron with shifted copies of the constraints in "cone",
730
 * i.e., a polyhedron with "cone" as its recession cone, such that each
731
 * point x in this polyhedron is such that the unit box positioned at x
732
 * lies entirely inside the affine cone 'vec + cone'.
733
 * Any rational point in this polyhedron may therefore be rounded up
734
 * to yield an integer point that lies inside said affine cone.
735
 *
736
 * Denote the constraints of cone by "<a_i, x> >= 0" and the rational
737
 * point "vec" by v/d.
738
 * Let b_i = <a_i, v>.  Then the affine cone 'vec + cone' is given
739
 * by <a_i, x> - b/d >= 0.
740
 * The polyhedron <a_i, x> - ceil{b/d} >= 0 is a subset of this affine cone.
741
 * We prefer this polyhedron over the actual affine cone because it doesn't
742
 * require a scaling of the constraints.
743
 * If each of the vertices of the unit cube positioned at x lies inside
744
 * this polyhedron, then the whole unit cube at x lies inside the affine cone.
745
 * We therefore impose that x' = x + \sum e_i, for any selection of unit
746
 * vectors lies inside the polyhedron, i.e.,
747
 *
748
 *  <a_i, x'> - ceil{b/d} = <a_i, x> + sum a_i - ceil{b/d} >= 0
749
 *
750
 * The most stringent of these constraints is the one that selects
751
 * all negative a_i, so the polyhedron we are looking for has constraints
752
 *
753
 *  <a_i, x> + sum_{a_i < 0} a_i - ceil{b/d} >= 0
754
 *
755
 * Note that if cone were known to have only non-negative rays
756
 * (which can be accomplished by a unimodular transformation),
757
 * then we would only have to check the points x' = x + e_i
758
 * and we only have to add the smallest negative a_i (if any)
759
 * instead of the sum of all negative a_i.
760
 */
761
static __isl_give isl_basic_set *shift_cone(__isl_take isl_basic_set *cone,
762
  __isl_take isl_vec *vec)
763
3.47k
{
764
3.47k
  int i, j, k;
765
3.47k
  unsigned total;
766
3.47k
767
3.47k
  struct isl_basic_set *shift = NULL;
768
3.47k
769
3.47k
  if (
!cone || 3.47k
!vec3.47k
)
770
0
    goto error;
771
3.47k
772
3.47k
  
isl_assert3.47k
(cone->ctx, cone->n_eq == 0, goto error);3.47k
773
3.47k
774
3.47k
  total = isl_basic_set_total_dim(cone);
775
3.47k
776
3.47k
  shift = isl_basic_set_alloc_space(isl_basic_set_get_space(cone),
777
3.47k
          0, 0, cone->n_ineq);
778
3.47k
779
12.1k
  for (i = 0; 
i < cone->n_ineq12.1k
;
++i8.66k
)
{8.66k
780
8.66k
    k = isl_basic_set_alloc_inequality(shift);
781
8.66k
    if (k < 0)
782
0
      goto error;
783
8.66k
    isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total);
784
8.66k
    isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total,
785
8.66k
              &shift->ineq[k][0]);
786
8.66k
    isl_int_cdiv_q(shift->ineq[k][0],
787
8.66k
             shift->ineq[k][0], vec->el[0]);
788
8.66k
    isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]);
789
37.3k
    for (j = 0; 
j < total37.3k
;
++j28.6k
)
{28.6k
790
28.6k
      if (isl_int_is_nonneg(shift->ineq[k][1 + j]))
791
21.2k
        continue;
792
7.42k
      
isl_int_add7.42k
(shift->ineq[k][0],7.42k
793
7.42k
            shift->ineq[k][0], shift->ineq[k][1 + j]);
794
7.42k
    }
795
8.66k
  }
796
3.47k
797
3.47k
  isl_basic_set_free(cone);
798
3.47k
  isl_vec_free(vec);
799
3.47k
800
3.47k
  return isl_basic_set_finalize(shift);
801
0
error:
802
0
  isl_basic_set_free(shift);
803
0
  isl_basic_set_free(cone);
804
0
  isl_vec_free(vec);
805
0
  return NULL;
806
3.47k
}
807
808
/* Given a rational point vec in a (transformed) basic set,
809
 * such that cone is the recession cone of the original basic set,
810
 * "round up" the rational point to an integer point.
811
 *
812
 * We first check if the rational point just happens to be integer.
813
 * If not, we transform the cone in the same way as the basic set,
814
 * pick a point x in this cone shifted to the rational point such that
815
 * the whole unit cube at x is also inside this affine cone.
816
 * Then we simply round up the coordinates of x and return the
817
 * resulting integer point.
818
 */
819
static __isl_give isl_vec *round_up_in_cone(__isl_take isl_vec *vec,
820
  __isl_take isl_basic_set *cone, __isl_take isl_mat *U)
821
77.8k
{
822
77.8k
  unsigned total;
823
77.8k
824
77.8k
  if (
!vec || 77.8k
!cone77.8k
||
!U77.8k
)
825
0
    goto error;
826
77.8k
827
77.8k
  
isl_assert77.8k
(vec->ctx, vec->size != 0, goto error);77.8k
828
77.8k
  
if (77.8k
isl_int_is_one77.8k
(vec->el[0]))
{74.3k
829
74.3k
    isl_mat_free(U);
830
74.3k
    isl_basic_set_free(cone);
831
74.3k
    return vec;
832
74.3k
  }
833
77.8k
834
3.47k
  total = isl_basic_set_total_dim(cone);
835
3.47k
  cone = isl_basic_set_preimage(cone, U);
836
3.47k
  cone = isl_basic_set_remove_dims(cone, isl_dim_set,
837
3.47k
           0, total - (vec->size - 1));
838
3.47k
839
3.47k
  cone = shift_cone(cone, vec);
840
3.47k
841
3.47k
  vec = rational_sample(cone);
842
3.47k
  vec = isl_vec_ceil(vec);
843
3.47k
  return vec;
844
0
error:
845
0
  isl_mat_free(U);
846
0
  isl_vec_free(vec);
847
0
  isl_basic_set_free(cone);
848
0
  return NULL;
849
77.8k
}
850
851
/* Concatenate two integer vectors, i.e., two vectors with denominator
852
 * (stored in element 0) equal to 1.
853
 */
854
static __isl_give isl_vec *vec_concat(__isl_take isl_vec *vec1,
855
  __isl_take isl_vec *vec2)
856
77.8k
{
857
77.8k
  struct isl_vec *vec;
858
77.8k
859
77.8k
  if (
!vec1 || 77.8k
!vec277.8k
)
860
0
    goto error;
861
77.8k
  
isl_assert77.8k
(vec1->ctx, vec1->size > 0, goto error);77.8k
862
77.8k
  
isl_assert77.8k
(vec2->ctx, vec2->size > 0, goto error);77.8k
863
77.8k
  
isl_assert77.8k
(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error);77.8k
864
77.8k
  
isl_assert77.8k
(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error);77.8k
865
77.8k
866
77.8k
  vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1);
867
77.8k
  if (!vec)
868
0
    goto error;
869
77.8k
870
77.8k
  isl_seq_cpy(vec->el, vec1->el, vec1->size);
871
77.8k
  isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1);
872
77.8k
873
77.8k
  isl_vec_free(vec1);
874
77.8k
  isl_vec_free(vec2);
875
77.8k
876
77.8k
  return vec;
877
0
error:
878
0
  isl_vec_free(vec1);
879
0
  isl_vec_free(vec2);
880
0
  return NULL;
881
77.8k
}
882
883
/* Give a basic set "bset" with recession cone "cone", compute and
884
 * return an integer point in bset, if any.
885
 *
886
 * If the recession cone is full-dimensional, then we know that
887
 * bset contains an infinite number of integer points and it is
888
 * fairly easy to pick one of them.
889
 * If the recession cone is not full-dimensional, then we first
890
 * transform bset such that the bounded directions appear as
891
 * the first dimensions of the transformed basic set.
892
 * We do this by using a unimodular transformation that transforms
893
 * the equalities in the recession cone to equalities on the first
894
 * dimensions.
895
 *
896
 * The transformed set is then projected onto its bounded dimensions.
897
 * Note that to compute this projection, we can simply drop all constraints
898
 * involving any of the unbounded dimensions since these constraints
899
 * cannot be combined to produce a constraint on the bounded dimensions.
900
 * To see this, assume that there is such a combination of constraints
901
 * that produces a constraint on the bounded dimensions.  This means
902
 * that some combination of the unbounded dimensions has both an upper
903
 * bound and a lower bound in terms of the bounded dimensions, but then
904
 * this combination would be a bounded direction too and would have been
905
 * transformed into a bounded dimensions.
906
 *
907
 * We then compute a sample value in the bounded dimensions.
908
 * If no such value can be found, then the original set did not contain
909
 * any integer points and we are done.
910
 * Otherwise, we plug in the value we found in the bounded dimensions,
911
 * project out these bounded dimensions and end up with a set with
912
 * a full-dimensional recession cone.
913
 * A sample point in this set is computed by "rounding up" any
914
 * rational point in the set.
915
 *
916
 * The sample points in the bounded and unbounded dimensions are
917
 * then combined into a single sample point and transformed back
918
 * to the original space.
919
 */
920
__isl_give isl_vec *isl_basic_set_sample_with_cone(
921
  __isl_take isl_basic_set *bset, __isl_take isl_basic_set *cone)
922
82.0k
{
923
82.0k
  unsigned total;
924
82.0k
  unsigned cone_dim;
925
82.0k
  struct isl_mat *M, *U;
926
82.0k
  struct isl_vec *sample;
927
82.0k
  struct isl_vec *cone_sample;
928
82.0k
  struct isl_ctx *ctx;
929
82.0k
  struct isl_basic_set *bounded;
930
82.0k
931
82.0k
  if (
!bset || 82.0k
!cone82.0k
)
932
0
    goto error;
933
82.0k
934
82.0k
  ctx = isl_basic_set_get_ctx(bset);
935
82.0k
  total = isl_basic_set_total_dim(cone);
936
82.0k
  cone_dim = total - cone->n_eq;
937
82.0k
938
82.0k
  M = isl_mat_sub_alloc6(ctx, cone->eq, 0, cone->n_eq, 1, total);
939
82.0k
  M = isl_mat_left_hermite(M, 0, &U, NULL);
940
82.0k
  if (!M)
941
0
    goto error;
942
82.0k
  isl_mat_free(M);
943
82.0k
944
82.0k
  U = isl_mat_lin_to_aff(U);
945
82.0k
  bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
946
82.0k
947
82.0k
  bounded = isl_basic_set_copy(bset);
948
82.0k
  bounded = isl_basic_set_drop_constraints_involving(bounded,
949
82.0k
               total - cone_dim, cone_dim);
950
82.0k
  bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim);
951
82.0k
  sample = sample_bounded(bounded);
952
82.0k
  if (
!sample || 82.0k
sample->size == 082.0k
)
{4.18k
953
4.18k
    isl_basic_set_free(bset);
954
4.18k
    isl_basic_set_free(cone);
955
4.18k
    isl_mat_free(U);
956
4.18k
    return sample;
957
4.18k
  }
958
77.8k
  bset = plug_in(bset, isl_vec_copy(sample));
959
77.8k
  cone_sample = rational_sample(bset);
960
77.8k
  cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(U));
961
77.8k
  sample = vec_concat(sample, cone_sample);
962
77.8k
  sample = isl_mat_vec_product(U, sample);
963
77.8k
  return sample;
964
0
error:
965
0
  isl_basic_set_free(cone);
966
0
  isl_basic_set_free(bset);
967
0
  return NULL;
968
82.0k
}
969
970
static void vec_sum_of_neg(struct isl_vec *v, isl_int *s)
971
1.79k
{
972
1.79k
  int i;
973
1.79k
974
1.79k
  isl_int_set_si(*s, 0);
975
1.79k
976
5.23k
  for (i = 0; 
i < v->size5.23k
;
++i3.44k
)
977
3.44k
    
if (3.44k
isl_int_is_neg3.44k
(v->el[i]))
978
531
      isl_int_add(*s, *s, v->el[i]);
979
1.79k
}
980
981
/* Given a tableau "tab", a tableau "tab_cone" that corresponds
982
 * to the recession cone and the inverse of a new basis U = inv(B),
983
 * with the unbounded directions in B last,
984
 * add constraints to "tab" that ensure any rational value
985
 * in the unbounded directions can be rounded up to an integer value.
986
 *
987
 * The new basis is given by x' = B x, i.e., x = U x'.
988
 * For any rational value of the last tab->n_unbounded coordinates
989
 * in the update tableau, the value that is obtained by rounding
990
 * up this value should be contained in the original tableau.
991
 * For any constraint "a x + c >= 0", we therefore need to add
992
 * a constraint "a x + c + s >= 0", with s the sum of all negative
993
 * entries in the last elements of "a U".
994
 *
995
 * Since we are not interested in the first entries of any of the "a U",
996
 * we first drop the columns of U that correpond to bounded directions.
997
 */
998
static int tab_shift_cone(struct isl_tab *tab,
999
  struct isl_tab *tab_cone, struct isl_mat *U)
1000
428
{
1001
428
  int i;
1002
428
  isl_int v;
1003
428
  struct isl_basic_set *bset = NULL;
1004
428
1005
428
  if (
tab && 428
tab->n_unbounded == 0428
)
{0
1006
0
    isl_mat_free(U);
1007
0
    return 0;
1008
0
  }
1009
428
  
isl_int_init428
(v);428
1010
428
  if (
!tab || 428
!tab_cone428
||
!U428
)
1011
0
    goto error;
1012
428
  bset = isl_tab_peek_bset(tab_cone);
1013
428
  U = isl_mat_drop_cols(U, 0, tab->n_var - tab->n_unbounded);
1014
6.04k
  for (i = 0; 
i < bset->n_ineq6.04k
;
++i5.61k
)
{5.61k
1015
5.61k
    int ok;
1016
5.61k
    struct isl_vec *row = NULL;
1017
5.61k
    if (isl_tab_is_equality(tab_cone, tab_cone->n_eq + i))
1018
3.82k
      continue;
1019
1.79k
    row = isl_vec_alloc(bset->ctx, tab_cone->n_var);
1020
1.79k
    if (!row)
1021
0
      goto error;
1022
1.79k
    isl_seq_cpy(row->el, bset->ineq[i] + 1, tab_cone->n_var);
1023
1.79k
    row = isl_vec_mat_product(row, isl_mat_copy(U));
1024
1.79k
    if (!row)
1025
0
      goto error;
1026
1.79k
    vec_sum_of_neg(row, &v);
1027
1.79k
    isl_vec_free(row);
1028
1.79k
    if (isl_int_is_zero(v))
1029
1.32k
      continue;
1030
465
    
if (465
isl_tab_extend_cons(tab, 1) < 0465
)
1031
0
      goto error;
1032
465
    
isl_int_add465
(bset->ineq[i][0], bset->ineq[i][0], v);465
1033
465
    ok = isl_tab_add_ineq(tab, bset->ineq[i]) >= 0;
1034
465
    isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], v);
1035
465
    if (!ok)
1036
0
      goto error;
1037
465
  }
1038
428
1039
428
  isl_mat_free(U);
1040
428
  isl_int_clear(v);
1041
428
  return 0;
1042
0
error:
1043
0
  isl_mat_free(U);
1044
0
  isl_int_clear(v);
1045
0
  return -1;
1046
428
}
1047
1048
/* Compute and return an initial basis for the possibly
1049
 * unbounded tableau "tab".  "tab_cone" is a tableau
1050
 * for the corresponding recession cone.
1051
 * Additionally, add constraints to "tab" that ensure
1052
 * that any rational value for the unbounded directions
1053
 * can be rounded up to an integer value.
1054
 *
1055
 * If the tableau is bounded, i.e., if the recession cone
1056
 * is zero-dimensional, then we just use inital_basis.
1057
 * Otherwise, we construct a basis whose first directions
1058
 * correspond to equalities, followed by bounded directions,
1059
 * i.e., equalities in the recession cone.
1060
 * The remaining directions are then unbounded.
1061
 */
1062
int isl_tab_set_initial_basis_with_cone(struct isl_tab *tab,
1063
  struct isl_tab *tab_cone)
1064
543
{
1065
543
  struct isl_mat *eq;
1066
543
  struct isl_mat *cone_eq;
1067
543
  struct isl_mat *U, *Q;
1068
543
1069
543
  if (
!tab || 543
!tab_cone543
)
1070
0
    return -1;
1071
543
1072
543
  
if (543
tab_cone->n_col == tab_cone->n_dead543
)
{115
1073
115
    tab->basis = initial_basis(tab);
1074
115
    return tab->basis ? 
0115
:
-10
;
1075
115
  }
1076
543
1077
428
  eq = tab_equalities(tab);
1078
428
  if (!eq)
1079
0
    return -1;
1080
428
  tab->n_zero = eq->n_row;
1081
428
  cone_eq = tab_equalities(tab_cone);
1082
428
  eq = isl_mat_concat(eq, cone_eq);
1083
428
  if (!eq)
1084
0
    return -1;
1085
428
  tab->n_unbounded = tab->n_var - (eq->n_row - tab->n_zero);
1086
428
  eq = isl_mat_left_hermite(eq, 0, &U, &Q);
1087
428
  if (!eq)
1088
0
    return -1;
1089
428
  isl_mat_free(eq);
1090
428
  tab->basis = isl_mat_lin_to_aff(Q);
1091
428
  if (tab_shift_cone(tab, tab_cone, U) < 0)
1092
0
    return -1;
1093
428
  
if (428
!tab->basis428
)
1094
0
    return -1;
1095
428
  return 0;
1096
428
}
1097
1098
/* Compute and return a sample point in bset using generalized basis
1099
 * reduction.  We first check if the input set has a non-trivial
1100
 * recession cone.  If so, we perform some extra preprocessing in
1101
 * sample_with_cone.  Otherwise, we directly perform generalized basis
1102
 * reduction.
1103
 */
1104
static __isl_give isl_vec *gbr_sample(__isl_take isl_basic_set *bset)
1105
95.3k
{
1106
95.3k
  unsigned dim;
1107
95.3k
  struct isl_basic_set *cone;
1108
95.3k
1109
95.3k
  dim = isl_basic_set_total_dim(bset);
1110
95.3k
1111
95.3k
  cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
1112
95.3k
  if (!cone)
1113
0
    goto error;
1114
95.3k
1115
95.3k
  
if (95.3k
cone->n_eq < dim95.3k
)
1116
81.2k
    return isl_basic_set_sample_with_cone(bset, cone);
1117
95.3k
1118
14.0k
  isl_basic_set_free(cone);
1119
14.0k
  return sample_bounded(bset);
1120
0
error:
1121
0
  isl_basic_set_free(bset);
1122
0
  return NULL;
1123
95.3k
}
1124
1125
static __isl_give isl_vec *basic_set_sample(__isl_take isl_basic_set *bset,
1126
  int bounded)
1127
184k
{
1128
184k
  struct isl_ctx *ctx;
1129
184k
  unsigned dim;
1130
184k
  if (!bset)
1131
0
    return NULL;
1132
184k
1133
184k
  ctx = bset->ctx;
1134
184k
  if (isl_basic_set_plain_is_empty(bset))
1135
323
    return empty_sample(bset);
1136
184k
1137
184k
  dim = isl_basic_set_n_dim(bset);
1138
184k
  isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
1139
184k
  
isl_assert184k
(ctx, bset->n_div == 0, goto error);184k
1140
184k
1141
184k
  
if (184k
bset->sample && 184k
bset->sample->size == 1 + dim171
)
{35
1142
35
    int contains = isl_basic_set_contains(bset, bset->sample);
1143
35
    if (contains < 0)
1144
0
      goto error;
1145
35
    
if (35
contains35
)
{0
1146
0
      struct isl_vec *sample = isl_vec_copy(bset->sample);
1147
0
      isl_basic_set_free(bset);
1148
0
      return sample;
1149
0
    }
1150
35
  }
1151
184k
  isl_vec_free(bset->sample);
1152
184k
  bset->sample = NULL;
1153
184k
1154
184k
  if (bset->n_eq > 0)
1155
59.8k
    
return sample_eq(bset, bounded ? 59.8k
isl_basic_set_sample_bounded0
1156
59.8k
                 : isl_basic_set_sample_vec);
1157
124k
  
if (124k
dim == 0124k
)
1158
3.70k
    return zero_sample(bset);
1159
120k
  
if (120k
dim == 1120k
)
1160
25.3k
    return interval_sample(bset);
1161
120k
1162
95.3k
  
return bounded ? 95.3k
sample_bounded(bset)0
:
gbr_sample(bset)95.3k
;
1163
0
error:
1164
0
  isl_basic_set_free(bset);
1165
0
  return NULL;
1166
120k
}
1167
1168
__isl_give isl_vec *isl_basic_set_sample_vec(__isl_take isl_basic_set *bset)
1169
184k
{
1170
184k
  return basic_set_sample(bset, 0);
1171
184k
}
1172
1173
/* Compute an integer sample in "bset", where the caller guarantees
1174
 * that "bset" is bounded.
1175
 */
1176
__isl_give isl_vec *isl_basic_set_sample_bounded(__isl_take isl_basic_set *bset)
1177
0
{
1178
0
  return basic_set_sample(bset, 1);
1179
0
}
1180
1181
__isl_give isl_basic_set *isl_basic_set_from_vec(__isl_take isl_vec *vec)
1182
154k
{
1183
154k
  int i;
1184
154k
  int k;
1185
154k
  struct isl_basic_set *bset = NULL;
1186
154k
  struct isl_ctx *ctx;
1187
154k
  unsigned dim;
1188
154k
1189
154k
  if (!vec)
1190
0
    return NULL;
1191
154k
  ctx = vec->ctx;
1192
154k
  isl_assert(ctx, vec->size != 0, goto error);
1193
154k
1194
154k
  bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
1195
154k
  if (!bset)
1196
0
    goto error;
1197
154k
  dim = isl_basic_set_n_dim(bset);
1198
597k
  for (i = dim - 1; 
i >= 0597k
;
--i442k
)
{442k
1199
442k
    k = isl_basic_set_alloc_equality(bset);
1200
442k
    if (k < 0)
1201
0
      goto error;
1202
442k
    isl_seq_clr(bset->eq[k], 1 + dim);
1203
442k
    isl_int_neg(bset->eq[k][0], vec->el[1 + i]);
1204
442k
    isl_int_set(bset->eq[k][1 + i], vec->el[0]);
1205
442k
  }
1206
154k
  bset->sample = vec;
1207
154k
1208
154k
  return bset;
1209
0
error:
1210
0
  isl_basic_set_free(bset);
1211
0
  isl_vec_free(vec);
1212
0
  return NULL;
1213
154k
}
1214
1215
__isl_give isl_basic_map *isl_basic_map_sample(__isl_take isl_basic_map *bmap)
1216
1
{
1217
1
  struct isl_basic_set *bset;
1218
1
  struct isl_vec *sample_vec;
1219
1
1220
1
  bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
1221
1
  sample_vec = isl_basic_set_sample_vec(bset);
1222
1
  if (!sample_vec)
1223
0
    goto error;
1224
1
  
if (1
sample_vec->size == 01
)
{0
1225
0
    isl_vec_free(sample_vec);
1226
0
    return isl_basic_map_set_to_empty(bmap);
1227
0
  }
1228
1
  isl_vec_free(bmap->sample);
1229
1
  bmap->sample = isl_vec_copy(sample_vec);
1230
1
  bset = isl_basic_set_from_vec(sample_vec);
1231
1
  return isl_basic_map_overlying_set(bset, bmap);
1232
0
error:
1233
0
  isl_basic_map_free(bmap);
1234
0
  return NULL;
1235
1
}
1236
1237
__isl_give isl_basic_set *isl_basic_set_sample(__isl_take isl_basic_set *bset)
1238
1
{
1239
1
  return isl_basic_map_sample(bset);
1240
1
}
1241
1242
__isl_give isl_basic_map *isl_map_sample(__isl_take isl_map *map)
1243
0
{
1244
0
  int i;
1245
0
  isl_basic_map *sample = NULL;
1246
0
1247
0
  if (!map)
1248
0
    goto error;
1249
0
1250
0
  
for (i = 0; 0
i < map->n0
;
++i0
)
{0
1251
0
    sample = isl_basic_map_sample(isl_basic_map_copy(map->p[i]));
1252
0
    if (!sample)
1253
0
      goto error;
1254
0
    
if (0
!0
ISL_F_ISSET0
(sample, ISL_BASIC_MAP_EMPTY))
1255
0
      break;
1256
0
    isl_basic_map_free(sample);
1257
0
  }
1258
0
  
if (0
i == map->n0
)
1259
0
    sample = isl_basic_map_empty(isl_map_get_space(map));
1260
0
  isl_map_free(map);
1261
0
  return sample;
1262
0
error:
1263
0
  isl_map_free(map);
1264
0
  return NULL;
1265
0
}
1266
1267
__isl_give isl_basic_set *isl_set_sample(__isl_take isl_set *set)
1268
0
{
1269
0
  return bset_from_bmap(isl_map_sample(set_to_map(set)));
1270
0
}
1271
1272
__isl_give isl_point *isl_basic_set_sample_point(__isl_take isl_basic_set *bset)
1273
1
{
1274
1
  isl_vec *vec;
1275
1
  isl_space *dim;
1276
1
1277
1
  dim = isl_basic_set_get_space(bset);
1278
1
  bset = isl_basic_set_underlying_set(bset);
1279
1
  vec = isl_basic_set_sample_vec(bset);
1280
1
1281
1
  return isl_point_alloc(dim, vec);
1282
1
}
1283
1284
__isl_give isl_point *isl_set_sample_point(__isl_take isl_set *set)
1285
2
{
1286
2
  int i;
1287
2
  isl_point *pnt;
1288
2
1289
2
  if (!set)
1290
0
    return NULL;
1291
2
1292
2
  
for (i = 0; 2
i < set->n2
;
++i0
)
{1
1293
1
    pnt = isl_basic_set_sample_point(isl_basic_set_copy(set->p[i]));
1294
1
    if (!pnt)
1295
0
      goto error;
1296
1
    
if (1
!isl_point_is_void(pnt)1
)
1297
1
      break;
1298
0
    isl_point_free(pnt);
1299
0
  }
1300
2
  
if (2
i == set->n2
)
1301
1
    pnt = isl_point_void(isl_set_get_space(set));
1302
2
1303
2
  isl_set_free(set);
1304
2
  return pnt;
1305
0
error:
1306
0
  isl_set_free(set);
1307
0
  return NULL;
1308
2
}