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 497843 : static __isl_give isl_vec *empty_sample(__isl_take isl_basic_set *bset)
28 : {
29 : struct isl_vec *vec;
30 :
31 497843 : vec = isl_vec_alloc(bset->ctx, 0);
32 497843 : isl_basic_set_free(bset);
33 497843 : 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 56296302 : static __isl_give isl_vec *zero_sample(__isl_take isl_basic_set *bset)
41 : {
42 : unsigned dim;
43 : struct isl_vec *sample;
44 :
45 56296302 : dim = isl_basic_set_total_dim(bset);
46 56296302 : sample = isl_vec_alloc(bset->ctx, 1 + dim);
47 56296302 : if (sample) {
48 56296302 : isl_int_set_si(sample->el[0], 1);
49 56296302 : isl_seq_clr(sample->el + 1, dim);
50 : }
51 56296302 : isl_basic_set_free(bset);
52 56296302 : return sample;
53 : }
54 :
55 3050448 : 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 3050448 : bset = isl_basic_set_simplify(bset);
62 3050448 : if (!bset)
63 0 : return NULL;
64 3050448 : if (isl_basic_set_plain_is_empty(bset))
65 0 : return empty_sample(bset);
66 3050448 : if (bset->n_eq == 0 && bset->n_ineq == 0)
67 12052 : return zero_sample(bset);
68 :
69 3038396 : sample = isl_vec_alloc(bset->ctx, 2);
70 3038396 : if (!sample)
71 0 : goto error;
72 3038396 : if (!bset)
73 0 : return NULL;
74 3038396 : isl_int_set_si(sample->block.data[0], 1);
75 :
76 3038396 : 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 3038396 : isl_int_init(t);
91 3038396 : if (isl_int_is_one(bset->ineq[0][1]))
92 1566889 : isl_int_neg(sample->block.data[1], bset->ineq[0][0]);
93 : else
94 1471507 : isl_int_set(sample->block.data[1], bset->ineq[0][0]);
95 4055411 : for (i = 1; i < bset->n_ineq; ++i) {
96 1017015 : isl_seq_inner_product(sample->block.data,
97 1017015 : bset->ineq[i], 2, &t);
98 1017015 : if (isl_int_is_neg(t))
99 0 : break;
100 : }
101 3038396 : isl_int_clear(t);
102 3038396 : if (i < bset->n_ineq) {
103 0 : isl_vec_free(sample);
104 0 : return empty_sample(bset);
105 : }
106 :
107 3038396 : isl_basic_set_free(bset);
108 3038396 : 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 318523352 : 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 318523352 : if (!bset)
129 0 : return NULL;
130 :
131 318523352 : bset = isl_basic_set_remove_equalities(bset, &T, NULL);
132 318523352 : sample = recurse(bset);
133 318523352 : if (!sample || sample->size == 0)
134 295195710 : isl_mat_free(T);
135 : else
136 23327642 : sample = isl_mat_vec_product(T, sample);
137 318523352 : 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 1105670 : 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 1105670 : if (!tab)
152 0 : return NULL;
153 :
154 1105670 : bset = isl_tab_peek_bset(tab);
155 1105670 : isl_assert(tab->mat->ctx, bset, return NULL);
156 :
157 1105670 : n_eq = tab->n_var - tab->n_col + tab->n_dead;
158 1105670 : if (tab->empty || n_eq == 0)
159 0 : return isl_mat_alloc(tab->mat->ctx, 0, tab->n_var);
160 1105670 : if (n_eq == tab->n_var)
161 0 : return isl_mat_identity(tab->mat->ctx, tab->n_var);
162 :
163 1105670 : eq = isl_mat_alloc(tab->mat->ctx, n_eq, tab->n_var);
164 1105670 : if (!eq)
165 0 : return NULL;
166 16050466 : for (i = 0, j = 0; i < tab->n_con; ++i) {
167 14944796 : if (tab->con[i].is_row)
168 9922293 : continue;
169 5022503 : if (tab->con[i].index >= 0 && tab->con[i].index >= tab->n_dead)
170 3408688 : continue;
171 1613815 : if (i < bset->n_eq)
172 109 : isl_seq_cpy(eq->row[j], bset->eq[i] + 1, tab->n_var);
173 : else
174 3227412 : isl_seq_cpy(eq->row[j],
175 1613706 : bset->ineq[i - bset->n_eq] + 1, tab->n_var);
176 1613815 : ++j;
177 : }
178 1105670 : isl_assert(bset->ctx, j == n_eq, goto error);
179 1105670 : 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 16667022 : 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 16667022 : tab->n_unbounded = 0;
199 16667022 : tab->n_zero = n_eq = tab->n_var - tab->n_col + tab->n_dead;
200 16667022 : if (tab->empty || n_eq == 0 || n_eq == tab->n_var)
201 15561352 : return isl_mat_identity(tab->mat->ctx, 1 + tab->n_var);
202 :
203 1105670 : eq = tab_equalities(tab);
204 1105670 : eq = isl_mat_left_hermite(eq, 0, NULL, &Q);
205 1105670 : if (!eq)
206 0 : return NULL;
207 1105670 : isl_mat_free(eq);
208 :
209 1105670 : Q = isl_mat_lin_to_aff(Q);
210 1105670 : 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 71903652 : static enum isl_lp_result compute_min(isl_ctx *ctx, struct isl_tab *tab,
220 : __isl_keep isl_vec *min, int level)
221 : {
222 215710956 : return isl_tab_min(tab, tab->basis->row[1 + level],
223 215710956 : 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 69069068 : 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 69069068 : unsigned dim = tab->n_var;
237 :
238 69069068 : isl_seq_neg(tab->basis->row[1 + level] + 1,
239 69069068 : tab->basis->row[1 + level] + 1, dim);
240 207207204 : res = isl_tab_min(tab, tab->basis->row[1 + level],
241 207207204 : ctx->one, &max->el[level], NULL, 0);
242 69069068 : isl_seq_neg(tab->basis->row[1 + level] + 1,
243 69069068 : tab->basis->row[1 + level] + 1, dim);
244 69069068 : isl_int_neg(max->el[level], max->el[level]);
245 :
246 69069068 : 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 15707993 : 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 15707993 : snap = isl_tab_snap(tab);
279 :
280 : do {
281 59364401 : isl_int_add(tab->basis->row[1 + level][0],
282 : min->el[level], max->el[level]);
283 59364401 : isl_int_fdiv_q_ui(tab->basis->row[1 + level][0],
284 : tab->basis->row[1 + level][0], 2);
285 59364401 : isl_int_neg(tab->basis->row[1 + level][0],
286 : tab->basis->row[1 + level][0]);
287 59364401 : if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
288 0 : return -1;
289 59364401 : isl_int_set_si(tab->basis->row[1 + level][0], 0);
290 :
291 59364401 : if (++level >= tab->n_var - tab->n_unbounded)
292 8266949 : return 1;
293 51097452 : if (isl_tab_sample_is_integer(tab))
294 3574162 : return 1;
295 :
296 47523290 : res = compute_min(ctx, tab, min, level);
297 47523290 : if (res == isl_lp_error)
298 0 : return -1;
299 47523290 : if (res != isl_lp_ok)
300 0 : isl_die(ctx, isl_error_internal,
301 : "expecting bounded rational solution",
302 : return -1);
303 47523290 : res = compute_max(ctx, tab, max, level);
304 47523290 : if (res == isl_lp_error)
305 0 : return -1;
306 47523290 : if (res != isl_lp_ok)
307 0 : isl_die(ctx, isl_error_internal,
308 : "expecting bounded rational solution",
309 : return -1);
310 47523290 : } while (isl_int_le(min->el[level], max->el[level]));
311 :
312 3866882 : if (isl_tab_rollback(tab, snap) < 0)
313 0 : return -1;
314 :
315 3866882 : 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 17015310 : 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 17015310 : if (!tab)
388 0 : return NULL;
389 17015310 : if (tab->empty)
390 260940 : return isl_vec_alloc(tab->mat->ctx, 0);
391 :
392 16754370 : if (!tab->basis)
393 16667022 : tab->basis = initial_basis(tab);
394 16754370 : if (!tab->basis)
395 0 : return NULL;
396 16754370 : isl_assert(tab->mat->ctx, tab->basis->n_row == tab->n_var + 1,
397 : return NULL);
398 16754370 : isl_assert(tab->mat->ctx, tab->basis->n_col == tab->n_var + 1,
399 : return NULL);
400 :
401 16754370 : ctx = tab->mat->ctx;
402 16754370 : dim = tab->n_var;
403 16754370 : gbr = ctx->opt->gbr;
404 :
405 16754370 : 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 16754370 : if (isl_tab_extend_cons(tab, dim + 1) < 0)
415 0 : return NULL;
416 :
417 16754370 : min = isl_vec_alloc(ctx, dim);
418 16754370 : max = isl_vec_alloc(ctx, dim);
419 16754370 : snap = isl_alloc_array(ctx, struct isl_tab_undo *, dim);
420 :
421 16754370 : if (!min || !max || !snap)
422 : goto error;
423 :
424 16754370 : level = 0;
425 16754370 : init = 1;
426 16754370 : reduced = 0;
427 :
428 44466976 : while (level >= 0) {
429 26639311 : if (init) {
430 : int choice;
431 :
432 24380362 : res = compute_min(ctx, tab, min, level);
433 24380362 : if (res == isl_lp_error)
434 0 : goto error;
435 24380362 : if (res != isl_lp_ok)
436 0 : isl_die(ctx, isl_error_internal,
437 : "expecting bounded rational solution",
438 : goto error);
439 24380362 : if (isl_tab_sample_is_integer(tab))
440 2834584 : break;
441 21545778 : res = compute_max(ctx, tab, max, level);
442 21545778 : if (res == isl_lp_error)
443 0 : goto error;
444 21545778 : if (res != isl_lp_ok)
445 0 : isl_die(ctx, isl_error_internal,
446 : "expecting bounded rational solution",
447 : goto error);
448 21545778 : if (isl_tab_sample_is_integer(tab))
449 955911 : break;
450 20589867 : choice = isl_int_lt(min->el[level], max->el[level]);
451 20589867 : if (choice) {
452 : int g;
453 15707993 : g = greedy_search(ctx, tab, min, max, level);
454 15707993 : if (g < 0)
455 0 : goto error;
456 15707993 : if (g)
457 11841111 : break;
458 : }
459 12050846 : if (!reduced && choice &&
460 3302090 : ctx->opt->gbr != ISL_GBR_NEVER) {
461 : unsigned gbr_only_first;
462 3302090 : if (ctx->opt->gbr == ISL_GBR_ONCE)
463 0 : ctx->opt->gbr = ISL_GBR_NEVER;
464 3302090 : tab->n_zero = level;
465 3302090 : gbr_only_first = ctx->opt->gbr_only_first;
466 6604180 : ctx->opt->gbr_only_first =
467 3302090 : ctx->opt->gbr == ISL_GBR_ALWAYS;
468 3302090 : tab = isl_tab_compute_reduced_basis(tab);
469 3302090 : ctx->opt->gbr_only_first = gbr_only_first;
470 3302090 : if (!tab || !tab->basis)
471 : goto error;
472 3302090 : reduced = 1;
473 3302090 : continue;
474 : }
475 5446666 : reduced = 0;
476 5446666 : snap[level] = isl_tab_snap(tab);
477 : } else
478 2258949 : isl_int_add_ui(min->el[level], min->el[level], 1);
479 :
480 7705615 : if (isl_int_gt(min->el[level], max->el[level])) {
481 3332244 : level--;
482 3332244 : init = 0;
483 3332244 : if (level >= 0)
484 2258949 : if (isl_tab_rollback(tab, snap[level]) < 0)
485 0 : goto error;
486 3332244 : continue;
487 : }
488 4373371 : isl_int_neg(tab->basis->row[1 + level][0], min->el[level]);
489 4373371 : if (isl_tab_add_valid_eq(tab, tab->basis->row[1 + level]) < 0)
490 0 : goto error;
491 4373371 : isl_int_set_si(tab->basis->row[1 + level][0], 0);
492 4373371 : if (level + tab->n_unbounded < dim - 1) {
493 4323902 : ++level;
494 4323902 : init = 1;
495 4323902 : continue;
496 : }
497 49469 : break;
498 : }
499 :
500 16754370 : if (level >= 0) {
501 15681075 : sample = isl_tab_get_sample_value(tab);
502 15681075 : if (!sample)
503 0 : goto error;
504 15681075 : 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 1073295 : sample = isl_vec_alloc(ctx, 0);
513 :
514 16754370 : ctx->opt->gbr = gbr;
515 16754370 : isl_vec_free(min);
516 16754370 : isl_vec_free(max);
517 16754370 : free(snap);
518 16754370 : 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 978268 : 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 978268 : isl_vec *sample = NULL;
537 : isl_ctx *ctx;
538 : unsigned nparam;
539 : unsigned nvar;
540 :
541 978268 : ctx = isl_basic_set_get_ctx(bset);
542 978268 : if (!ctx)
543 0 : goto error;
544 :
545 978268 : nparam = isl_basic_set_dim(bset, isl_dim_param);
546 978268 : nvar = isl_basic_set_dim(bset, isl_dim_set);
547 :
548 978268 : sample = isl_vec_alloc(ctx, 1 + isl_basic_set_total_dim(bset));
549 978268 : if (!sample)
550 0 : goto error;
551 978268 : isl_int_set_si(sample->el[0], 1);
552 :
553 978268 : bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
554 :
555 1183414 : for (i = 0, n = 0; i < f->n_group; ++i) {
556 : isl_basic_set *bset_i;
557 : isl_vec *sample_i;
558 :
559 1092760 : bset_i = isl_basic_set_copy(bset);
560 2185520 : bset_i = isl_basic_set_drop_constraints_involving(bset_i,
561 2185520 : nparam + n + f->len[i], nvar - n - f->len[i]);
562 1092760 : bset_i = isl_basic_set_drop_constraints_involving(bset_i,
563 : nparam, n);
564 2185520 : bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
565 2185520 : n + f->len[i], nvar - n - f->len[i]);
566 1092760 : bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
567 :
568 1092760 : sample_i = sample_bounded(bset_i);
569 1092760 : if (!sample_i)
570 0 : goto error;
571 1092760 : if (sample_i->size == 0) {
572 887614 : isl_basic_set_free(bset);
573 887614 : isl_factorizer_free(f);
574 887614 : isl_vec_free(sample);
575 887614 : return sample_i;
576 : }
577 410292 : isl_seq_cpy(sample->el + 1 + nparam + n,
578 410292 : sample_i->el + 1, f->len[i]);
579 205146 : isl_vec_free(sample_i);
580 :
581 205146 : n += f->len[i];
582 : }
583 :
584 90654 : f->morph = isl_morph_inverse(f->morph);
585 90654 : sample = isl_morph_vec(isl_morph_copy(f->morph), sample);
586 :
587 90654 : isl_basic_set_free(bset);
588 90654 : isl_factorizer_free(f);
589 90654 : 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 1396119374 : static __isl_give isl_vec *sample_bounded(__isl_take isl_basic_set *bset)
605 : {
606 : unsigned dim;
607 : struct isl_vec *sample;
608 1396119374 : struct isl_tab *tab = NULL;
609 : isl_factorizer *f;
610 :
611 1396119374 : if (!bset)
612 0 : return NULL;
613 :
614 1396119374 : if (isl_basic_set_plain_is_empty(bset))
615 0 : return empty_sample(bset);
616 :
617 1396119374 : dim = isl_basic_set_total_dim(bset);
618 1396119374 : if (dim == 0)
619 56235738 : return zero_sample(bset);
620 1339883636 : if (dim == 1)
621 985340 : return interval_sample(bset);
622 1338898296 : if (bset->n_eq > 0)
623 0 : return sample_eq(bset, sample_bounded);
624 :
625 1338898296 : f = isl_basic_set_factorizer(bset);
626 1338898296 : if (!f)
627 0 : goto error;
628 1338898296 : if (f->n_group != 0)
629 978268 : return factored_sample(bset, f);
630 1337920028 : isl_factorizer_free(f);
631 :
632 1337920028 : tab = isl_tab_from_basic_set(bset, 1);
633 1337920028 : if (tab && tab->empty) {
634 1321233234 : isl_tab_free(tab);
635 1321233234 : ISL_F_SET(bset, ISL_BASIC_SET_EMPTY);
636 1321233234 : sample = isl_vec_alloc(isl_basic_set_get_ctx(bset), 0);
637 1321233234 : isl_basic_set_free(bset);
638 1321233234 : return sample;
639 : }
640 :
641 16686794 : if (!ISL_F_ISSET(bset, ISL_BASIC_SET_NO_IMPLICIT))
642 16643366 : if (isl_tab_detect_implicit_equalities(tab) < 0)
643 0 : goto error;
644 :
645 16686794 : sample = isl_tab_sample(tab);
646 16686794 : if (!sample)
647 0 : goto error;
648 :
649 16686794 : if (sample->size > 0) {
650 15497872 : isl_vec_free(bset->sample);
651 15497872 : bset->sample = isl_vec_copy(sample);
652 : }
653 :
654 16686794 : isl_basic_set_free(bset);
655 16686794 : isl_tab_free(tab);
656 16686794 : 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 60440261 : 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 60440261 : if (!bset || !sample)
683 : goto error;
684 :
685 60440261 : total = isl_basic_set_total_dim(bset);
686 60440261 : T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1));
687 60440261 : if (!T)
688 0 : goto error;
689 :
690 132770719 : for (i = 0; i < sample->size; ++i) {
691 72330458 : isl_int_set(T->row[i][0], sample->el[i]);
692 72330458 : isl_seq_clr(T->row[i] + 1, T->n_col - 1);
693 : }
694 360689842 : for (i = 0; i < T->n_col - 1; ++i) {
695 300249581 : isl_seq_clr(T->row[sample->size + i], T->n_col);
696 300249581 : isl_int_set_si(T->row[sample->size + i][1 + i], 1);
697 : }
698 60440261 : isl_vec_free(sample);
699 :
700 60440261 : bset = isl_basic_set_preimage(bset, T);
701 60440261 : 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 81092712 : 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 81092712 : if (!bset)
717 0 : return NULL;
718 :
719 81092712 : tab = isl_tab_from_basic_set(bset, 0);
720 81092712 : sample = isl_tab_get_sample_value(tab);
721 81092712 : isl_tab_free(tab);
722 :
723 81092712 : isl_basic_set_free(bset);
724 :
725 81092712 : 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 20652451 : 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 20652451 : struct isl_basic_set *shift = NULL;
768 :
769 20652451 : if (!cone || !vec)
770 : goto error;
771 :
772 20652451 : isl_assert(cone->ctx, cone->n_eq == 0, goto error);
773 :
774 20652451 : total = isl_basic_set_total_dim(cone);
775 :
776 20652451 : shift = isl_basic_set_alloc_space(isl_basic_set_get_space(cone),
777 : 0, 0, cone->n_ineq);
778 :
779 134141551 : for (i = 0; i < cone->n_ineq; ++i) {
780 113489100 : k = isl_basic_set_alloc_inequality(shift);
781 113489100 : if (k < 0)
782 0 : goto error;
783 113489100 : isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total);
784 113489100 : isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total,
785 113489100 : &shift->ineq[k][0]);
786 113489100 : isl_int_cdiv_q(shift->ineq[k][0],
787 : shift->ineq[k][0], vec->el[0]);
788 113489100 : isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]);
789 723883911 : for (j = 0; j < total; ++j) {
790 610394811 : if (isl_int_is_nonneg(shift->ineq[k][1 + j]))
791 459929358 : continue;
792 150465453 : isl_int_add(shift->ineq[k][0],
793 : shift->ineq[k][0], shift->ineq[k][1 + j]);
794 : }
795 : }
796 :
797 20652451 : isl_basic_set_free(cone);
798 20652451 : isl_vec_free(vec);
799 :
800 20652451 : 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 60440261 : 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 60440261 : if (!vec || !cone || !U)
825 : goto error;
826 :
827 60440261 : isl_assert(vec->ctx, vec->size != 0, goto error);
828 60440261 : if (isl_int_is_one(vec->el[0])) {
829 39787810 : isl_mat_free(U);
830 39787810 : isl_basic_set_free(cone);
831 39787810 : return vec;
832 : }
833 :
834 20652451 : total = isl_basic_set_total_dim(cone);
835 20652451 : cone = isl_basic_set_preimage(cone, U);
836 20652451 : cone = isl_basic_set_remove_dims(cone, isl_dim_set,
837 20652451 : 0, total - (vec->size - 1));
838 :
839 20652451 : cone = shift_cone(cone, vec);
840 :
841 20652451 : vec = rational_sample(cone);
842 20652451 : vec = isl_vec_ceil(vec);
843 20652451 : 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 60440261 : 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 60440261 : if (!vec1 || !vec2)
860 : goto error;
861 60440261 : isl_assert(vec1->ctx, vec1->size > 0, goto error);
862 60440261 : isl_assert(vec2->ctx, vec2->size > 0, goto error);
863 60440261 : isl_assert(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error);
864 60440261 : isl_assert(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error);
865 :
866 60440261 : vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1);
867 60440261 : if (!vec)
868 0 : goto error;
869 :
870 60440261 : isl_seq_cpy(vec->el, vec1->el, vec1->size);
871 60440261 : isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1);
872 :
873 60440261 : isl_vec_free(vec1);
874 60440261 : isl_vec_free(vec2);
875 :
876 60440261 : 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 787316449 : __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 787316449 : if (!bset || !cone)
932 : goto error;
933 :
934 787316449 : ctx = isl_basic_set_get_ctx(bset);
935 787316449 : total = isl_basic_set_total_dim(cone);
936 787316449 : cone_dim = total - cone->n_eq;
937 :
938 787316449 : M = isl_mat_sub_alloc6(ctx, cone->eq, 0, cone->n_eq, 1, total);
939 787316449 : M = isl_mat_left_hermite(M, 0, &U, NULL);
940 787316449 : if (!M)
941 0 : goto error;
942 787316449 : isl_mat_free(M);
943 :
944 787316449 : U = isl_mat_lin_to_aff(U);
945 787316449 : bset = isl_basic_set_preimage(bset, isl_mat_copy(U));
946 :
947 787316449 : bounded = isl_basic_set_copy(bset);
948 787316449 : bounded = isl_basic_set_drop_constraints_involving(bounded,
949 : total - cone_dim, cone_dim);
950 787316449 : bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim);
951 787316449 : sample = sample_bounded(bounded);
952 787316449 : if (!sample || sample->size == 0) {
953 726876188 : isl_basic_set_free(bset);
954 726876188 : isl_basic_set_free(cone);
955 726876188 : isl_mat_free(U);
956 726876188 : return sample;
957 : }
958 60440261 : bset = plug_in(bset, isl_vec_copy(sample));
959 60440261 : cone_sample = rational_sample(bset);
960 60440261 : cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(U));
961 60440261 : sample = vec_concat(sample, cone_sample);
962 60440261 : sample = isl_mat_vec_product(U, sample);
963 60440261 : 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 1395026614 : 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 1395026614 : dim = isl_basic_set_total_dim(bset);
1110 :
1111 1395026614 : cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset));
1112 1395026614 : if (!cone)
1113 0 : goto error;
1114 :
1115 1395026614 : if (cone->n_eq < dim)
1116 787316449 : return isl_basic_set_sample_with_cone(bset, cone);
1117 :
1118 607710165 : isl_basic_set_free(cone);
1119 607710165 : return sample_bounded(bset);
1120 : error:
1121 0 : isl_basic_set_free(bset);
1122 0 : return NULL;
1123 : }
1124 :
1125 1716161429 : 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 1716161429 : if (!bset)
1131 0 : return NULL;
1132 :
1133 1716161429 : ctx = bset->ctx;
1134 1716161429 : if (isl_basic_set_plain_is_empty(bset))
1135 497843 : return empty_sample(bset);
1136 :
1137 1715663586 : dim = isl_basic_set_n_dim(bset);
1138 1715663586 : isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
1139 1715663586 : isl_assert(ctx, bset->n_div == 0, goto error);
1140 :
1141 1715663586 : 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 1715663586 : isl_vec_free(bset->sample);
1152 1715663586 : bset->sample = NULL;
1153 :
1154 1715663586 : if (bset->n_eq > 0)
1155 318523352 : return sample_eq(bset, bounded ? isl_basic_set_sample_bounded
1156 : : isl_basic_set_sample_vec);
1157 1397140234 : if (dim == 0)
1158 48512 : return zero_sample(bset);
1159 1397091722 : if (dim == 1)
1160 2065108 : return interval_sample(bset);
1161 :
1162 1395026614 : return bounded ? sample_bounded(bset) : gbr_sample(bset);
1163 : error:
1164 0 : isl_basic_set_free(bset);
1165 0 : return NULL;
1166 : }
1167 :
1168 1716161429 : __isl_give isl_vec *isl_basic_set_sample_vec(__isl_take isl_basic_set *bset)
1169 : {
1170 1716161429 : 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 2090921 : __isl_give isl_basic_set *isl_basic_set_from_vec(__isl_take isl_vec *vec)
1182 : {
1183 : int i;
1184 : int k;
1185 2090921 : struct isl_basic_set *bset = NULL;
1186 : struct isl_ctx *ctx;
1187 : unsigned dim;
1188 :
1189 2090921 : if (!vec)
1190 0 : return NULL;
1191 2090921 : ctx = vec->ctx;
1192 2090921 : isl_assert(ctx, vec->size != 0, goto error);
1193 :
1194 2090921 : bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
1195 2090921 : if (!bset)
1196 0 : goto error;
1197 2090921 : dim = isl_basic_set_n_dim(bset);
1198 13227509 : for (i = dim - 1; i >= 0; --i) {
1199 11136588 : k = isl_basic_set_alloc_equality(bset);
1200 11136588 : if (k < 0)
1201 0 : goto error;
1202 11136588 : isl_seq_clr(bset->eq[k], 1 + dim);
1203 11136588 : isl_int_neg(bset->eq[k][0], vec->el[1 + i]);
1204 11136588 : isl_int_set(bset->eq[k][1 + i], vec->el[0]);
1205 : }
1206 2090921 : bset->sample = vec;
1207 :
1208 2090921 : 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 : }
|
