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

source code of polly/lib/External/isl/isl_sample.c