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