Line data Source code
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 166168 : static __isl_give isl_vec *empty_sample(__isl_take isl_basic_set *bset)
28 : {
29 : struct isl_vec *vec;
30 :
31 166168 : vec = isl_vec_alloc(bset->ctx, 0);
32 166168 : isl_basic_set_free(bset);
33 166168 : return vec;
34 : }
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 10587440 : static __isl_give isl_vec *zero_sample(__isl_take isl_basic_set *bset)
41 : {
42 : unsigned dim;
43 : struct isl_vec *sample;
44 :
45 10587440 : dim = isl_basic_set_total_dim(bset);
46 10587440 : sample = isl_vec_alloc(bset->ctx, 1 + dim);
47 10587440 : if (sample) {
48 10587440 : isl_int_set_si(sample->el[0], 1);
49 10587440 : isl_seq_clr(sample->el + 1, dim);
50 : }
51 10587440 : isl_basic_set_free(bset);
52 10587440 : return sample;
53 : }
54 :
55 1058979 : static __isl_give isl_vec *interval_sample(__isl_take isl_basic_set *bset)
56 : {
57 : int i;
58 : isl_int t;
59 : struct isl_vec *sample;
60 :
61 1058979 : bset = isl_basic_set_simplify(bset);
62 1058979 : if (!bset)
63 0 : return NULL;
64 1058979 : if (isl_basic_set_plain_is_empty(bset))
65 0 : return empty_sample(bset);
66 1058979 : if (bset->n_eq == 0 && bset->n_ineq == 0)
67 1083 : return zero_sample(bset);
68 :
69 1057896 : sample = isl_vec_alloc(bset->ctx, 2);
70 1057896 : if (!sample)
71 0 : goto error;
72 1057896 : if (!bset)
73 0 : return NULL;
74 1057896 : isl_int_set_si(sample->block.data[0], 1);
75 :
76 1057896 : if (bset->n_eq > 0) {
77 0 : isl_assert(bset->ctx, bset->n_eq == 1, goto error);
78 0 : isl_assert(bset->ctx, bset->n_ineq == 0, goto error);
79 0 : if (isl_int_is_one(bset->eq[0][1]))
80 0 : isl_int_neg(sample->el[1], bset->eq[0][0]);
81 : else {
82 0 : isl_assert(bset->ctx, isl_int_is_negone(bset->eq[0][1]),
83 : goto error);
84 0 : isl_int_set(sample->el[1], bset->eq[0][0]);
85 : }
86 0 : isl_basic_set_free(bset);
87 0 : return sample;
88 : }
89 :
90 1057896 : isl_int_init(t);
91 1057896 : if (isl_int_is_one(bset->ineq[0][1]))
92 575824 : isl_int_neg(sample->block.data[1], bset->ineq[0][0]);
93 : else
94 482072 : isl_int_set(sample->block.data[1], bset->ineq[0][0]);
95 2085644 : for (i = 1; i < bset->n_ineq; ++i) {
96 1027748 : isl_seq_inner_product(sample->block.data,
97 1027748 : bset->ineq[i], 2, &t);
98 1027748 : if (isl_int_is_neg(t))
99 0 : break;
100 : }
101 1057896 : isl_int_clear(t);
102 1057896 : if (i < bset->n_ineq) {
103 0 : isl_vec_free(sample);
104 0 : return empty_sample(bset);
105 : }
106 :
107 1057896 : isl_basic_set_free(bset);
108 1057896 : return sample;
109 : error:
110 0 : isl_basic_set_free(bset);
111 0 : isl_vec_free(sample);
112 0 : return NULL;
113 : }
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 22116812 : 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 : {
125 : struct isl_mat *T;
126 : struct isl_vec *sample;
127 :
128 22116812 : if (!bset)
129 0 : return NULL;
130 :
131 22116812 : bset = isl_basic_set_remove_equalities(bset, &T, NULL);
132 22116812 : sample = recurse(bset);
133 22116812 : if (!sample || sample->size == 0)
134 10345793 : isl_mat_free(T);
135 : else
136 11771019 : sample = isl_mat_vec_product(T, sample);
137 22116812 : return sample;
138 : }
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 1123655 : static struct isl_mat *tab_equalities(struct isl_tab *tab)
145 : {
146 : int i, j;
147 : int n_eq;
148 : struct isl_mat *eq;
149 : struct isl_basic_set *bset;
150 :
151 1123655 : if (!tab)
152 0 : return NULL;
153 :
154 1123655 : bset = isl_tab_peek_bset(tab);
155 1123655 : isl_assert(tab->mat->ctx, bset, return NULL);
156 :
157 1123655 : n_eq = tab->n_var - tab->n_col + tab->n_dead;
158 1123655 : if (tab->empty || n_eq == 0)
159 0 : return isl_mat_alloc(tab->mat->ctx, 0, tab->n_var);
160 1123655 : if (n_eq == tab->n_var)
161 0 : return isl_mat_identity(tab->mat->ctx, tab->n_var);
162 :
163 1123655 : eq = isl_mat_alloc(tab->mat->ctx, n_eq, tab->n_var);
164 1123655 : if (!eq)
165 0 : return NULL;
166 16327912 : for (i = 0, j = 0; i < tab->n_con; ++i) {
167 15204257 : if (tab->con[i].is_row)
168 10096248 : continue;
169 5108009 : if (tab->con[i].index >= 0 && tab->con[i].index >= tab->n_dead)
170 3471518 : continue;
171 1636491 : if (i < bset->n_eq)
172 109 : isl_seq_cpy(eq->row[j], bset->eq[i] + 1, tab->n_var);
173 : else
174 3272764 : isl_seq_cpy(eq->row[j],
175 1636382 : bset->ineq[i - bset->n_eq] + 1, tab->n_var);
176 1636491 : ++j;
177 : }
178 1123655 : isl_assert(bset->ctx, j == n_eq, goto error);
179 1123655 : return eq;
180 : error:
181 0 : isl_mat_free(eq);
182 0 : return NULL;
183 : }
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 16559854 : static struct isl_mat *initial_basis(struct isl_tab *tab)
193 : {
194 : int n_eq;
195 : struct isl_mat *eq;
196 : struct isl_mat *Q;
197 :
198 16559854 : tab->n_unbounded = 0;
199 16559854 : tab->n_zero = n_eq = tab->n_var - tab->n_col + tab->n_dead;
200 16559854 : if (tab->empty || n_eq == 0 || n_eq == tab->n_var)
201 15436199 : return isl_mat_identity(tab->mat->ctx, 1 + tab->n_var);
202 :
203 1123655 : eq = tab_equalities(tab);
204 1123655 : eq = isl_mat_left_hermite(eq, 0, NULL, &Q);
205 1123655 : if (!eq)
206 0 : return NULL;
207 1123655 : isl_mat_free(eq);
208 :
209 1123655 : Q = isl_mat_lin_to_aff(Q);
210 1123655 : return Q;
211 : }
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 72811324 : static enum isl_lp_result compute_min(isl_ctx *ctx, struct isl_tab *tab,
220 : __isl_keep isl_vec *min, int level)
221 : {
222 218433972 : return isl_tab_min(tab, tab->basis->row[1 + level],
223 218433972 : ctx->one, &min->el[level], NULL, 0);
224 : }
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 70348156 : static enum isl_lp_result compute_max(isl_ctx *ctx, struct isl_tab *tab,
233 : __isl_keep isl_vec *max, int level)
234 : {
235 : enum isl_lp_result res;
236 70348156 : unsigned dim = tab->n_var;
237 :
238 70348156 : isl_seq_neg(tab->basis->row[1 + level] + 1,
239 70348156 : tab->basis->row[1 + level] + 1, dim);
240 211044468 : res = isl_tab_min(tab, tab->basis->row[1 + level],
241 211044468 : ctx->one, &max->el[level], NULL, 0);
242 70348156 : isl_seq_neg(tab->basis->row[1 + level] + 1,
243 70348156 : tab->basis->row[1 + level] + 1, dim);
244 70348156 : isl_int_neg(max->el[level], max->el[level]);
245 :
246 70348156 : return res;
247 : }
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 16002777 : 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 : {
275 : struct isl_tab_undo *snap;
276 : enum isl_lp_result res;
277 :
278 16002777 : snap = isl_tab_snap(tab);
279 :
280 : do {
281 60474470 : isl_int_add(tab->basis->row[1 + level][0],
282 : min->el[level], max->el[level]);
283 60474470 : isl_int_fdiv_q_ui(tab->basis->row[1 + level][0],
284 : tab->basis->row[1 + level][0], 2);
285 60474470 : isl_int_neg(tab->basis->row[1 + level][0],
286 : tab->basis->row[1 + level][0]);
287 60474470 : if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
288 0 : return -1;
289 60474470 : isl_int_set_si(tab->basis->row[1 + level][0], 0);
290 :
291 60474470 : if (++level >= tab->n_var - tab->n_unbounded)
292 8422653 : return 1;
293 52051817 : if (isl_tab_sample_is_integer(tab))
294 3643765 : return 1;
295 :
296 48408052 : res = compute_min(ctx, tab, min, level);
297 48408052 : if (res == isl_lp_error)
298 0 : return -1;
299 48408052 : if (res != isl_lp_ok)
300 0 : isl_die(ctx, isl_error_internal,
301 : "expecting bounded rational solution",
302 : return -1);
303 48408052 : res = compute_max(ctx, tab, max, level);
304 48408052 : if (res == isl_lp_error)
305 0 : return -1;
306 48408052 : if (res != isl_lp_ok)
307 0 : isl_die(ctx, isl_error_internal,
308 : "expecting bounded rational solution",
309 : return -1);
310 48408052 : } while (isl_int_le(min->el[level], max->el[level]));
311 :
312 3936359 : if (isl_tab_rollback(tab, snap) < 0)
313 0 : return -1;
314 :
315 3936359 : return 0;
316 : }
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 16841069 : struct isl_vec *isl_tab_sample(struct isl_tab *tab)
374 : {
375 : unsigned dim;
376 : unsigned gbr;
377 : struct isl_ctx *ctx;
378 : struct isl_vec *sample;
379 : struct isl_vec *min;
380 : struct isl_vec *max;
381 : enum isl_lp_result res;
382 : int level;
383 : int init;
384 : int reduced;
385 : struct isl_tab_undo **snap;
386 :
387 16841069 : if (!tab)
388 0 : return NULL;
389 16841069 : if (tab->empty)
390 198270 : return isl_vec_alloc(tab->mat->ctx, 0);
391 :
392 16642799 : if (!tab->basis)
393 16559854 : tab->basis = initial_basis(tab);
394 16642799 : if (!tab->basis)
395 0 : return NULL;
396 16642799 : isl_assert(tab->mat->ctx, tab->basis->n_row == tab->n_var + 1,
397 : return NULL);
398 16642799 : isl_assert(tab->mat->ctx, tab->basis->n_col == tab->n_var + 1,
399 : return NULL);
400 :
401 16642799 : ctx = tab->mat->ctx;
402 16642799 : dim = tab->n_var;
403 16642799 : gbr = ctx->opt->gbr;
404 :
405 16642799 : if (tab->n_unbounded == tab->n_var) {
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 : sample);
411 0 : return sample;
412 : }
413 :
414 16642799 : if (isl_tab_extend_cons(tab, dim + 1) < 0)
415 0 : return NULL;
416 :
417 16642799 : min = isl_vec_alloc(ctx, dim);
418 16642799 : max = isl_vec_alloc(ctx, dim);
419 16642799 : snap = isl_alloc_array(ctx, struct isl_tab_undo *, dim);
420 :
421 16642799 : if (!min || !max || !snap)
422 : goto error;
423 :
424 16642799 : level = 0;
425 16642799 : init = 1;
426 16642799 : reduced = 0;
427 :
428 44431626 : while (level >= 0) {
429 26698549 : if (init) {
430 : int choice;
431 :
432 24403272 : res = compute_min(ctx, tab, min, level);
433 24403272 : if (res == isl_lp_error)
434 0 : goto error;
435 24403272 : if (res != isl_lp_ok)
436 0 : isl_die(ctx, isl_error_internal,
437 : "expecting bounded rational solution",
438 : goto error);
439 24403272 : if (isl_tab_sample_is_integer(tab))
440 2463168 : break;
441 21940104 : res = compute_max(ctx, tab, max, level);
442 21940104 : if (res == isl_lp_error)
443 0 : goto error;
444 21940104 : if (res != isl_lp_ok)
445 0 : isl_die(ctx, isl_error_internal,
446 : "expecting bounded rational solution",
447 : goto error);
448 21940104 : if (isl_tab_sample_is_integer(tab))
449 972471 : break;
450 20967633 : choice = isl_int_lt(min->el[level], max->el[level]);
451 20967633 : if (choice) {
452 : int g;
453 16002777 : g = greedy_search(ctx, tab, min, max, level);
454 16002777 : if (g < 0)
455 0 : goto error;
456 16002777 : if (g)
457 12066418 : break;
458 : }
459 12262174 : if (!reduced && choice &&
460 3360959 : ctx->opt->gbr != ISL_GBR_NEVER) {
461 : unsigned gbr_only_first;
462 3360959 : if (ctx->opt->gbr == ISL_GBR_ONCE)
463 0 : ctx->opt->gbr = ISL_GBR_NEVER;
464 3360959 : tab->n_zero = level;
465 3360959 : gbr_only_first = ctx->opt->gbr_only_first;
466 6721918 : ctx->opt->gbr_only_first =
467 3360959 : ctx->opt->gbr == ISL_GBR_ALWAYS;
468 3360959 : tab = isl_tab_compute_reduced_basis(tab);
469 3360959 : ctx->opt->gbr_only_first = gbr_only_first;
470 3360959 : if (!tab || !tab->basis)
471 : goto error;
472 3360959 : reduced = 1;
473 3360959 : continue;
474 : }
475 5540256 : reduced = 0;
476 5540256 : snap[level] = isl_tab_snap(tab);
477 : } else
478 2295277 : isl_int_add_ui(min->el[level], min->el[level], 1);
479 :
480 7835533 : if (isl_int_gt(min->el[level], max->el[level])) {
481 3385555 : level--;
482 3385555 : init = 0;
483 3385555 : if (level >= 0)
484 2295277 : if (isl_tab_rollback(tab, snap[level]) < 0)
485 0 : goto error;
486 3385555 : continue;
487 : }
488 4449978 : isl_int_neg(tab->basis->row[1 + level][0], min->el[level]);
489 4449978 : if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
490 0 : goto error;
491 4449978 : isl_int_set_si(tab->basis->row[1 + level][0], 0);
492 4449978 : if (level + tab->n_unbounded < dim - 1) {
493 4399514 : ++level;
494 4399514 : init = 1;
495 4399514 : continue;
496 : }
497 50464 : break;
498 : }
499 :
500 16642799 : if (level >= 0) {
501 15552521 : sample = isl_tab_get_sample_value(tab);
502 15552521 : if (!sample)
503 0 : goto error;
504 15552521 : if (tab->n_unbounded && !isl_int_is_one(sample->el[0])) {
505 0 : sample = isl_mat_vec_product(isl_mat_copy(tab->basis),
506 : sample);
507 0 : sample = isl_vec_ceil(sample);
508 0 : sample = isl_mat_vec_inverse_product(
509 0 : isl_mat_copy(tab->basis), sample);
510 : }
511 : } else
512 1090278 : sample = isl_vec_alloc(ctx, 0);
513 :
514 16642799 : ctx->opt->gbr = gbr;
515 16642799 : isl_vec_free(min);
516 16642799 : isl_vec_free(max);
517 16642799 : free(snap);
518 16642799 : return sample;
519 : 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 : }
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 129106 : static __isl_give isl_vec *factored_sample(__isl_take isl_basic_set *bset,
533 : __isl_take isl_factorizer *f)
534 : {
535 : int i, n;
536 129106 : isl_vec *sample = NULL;
537 : isl_ctx *ctx;
538 : unsigned nparam;
539 : unsigned nvar;
540 :
541 129106 : ctx = isl_basic_set_get_ctx(bset);
542 129106 : if (!ctx)
543 0 : goto error;
544 :
545 129106 : nparam = isl_basic_set_dim(bset, isl_dim_param);
546 129106 : nvar = isl_basic_set_dim(bset, isl_dim_set);
547 :
548 129106 : sample = isl_vec_alloc(ctx, 1 + isl_basic_set_total_dim(bset));
549 129106 : if (!sample)
550 0 : goto error;
551 129106 : isl_int_set_si(sample->el[0], 1);
552 :
553 129106 : bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
554 :
555 330407 : for (i = 0, n = 0; i < f->n_group; ++i) {
556 : isl_basic_set *bset_i;
557 : isl_vec *sample_i;
558 :
559 238314 : bset_i = isl_basic_set_copy(bset);
560 476628 : bset_i = isl_basic_set_drop_constraints_involving(bset_i,
561 476628 : nparam + n + f->len[i], nvar - n - f->len[i]);
562 238314 : bset_i = isl_basic_set_drop_constraints_involving(bset_i,
563 : nparam, n);
564 476628 : bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
565 476628 : n + f->len[i], nvar - n - f->len[i]);
566 238314 : bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
567 :
568 238314 : sample_i = sample_bounded(bset_i);
569 238314 : if (!sample_i)
570 0 : goto error;
571 238314 : if (sample_i->size == 0) {
572 37013 : isl_basic_set_free(bset);
573 37013 : isl_factorizer_free(f);
574 37013 : isl_vec_free(sample);
575 37013 : return sample_i;
576 : }
577 402602 : isl_seq_cpy(sample->el + 1 + nparam + n,
578 402602 : sample_i->el + 1, f->len[i]);
579 201301 : isl_vec_free(sample_i);
580 :
581 201301 : n += f->len[i];
582 : }
583 :
584 92093 : f->morph = isl_morph_inverse(f->morph);
585 92093 : sample = isl_morph_vec(isl_morph_copy(f->morph), sample);
586 :
587 92093 : isl_basic_set_free(bset);
588 92093 : isl_factorizer_free(f);
589 92093 : return sample;
590 : 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 : }
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 55523765 : static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset)
605 : {
606 : unsigned dim;
607 : struct isl_vec *sample;
608 55523765 : struct isl_tab *tab = NULL;
609 : isl_factorizer *f;
610 :
611 55523765 : if (!bset)
612 0 : return NULL;
613 :
614 55523765 : if (isl_basic_set_plain_is_empty(bset))
615 0 : return empty_sample(bset);
616 :
617 55523765 : dim = isl_basic_set_total_dim(bset);
618 55523765 : if (dim == 0)
619 10569884 : return zero_sample(bset);
620 44953881 : if (dim == 1)
621 995698 : return interval_sample(bset);
622 43958183 : if (bset->n_eq > 0)
623 0 : return sample_eq(bset, sample_bounded);
624 :
625 43958183 : f = isl_basic_set_factorizer(bset);
626 43958183 : if (!f)
627 0 : goto error;
628 43958183 : if (f->n_group != 0)
629 129106 : return factored_sample(bset, f);
630 43829077 : isl_factorizer_free(f);
631 :
632 43829077 : tab = isl_tab_from_basic_set(bset, 1);
633 43829077 : if (tab && tab->empty) {
634 27237109 : isl_tab_free(tab);
635 27237109 : ISL_F_SET(bset, ISL_BASIC_SET_EMPTY);
636 27237109 : sample = isl_vec_alloc(isl_basic_set_get_ctx(bset), 0);
637 27237109 : isl_basic_set_free(bset);
638 27237109 : return sample;
639 : }
640 :
641 16591968 : if (!ISL_F_ISSET(bset, ISL_BASIC_SET_NO_IMPLICIT))
642 16547718 : if (isl_tab_detect_implicit_equalities(tab) < 0)
643 0 : goto error;
644 :
645 16591968 : sample = isl_tab_sample(tab);
646 16591968 : if (!sample)
647 0 : goto error;
648 :
649 16591968 : if (sample->size > 0) {
650 15388094 : isl_vec_free(bset->sample);
651 15388094 : bset->sample = isl_vec_copy(sample);
652 : }
653 :
654 16591968 : isl_basic_set_free(bset);
655 16591968 : isl_tab_free(tab);
656 16591968 : return sample;
657 : error:
658 0 : isl_basic_set_free(bset);
659 0 : isl_tab_free(tab);
660 0 : return NULL;
661 : }
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 14447846 : static __isl_give isl_basic_set *plug_in(__isl_take isl_basic_set *bset,
676 : __isl_take isl_vec *sample)
677 : {
678 : int i;
679 : unsigned total;
680 : struct isl_mat *T;
681 :
682 14447846 : if (!bset || !sample)
683 : goto error;
684 :
685 14447846 : total = isl_basic_set_total_dim(bset);
686 14447846 : T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1));
687 14447846 : if (!T)
688 0 : goto error;
689 :
690 40112358 : for (i = 0; i < sample->size; ++i) {
691 25664512 : isl_int_set(T->row[i][0], sample->el[i]);
692 25664512 : isl_seq_clr(T->row[i] + 1, T->n_col - 1);
693 : }
694 76221758 : for (i = 0; i < T->n_col - 1; ++i) {
695 61773912 : isl_seq_clr(T->row[sample->size + i], T->n_col);
696 61773912 : isl_int_set_si(T->row[sample->size + i][1 + i], 1);
697 : }
698 14447846 : isl_vec_free(sample);
699 :
700 14447846 : bset = isl_basic_set_preimage(bset, T);
701 14447846 : return bset;
702 : error:
703 0 : isl_basic_set_free(bset);
704 0 : isl_vec_free(sample);
705 0 : return NULL;
706 : }
707 :
708 : /* Given a basic set "bset", return any (possibly non-integer) point
709 : * in the basic set.
710 : */
711 22151852 : static __isl_give isl_vec *rational_sample(__isl_take isl_basic_set *bset)
712 : {
713 : struct isl_tab *tab;
714 : struct isl_vec *sample;
715 :
716 22151852 : if (!bset)
717 0 : return NULL;
718 :
719 22151852 : tab = isl_tab_from_basic_set(bset, 0);
720 22151852 : sample = isl_tab_get_sample_value(tab);
721 22151852 : isl_tab_free(tab);
722 :
723 22151852 : isl_basic_set_free(bset);
724 :
725 22151852 : return sample;
726 : }
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 7704006 : static __isl_give isl_basic_set *shift_cone(__isl_take isl_basic_set *cone,
762 : __isl_take isl_vec *vec)
763 : {
764 : int i, j, k;
765 : unsigned total;
766 :
767 7704006 : struct isl_basic_set *shift = NULL;
768 :
769 7704006 : if (!cone || !vec)
770 : goto error;
771 :
772 7704006 : isl_assert(cone->ctx, cone->n_eq == 0, goto error);
773 :
774 7704006 : total = isl_basic_set_total_dim(cone);
775 :
776 7704006 : shift = isl_basic_set_alloc_space(isl_basic_set_get_space(cone),
777 : 0, 0, cone->n_ineq);
778 :
779 48743050 : for (i = 0; i < cone->n_ineq; ++i) {
780 41039044 : k = isl_basic_set_alloc_inequality(shift);
781 41039044 : if (k < 0)
782 0 : goto error;
783 41039044 : isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total);
784 41039044 : isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total,
785 41039044 : &shift->ineq[k][0]);
786 41039044 : isl_int_cdiv_q(shift->ineq[k][0],
787 : shift->ineq[k][0], vec->el[0]);
788 41039044 : isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]);
789 253120520 : for (j = 0; j < total; ++j) {
790 212081476 : if (isl_int_is_nonneg(shift->ineq[k][1 + j]))
791 155296258 : continue;
792 56785218 : isl_int_add(shift->ineq[k][0],
793 : shift->ineq[k][0], shift->ineq[k][1 + j]);
794 : }
795 : }
796 :
797 7704006 : isl_basic_set_free(cone);
798 7704006 : isl_vec_free(vec);
799 :
800 7704006 : return isl_basic_set_finalize(shift);
801 : 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 : }
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 14447846 : 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 : {
822 : unsigned total;
823 :
824 14447846 : if (!vec || !cone || !U)
825 : goto error;
826 :
827 14447846 : isl_assert(vec->ctx, vec->size != 0, goto error);
828 14447846 : if (isl_int_is_one(vec->el[0])) {
829 6743840 : isl_mat_free(U);
830 6743840 : isl_basic_set_free(cone);
831 6743840 : return vec;
832 : }
833 :
834 7704006 : total = isl_basic_set_total_dim(cone);
835 7704006 : cone = isl_basic_set_preimage(cone, U);
836 7704006 : cone = isl_basic_set_remove_dims(cone, isl_dim_set,
837 7704006 : 0, total - (vec->size - 1));
838 :
839 7704006 : cone = shift_cone(cone, vec);
840 :
841 7704006 : vec = rational_sample(cone);
842 7704006 : vec = isl_vec_ceil(vec);
843 7704006 : return vec;
844 : 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 : }
850 :
851 : /* Concatenate two integer vectors, i.e., two vectors with denominator
852 : * (stored in element 0) equal to 1.
853 : */
854 14447846 : static __isl_give isl_vec *vec_concat(__isl_take isl_vec *vec1,
855 : __isl_take isl_vec *vec2)
856 : {
857 : struct isl_vec *vec;
858 :
859 14447846 : if (!vec1 || !vec2)
860 : goto error;
861 14447846 : isl_assert(vec1->ctx, vec1->size > 0, goto error);
862 14447846 : isl_assert(vec2->ctx, vec2->size > 0, goto error);
863 14447846 : isl_assert(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error);
864 14447846 : isl_assert(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error);
865 :
866 14447846 : vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1);
867 14447846 : if (!vec)
868 0 : goto error;
869 :
870 14447846 : isl_seq_cpy(vec->el, vec1->el, vec1->size);
871 14447846 : isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1);
872 :
873 14447846 : isl_vec_free(vec1);
874 14447846 : isl_vec_free(vec2);
875 :
876 14447846 : return vec;
877 : error:
878 0 : isl_vec_free(vec1);
879 0 : isl_vec_free(vec2);
880 0 : return NULL;
881 : }
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 16016813 : __isl_give isl_vec *isl_basic_set_sample_with_cone(
921 : __isl_take isl_basic_set *bset, __isl_take isl_basic_set *cone)
922 : {
923 : unsigned total;
924 : unsigned cone_dim;
925 : struct isl_mat *M, *U;
926 : struct isl_vec *sample;
927 : struct isl_vec *cone_sample;
928 : struct isl_ctx *ctx;
929 : struct isl_basic_set *bounded;
930 :
931 16016813 : if (!bset || !cone)
932 : goto error;
933 :
934 16016813 : ctx = isl_basic_set_get_ctx(bset);
935 16016813 : total = isl_basic_set_total_dim(cone);
936 16016813 : cone_dim = total - cone->n_eq;
937 :
938 16016813 : M = isl_mat_sub_alloc6(ctx, cone->eq, 0, cone->n_eq, 1, total);
939 16016813 : M = isl_mat_left_hermite(M, 0, &U, NULL);
940 16016813 : if (!M)
941 0 : goto error;
942 16016813 : isl_mat_free(M);
943 :
944 16016813 : U = isl_mat_lin_to_aff(U);
945 16016813 : bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
946 :
947 16016813 : bounded = isl_basic_set_copy(bset);
948 16016813 : bounded = isl_basic_set_drop_constraints_involving(bounded,
949 : total - cone_dim, cone_dim);
950 16016813 : bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim);
951 16016813 : sample = sample_bounded(bounded);
952 16016813 : if (!sample || sample->size == 0) {
953 1568967 : isl_basic_set_free(bset);
954 1568967 : isl_basic_set_free(cone);
955 1568967 : isl_mat_free(U);
956 1568967 : return sample;
957 : }
958 14447846 : bset = plug_in(bset, isl_vec_copy(sample));
959 14447846 : cone_sample = rational_sample(bset);
960 14447846 : cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(U));
961 14447846 : sample = vec_concat(sample, cone_sample);
962 14447846 : sample = isl_mat_vec_product(U, sample);
963 14447846 : return sample;
964 : error:
965 0 : isl_basic_set_free(cone);
966 0 : isl_basic_set_free(bset);
967 0 : return NULL;
968 : }
969 :
970 0 : static void vec_sum_of_neg(struct isl_vec *v, isl_int *s)
971 : {
972 : int i;
973 :
974 0 : isl_int_set_si(*s, 0);
975 :
976 0 : for (i = 0; i < v->size; ++i)
977 0 : if (isl_int_is_neg(v->el[i]))
978 0 : isl_int_add(*s, *s, v->el[i]);
979 0 : }
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 0 : static int tab_shift_cone(struct isl_tab *tab,
999 : struct isl_tab *tab_cone, struct isl_mat *U)
1000 : {
1001 : int i;
1002 : isl_int v;
1003 0 : struct isl_basic_set *bset = NULL;
1004 :
1005 0 : if (tab && tab->n_unbounded == 0) {
1006 0 : isl_mat_free(U);
1007 0 : return 0;
1008 : }
1009 0 : isl_int_init(v);
1010 0 : if (!tab || !tab_cone || !U)
1011 : goto error;
1012 0 : bset = isl_tab_peek_bset(tab_cone);
1013 0 : U = isl_mat_drop_cols(U, 0, tab->n_var - tab->n_unbounded);
1014 0 : for (i = 0; i < bset->n_ineq; ++i) {
1015 : int ok;
1016 0 : struct isl_vec *row = NULL;
1017 0 : if (isl_tab_is_equality(tab_cone, tab_cone->n_eq + i))
1018 0 : continue;
1019 0 : row = isl_vec_alloc(bset->ctx, tab_cone->n_var);
1020 0 : if (!row)
1021 0 : goto error;
1022 0 : isl_seq_cpy(row->el, bset->ineq[i] + 1, tab_cone->n_var);
1023 0 : row = isl_vec_mat_product(row, isl_mat_copy(U));
1024 0 : if (!row)
1025 0 : goto error;
1026 0 : vec_sum_of_neg(row, &v);
1027 0 : isl_vec_free(row);
1028 0 : if (isl_int_is_zero(v))
1029 0 : continue;
1030 0 : if (isl_tab_extend_cons(tab, 1) < 0)
1031 0 : goto error;
1032 0 : isl_int_add(bset->ineq[i][0], bset->ineq[i][0], v);
1033 0 : ok = isl_tab_add_ineq(tab, bset->ineq[i]) >= 0;
1034 0 : isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], v);
1035 0 : if (!ok)
1036 0 : goto error;
1037 : }
1038 :
1039 0 : isl_mat_free(U);
1040 0 : isl_int_clear(v);
1041 0 : return 0;
1042 : error:
1043 0 : isl_mat_free(U);
1044 0 : isl_int_clear(v);
1045 0 : return -1;
1046 : }
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 0 : int isl_tab_set_initial_basis_with_cone(struct isl_tab *tab,
1063 : struct isl_tab *tab_cone)
1064 : {
1065 : struct isl_mat *eq;
1066 : struct isl_mat *cone_eq;
1067 : struct isl_mat *U, *Q;
1068 :
1069 0 : if (!tab || !tab_cone)
1070 0 : return -1;
1071 :
1072 0 : if (tab_cone->n_col == tab_cone->n_dead) {
1073 0 : tab->basis = initial_basis(tab);
1074 0 : return tab->basis ? 0 : -1;
1075 : }
1076 :
1077 0 : eq = tab_equalities(tab);
1078 0 : if (!eq)
1079 0 : return -1;
1080 0 : tab->n_zero = eq->n_row;
1081 0 : cone_eq = tab_equalities(tab_cone);
1082 0 : eq = isl_mat_concat(eq, cone_eq);
1083 0 : if (!eq)
1084 0 : return -1;
1085 0 : tab->n_unbounded = tab->n_var - (eq->n_row - tab->n_zero);
1086 0 : eq = isl_mat_left_hermite(eq, 0, &U, &Q);
1087 0 : if (!eq)
1088 0 : return -1;
1089 0 : isl_mat_free(eq);
1090 0 : tab->basis = isl_mat_lin_to_aff(Q);
1091 0 : if (tab_shift_cone(tab, tab_cone, U) < 0)
1092 0 : return -1;
1093 0 : if (!tab->basis)
1094 0 : return -1;
1095 0 : return 0;
1096 : }
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 55285451 : static __isl_give isl_vec *gbr_sample(__isl_take isl_basic_set *bset)
1105 : {
1106 : unsigned dim;
1107 : struct isl_basic_set *cone;
1108 :
1109 55285451 : dim = isl_basic_set_total_dim(bset);
1110 :
1111 55285451 : cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
1112 55285451 : if (!cone)
1113 0 : goto error;
1114 :
1115 55285451 : if (cone->n_eq < dim)
1116 16016813 : return isl_basic_set_sample_with_cone(bset, cone);
1117 :
1118 39268638 : isl_basic_set_free(cone);
1119 39268638 : return sample_bounded(bset);
1120 : error:
1121 0 : isl_basic_set_free(bset);
1122 0 : return NULL;
1123 : }
1124 :
1125 77648185 : static __isl_give isl_vec *basic_set_sample(__isl_take isl_basic_set *bset,
1126 : int bounded)
1127 : {
1128 : struct isl_ctx *ctx;
1129 : unsigned dim;
1130 77648185 : if (!bset)
1131 0 : return NULL;
1132 :
1133 77648185 : ctx = bset->ctx;
1134 77648185 : if (isl_basic_set_plain_is_empty(bset))
1135 166168 : return empty_sample(bset);
1136 :
1137 77482017 : dim = isl_basic_set_n_dim(bset);
1138 77482017 : isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
1139 77482017 : isl_assert(ctx, bset->n_div == 0, goto error);
1140 :
1141 77482017 : if (bset->sample && bset->sample->size == 1 + dim) {
1142 0 : int contains = isl_basic_set_contains(bset, bset->sample);
1143 0 : if (contains < 0)
1144 0 : goto error;
1145 0 : if (contains) {
1146 0 : struct isl_vec *sample = isl_vec_copy(bset->sample);
1147 0 : isl_basic_set_free(bset);
1148 0 : return sample;
1149 : }
1150 : }
1151 77482017 : isl_vec_free(bset->sample);
1152 77482017 : bset->sample = NULL;
1153 :
1154 77482017 : if (bset->n_eq > 0)
1155 22116812 : return sample_eq(bset, bounded ? isl_basic_set_sample_bounded
1156 : : isl_basic_set_sample_vec);
1157 55365205 : if (dim == 0)
1158 16473 : return zero_sample(bset);
1159 55348732 : if (dim == 1)
1160 63281 : return interval_sample(bset);
1161 :
1162 55285451 : return bounded ? sample_bounded(bset) : gbr_sample(bset);
1163 : error:
1164 0 : isl_basic_set_free(bset);
1165 0 : return NULL;
1166 : }
1167 :
1168 77648185 : __isl_give isl_vec *isl_basic_set_sample_vec(__isl_take isl_basic_set *bset)
1169 : {
1170 77648185 : return basic_set_sample(bset, 0);
1171 : }
1172 :
1173 : /* Compute an integer sample in "bset", where the caller guarantees
1174 : * that "bset" is bounded.
1175 : */
1176 0 : __isl_give isl_vec *isl_basic_set_sample_bounded(__isl_take isl_basic_set *bset)
1177 : {
1178 0 : return basic_set_sample(bset, 1);
1179 : }
1180 :
1181 2048330 : __isl_give isl_basic_set *isl_basic_set_from_vec(__isl_take isl_vec *vec)
1182 : {
1183 : int i;
1184 : int k;
1185 2048330 : struct isl_basic_set *bset = NULL;
1186 : struct isl_ctx *ctx;
1187 : unsigned dim;
1188 :
1189 2048330 : if (!vec)
1190 0 : return NULL;
1191 2048330 : ctx = vec->ctx;
1192 2048330 : isl_assert(ctx, vec->size != 0, goto error);
1193 :
1194 2048330 : bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
1195 2048330 : if (!bset)
1196 0 : goto error;
1197 2048330 : dim = isl_basic_set_n_dim(bset);
1198 13203600 : for (i = dim - 1; i >= 0; --i) {
1199 11155270 : k = isl_basic_set_alloc_equality(bset);
1200 11155270 : if (k < 0)
1201 0 : goto error;
1202 11155270 : isl_seq_clr(bset->eq[k], 1 + dim);
1203 11155270 : isl_int_neg(bset->eq[k][0], vec->el[1 + i]);
1204 11155270 : isl_int_set(bset->eq[k][1 + i], vec->el[0]);
1205 : }
1206 2048330 : bset->sample = vec;
1207 :
1208 2048330 : return bset;
1209 : error:
1210 0 : isl_basic_set_free(bset);
1211 0 : isl_vec_free(vec);
1212 0 : return NULL;
1213 : }
1214 :
1215 0 : __isl_give isl_basic_map *isl_basic_map_sample(__isl_take isl_basic_map *bmap)
1216 : {
1217 : struct isl_basic_set *bset;
1218 : struct isl_vec *sample_vec;
1219 :
1220 0 : bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
1221 0 : sample_vec = isl_basic_set_sample_vec(bset);
1222 0 : if (!sample_vec)
1223 0 : goto error;
1224 0 : if (sample_vec->size == 0) {
1225 0 : isl_vec_free(sample_vec);
1226 0 : return isl_basic_map_set_to_empty(bmap);
1227 : }
1228 0 : isl_vec_free(bmap->sample);
1229 0 : bmap->sample = isl_vec_copy(sample_vec);
1230 0 : bset = isl_basic_set_from_vec(sample_vec);
1231 0 : return isl_basic_map_overlying_set(bset, bmap);
1232 : error:
1233 0 : isl_basic_map_free(bmap);
1234 0 : return NULL;
1235 : }
1236 :
1237 0 : __isl_give isl_basic_set *isl_basic_set_sample(__isl_take isl_basic_set *bset)
1238 : {
1239 0 : return isl_basic_map_sample(bset);
1240 : }
1241 :
1242 0 : __isl_give isl_basic_map *isl_map_sample(__isl_take isl_map *map)
1243 : {
1244 : int i;
1245 0 : isl_basic_map *sample = NULL;
1246 :
1247 0 : if (!map)
1248 0 : goto error;
1249 :
1250 0 : for (i = 0; i < map->n; ++i) {
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 (!ISL_F_ISSET(sample, ISL_BASIC_MAP_EMPTY))
1255 0 : break;
1256 0 : isl_basic_map_free(sample);
1257 : }
1258 0 : if (i == map->n)
1259 0 : sample = isl_basic_map_empty(isl_map_get_space(map));
1260 0 : isl_map_free(map);
1261 0 : return sample;
1262 : error:
1263 0 : isl_map_free(map);
1264 0 : return NULL;
1265 : }
1266 :
1267 0 : __isl_give isl_basic_set *isl_set_sample(__isl_take isl_set *set)
1268 : {
1269 0 : return bset_from_bmap(isl_map_sample(set_to_map(set)));
1270 : }
1271 :
1272 0 : __isl_give isl_point *isl_basic_set_sample_point(__isl_take isl_basic_set *bset)
1273 : {
1274 : isl_vec *vec;
1275 : isl_space *dim;
1276 :
1277 0 : dim = isl_basic_set_get_space(bset);
1278 0 : bset = isl_basic_set_underlying_set(bset);
1279 0 : vec = isl_basic_set_sample_vec(bset);
1280 :
1281 0 : return isl_point_alloc(dim, vec);
1282 : }
1283 :
1284 0 : __isl_give isl_point *isl_set_sample_point(__isl_take isl_set *set)
1285 : {
1286 : int i;
1287 : isl_point *pnt;
1288 :
1289 0 : if (!set)
1290 0 : return NULL;
1291 :
1292 0 : for (i = 0; i < set->n; ++i) {
1293 0 : pnt = isl_basic_set_sample_point(isl_basic_set_copy(set->p[i]));
1294 0 : if (!pnt)
1295 0 : goto error;
1296 0 : if (!isl_point_is_void(pnt))
1297 0 : break;
1298 0 : isl_point_free(pnt);
1299 : }
1300 0 : if (i == set->n)
1301 0 : pnt = isl_point_void(isl_set_get_space(set));
1302 :
1303 0 : isl_set_free(set);
1304 0 : return pnt;
1305 : error:
1306 0 : isl_set_free(set);
1307 0 : return NULL;
1308 : }
|
