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/lp.h> 13#include <isl/map.h> 14#include <isl_mat_private.h> 15#include <isl/set.h> 16#include <isl/seq.h> 17#include <isl_options_private.h> 18#include "isl_equalities.h" 19#include "isl_tab.h" 20 21static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set); 22 23/* Return 1 if constraint c is redundant with respect to the constraints 24 * in bmap. If c is a lower [upper] bound in some variable and bmap 25 * does not have a lower [upper] bound in that variable, then c cannot 26 * be redundant and we do not need solve any lp. 27 */ 28int isl_basic_map_constraint_is_redundant(struct isl_basic_map **bmap, 29 isl_int *c, isl_int *opt_n, isl_int *opt_d) 30{ 31 enum isl_lp_result res; 32 unsigned total; 33 int i, j; 34 35 if (!bmap) 36 return -1; 37 38 total = isl_basic_map_total_dim(*bmap); 39 for (i = 0; i < total; ++i) { 40 int sign; 41 if (isl_int_is_zero(c[1+i])) 42 continue; 43 sign = isl_int_sgn(c[1+i]); 44 for (j = 0; j < (*bmap)->n_ineq; ++j) 45 if (sign == isl_int_sgn((*bmap)->ineq[j][1+i])) 46 break; 47 if (j == (*bmap)->n_ineq) 48 break; 49 } 50 if (i < total) 51 return 0; 52 53 res = isl_basic_map_solve_lp(*bmap, 0, c, (*bmap)->ctx->one, 54 opt_n, opt_d, NULL); 55 if (res == isl_lp_unbounded) 56 return 0; 57 if (res == isl_lp_error) 58 return -1; 59 if (res == isl_lp_empty) { 60 *bmap = isl_basic_map_set_to_empty(*bmap); 61 return 0; 62 } 63 return !isl_int_is_neg(*opt_n); 64} 65 66int isl_basic_set_constraint_is_redundant(struct isl_basic_set **bset, 67 isl_int *c, isl_int *opt_n, isl_int *opt_d) 68{ 69 return isl_basic_map_constraint_is_redundant( 70 (struct isl_basic_map **)bset, c, opt_n, opt_d); 71} 72 73/* Remove redundant 74 * constraints. If the minimal value along the normal of a constraint 75 * is the same if the constraint is removed, then the constraint is redundant. 76 * 77 * Alternatively, we could have intersected the basic map with the 78 * corresponding equality and the checked if the dimension was that 79 * of a facet. 80 */ 81__isl_give isl_basic_map *isl_basic_map_remove_redundancies( 82 __isl_take isl_basic_map *bmap) 83{ 84 struct isl_tab *tab; 85 86 if (!bmap) 87 return NULL; 88 89 bmap = isl_basic_map_gauss(bmap, NULL); 90 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) 91 return bmap; 92 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NO_REDUNDANT)) 93 return bmap; 94 if (bmap->n_ineq <= 1) 95 return bmap; 96 97 tab = isl_tab_from_basic_map(bmap, 0); 98 if (isl_tab_detect_implicit_equalities(tab) < 0) 99 goto error; 100 if (isl_tab_detect_redundant(tab) < 0) 101 goto error; 102 bmap = isl_basic_map_update_from_tab(bmap, tab); 103 isl_tab_free(tab); 104 ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT); 105 ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT); 106 return bmap; 107error: 108 isl_tab_free(tab); 109 isl_basic_map_free(bmap); 110 return NULL; 111} 112 113__isl_give isl_basic_set *isl_basic_set_remove_redundancies( 114 __isl_take isl_basic_set *bset) 115{ 116 return (struct isl_basic_set *) 117 isl_basic_map_remove_redundancies((struct isl_basic_map *)bset); 118} 119 120/* Remove redundant constraints in each of the basic maps. 121 */ 122__isl_give isl_map *isl_map_remove_redundancies(__isl_take isl_map *map) 123{ 124 return isl_map_inline_foreach_basic_map(map, 125 &isl_basic_map_remove_redundancies); 126} 127 128__isl_give isl_set *isl_set_remove_redundancies(__isl_take isl_set *set) 129{ 130 return isl_map_remove_redundancies(set); 131} 132 133/* Check if the set set is bound in the direction of the affine 134 * constraint c and if so, set the constant term such that the 135 * resulting constraint is a bounding constraint for the set. 136 */ 137static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len) 138{ 139 int first; 140 int j; 141 isl_int opt; 142 isl_int opt_denom; 143 144 isl_int_init(opt); 145 isl_int_init(opt_denom); 146 first = 1; 147 for (j = 0; j < set->n; ++j) { 148 enum isl_lp_result res; 149 150 if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY)) 151 continue; 152 153 res = isl_basic_set_solve_lp(set->p[j], 154 0, c, set->ctx->one, &opt, &opt_denom, NULL); 155 if (res == isl_lp_unbounded) 156 break; 157 if (res == isl_lp_error) 158 goto error; 159 if (res == isl_lp_empty) { 160 set->p[j] = isl_basic_set_set_to_empty(set->p[j]); 161 if (!set->p[j]) 162 goto error; 163 continue; 164 } 165 if (first || isl_int_is_neg(opt)) { 166 if (!isl_int_is_one(opt_denom)) 167 isl_seq_scale(c, c, opt_denom, len); 168 isl_int_sub(c[0], c[0], opt); 169 } 170 first = 0; 171 } 172 isl_int_clear(opt); 173 isl_int_clear(opt_denom); 174 return j >= set->n; 175error: 176 isl_int_clear(opt); 177 isl_int_clear(opt_denom); 178 return -1; 179} 180 181__isl_give isl_basic_map *isl_basic_map_set_rational( 182 __isl_take isl_basic_set *bmap) 183{ 184 if (!bmap) 185 return NULL; 186 187 if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) 188 return bmap; 189 190 bmap = isl_basic_map_cow(bmap); 191 if (!bmap) 192 return NULL; 193 194 ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL); 195 196 return isl_basic_map_finalize(bmap); 197} 198 199__isl_give isl_basic_set *isl_basic_set_set_rational( 200 __isl_take isl_basic_set *bset) 201{ 202 return isl_basic_map_set_rational(bset); 203} 204 205__isl_give isl_map *isl_map_set_rational(__isl_take isl_map *map) 206{ 207 int i; 208 209 map = isl_map_cow(map); 210 if (!map) 211 return NULL; 212 for (i = 0; i < map->n; ++i) { 213 map->p[i] = isl_basic_map_set_rational(map->p[i]); 214 if (!map->p[i]) 215 goto error; 216 } 217 return map; 218error: 219 isl_map_free(map); 220 return NULL; 221} 222 223__isl_give isl_set *isl_set_set_rational(__isl_take isl_set *set) 224{ 225 return isl_map_set_rational(set); 226} 227 228static struct isl_basic_set *isl_basic_set_add_equality( 229 struct isl_basic_set *bset, isl_int *c) 230{ 231 int i; 232 unsigned dim; 233 234 if (!bset) 235 return NULL; 236 237 if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY)) 238 return bset; 239 240 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 241 isl_assert(bset->ctx, bset->n_div == 0, goto error); 242 dim = isl_basic_set_n_dim(bset); 243 bset = isl_basic_set_cow(bset); 244 bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0); 245 i = isl_basic_set_alloc_equality(bset); 246 if (i < 0) 247 goto error; 248 isl_seq_cpy(bset->eq[i], c, 1 + dim); 249 return bset; 250error: 251 isl_basic_set_free(bset); 252 return NULL; 253} 254 255static struct isl_set *isl_set_add_basic_set_equality(struct isl_set *set, isl_int *c) 256{ 257 int i; 258 259 set = isl_set_cow(set); 260 if (!set) 261 return NULL; 262 for (i = 0; i < set->n; ++i) { 263 set->p[i] = isl_basic_set_add_equality(set->p[i], c); 264 if (!set->p[i]) 265 goto error; 266 } 267 return set; 268error: 269 isl_set_free(set); 270 return NULL; 271} 272 273/* Given a union of basic sets, construct the constraints for wrapping 274 * a facet around one of its ridges. 275 * In particular, if each of n the d-dimensional basic sets i in "set" 276 * contains the origin, satisfies the constraints x_1 >= 0 and x_2 >= 0 277 * and is defined by the constraints 278 * [ 1 ] 279 * A_i [ x ] >= 0 280 * 281 * then the resulting set is of dimension n*(1+d) and has as constraints 282 * 283 * [ a_i ] 284 * A_i [ x_i ] >= 0 285 * 286 * a_i >= 0 287 * 288 * \sum_i x_{i,1} = 1 289 */ 290static struct isl_basic_set *wrap_constraints(struct isl_set *set) 291{ 292 struct isl_basic_set *lp; 293 unsigned n_eq; 294 unsigned n_ineq; 295 int i, j, k; 296 unsigned dim, lp_dim; 297 298 if (!set) 299 return NULL; 300 301 dim = 1 + isl_set_n_dim(set); 302 n_eq = 1; 303 n_ineq = set->n; 304 for (i = 0; i < set->n; ++i) { 305 n_eq += set->p[i]->n_eq; 306 n_ineq += set->p[i]->n_ineq; 307 } 308 lp = isl_basic_set_alloc(set->ctx, 0, dim * set->n, 0, n_eq, n_ineq); 309 lp = isl_basic_set_set_rational(lp); 310 if (!lp) 311 return NULL; 312 lp_dim = isl_basic_set_n_dim(lp); 313 k = isl_basic_set_alloc_equality(lp); 314 isl_int_set_si(lp->eq[k][0], -1); 315 for (i = 0; i < set->n; ++i) { 316 isl_int_set_si(lp->eq[k][1+dim*i], 0); 317 isl_int_set_si(lp->eq[k][1+dim*i+1], 1); 318 isl_seq_clr(lp->eq[k]+1+dim*i+2, dim-2); 319 } 320 for (i = 0; i < set->n; ++i) { 321 k = isl_basic_set_alloc_inequality(lp); 322 isl_seq_clr(lp->ineq[k], 1+lp_dim); 323 isl_int_set_si(lp->ineq[k][1+dim*i], 1); 324 325 for (j = 0; j < set->p[i]->n_eq; ++j) { 326 k = isl_basic_set_alloc_equality(lp); 327 isl_seq_clr(lp->eq[k], 1+dim*i); 328 isl_seq_cpy(lp->eq[k]+1+dim*i, set->p[i]->eq[j], dim); 329 isl_seq_clr(lp->eq[k]+1+dim*(i+1), dim*(set->n-i-1)); 330 } 331 332 for (j = 0; j < set->p[i]->n_ineq; ++j) { 333 k = isl_basic_set_alloc_inequality(lp); 334 isl_seq_clr(lp->ineq[k], 1+dim*i); 335 isl_seq_cpy(lp->ineq[k]+1+dim*i, set->p[i]->ineq[j], dim); 336 isl_seq_clr(lp->ineq[k]+1+dim*(i+1), dim*(set->n-i-1)); 337 } 338 } 339 return lp; 340} 341 342/* Given a facet "facet" of the convex hull of "set" and a facet "ridge" 343 * of that facet, compute the other facet of the convex hull that contains 344 * the ridge. 345 * 346 * We first transform the set such that the facet constraint becomes 347 * 348 * x_1 >= 0 349 * 350 * I.e., the facet lies in 351 * 352 * x_1 = 0 353 * 354 * and on that facet, the constraint that defines the ridge is 355 * 356 * x_2 >= 0 357 * 358 * (This transformation is not strictly needed, all that is needed is 359 * that the ridge contains the origin.) 360 * 361 * Since the ridge contains the origin, the cone of the convex hull 362 * will be of the form 363 * 364 * x_1 >= 0 365 * x_2 >= a x_1 366 * 367 * with this second constraint defining the new facet. 368 * The constant a is obtained by settting x_1 in the cone of the 369 * convex hull to 1 and minimizing x_2. 370 * Now, each element in the cone of the convex hull is the sum 371 * of elements in the cones of the basic sets. 372 * If a_i is the dilation factor of basic set i, then the problem 373 * we need to solve is 374 * 375 * min \sum_i x_{i,2} 376 * st 377 * \sum_i x_{i,1} = 1 378 * a_i >= 0 379 * [ a_i ] 380 * A [ x_i ] >= 0 381 * 382 * with 383 * [ 1 ] 384 * A_i [ x_i ] >= 0 385 * 386 * the constraints of each (transformed) basic set. 387 * If a = n/d, then the constraint defining the new facet (in the transformed 388 * space) is 389 * 390 * -n x_1 + d x_2 >= 0 391 * 392 * In the original space, we need to take the same combination of the 393 * corresponding constraints "facet" and "ridge". 394 * 395 * If a = -infty = "-1/0", then we just return the original facet constraint. 396 * This means that the facet is unbounded, but has a bounded intersection 397 * with the union of sets. 398 */ 399isl_int *isl_set_wrap_facet(__isl_keep isl_set *set, 400 isl_int *facet, isl_int *ridge) 401{ 402 int i; 403 isl_ctx *ctx; 404 struct isl_mat *T = NULL; 405 struct isl_basic_set *lp = NULL; 406 struct isl_vec *obj; 407 enum isl_lp_result res; 408 isl_int num, den; 409 unsigned dim; 410 411 if (!set) 412 return NULL; 413 ctx = set->ctx; 414 set = isl_set_copy(set); 415 set = isl_set_set_rational(set); 416 417 dim = 1 + isl_set_n_dim(set); 418 T = isl_mat_alloc(ctx, 3, dim); 419 if (!T) 420 goto error; 421 isl_int_set_si(T->row[0][0], 1); 422 isl_seq_clr(T->row[0]+1, dim - 1); 423 isl_seq_cpy(T->row[1], facet, dim); 424 isl_seq_cpy(T->row[2], ridge, dim); 425 T = isl_mat_right_inverse(T); 426 set = isl_set_preimage(set, T); 427 T = NULL; 428 if (!set) 429 goto error; 430 lp = wrap_constraints(set); 431 obj = isl_vec_alloc(ctx, 1 + dim*set->n); 432 if (!obj) 433 goto error; 434 isl_int_set_si(obj->block.data[0], 0); 435 for (i = 0; i < set->n; ++i) { 436 isl_seq_clr(obj->block.data + 1 + dim*i, 2); 437 isl_int_set_si(obj->block.data[1 + dim*i+2], 1); 438 isl_seq_clr(obj->block.data + 1 + dim*i+3, dim-3); 439 } 440 isl_int_init(num); 441 isl_int_init(den); 442 res = isl_basic_set_solve_lp(lp, 0, 443 obj->block.data, ctx->one, &num, &den, NULL); 444 if (res == isl_lp_ok) { 445 isl_int_neg(num, num); 446 isl_seq_combine(facet, num, facet, den, ridge, dim); 447 isl_seq_normalize(ctx, facet, dim); 448 } 449 isl_int_clear(num); 450 isl_int_clear(den); 451 isl_vec_free(obj); 452 isl_basic_set_free(lp); 453 isl_set_free(set); 454 if (res == isl_lp_error) 455 return NULL; 456 isl_assert(ctx, res == isl_lp_ok || res == isl_lp_unbounded, 457 return NULL); 458 return facet; 459error: 460 isl_basic_set_free(lp); 461 isl_mat_free(T); 462 isl_set_free(set); 463 return NULL; 464} 465 466/* Compute the constraint of a facet of "set". 467 * 468 * We first compute the intersection with a bounding constraint 469 * that is orthogonal to one of the coordinate axes. 470 * If the affine hull of this intersection has only one equality, 471 * we have found a facet. 472 * Otherwise, we wrap the current bounding constraint around 473 * one of the equalities of the face (one that is not equal to 474 * the current bounding constraint). 475 * This process continues until we have found a facet. 476 * The dimension of the intersection increases by at least 477 * one on each iteration, so termination is guaranteed. 478 */ 479static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set) 480{ 481 struct isl_set *slice = NULL; 482 struct isl_basic_set *face = NULL; 483 int i; 484 unsigned dim = isl_set_n_dim(set); 485 int is_bound; 486 isl_mat *bounds = NULL; 487 488 isl_assert(set->ctx, set->n > 0, goto error); 489 bounds = isl_mat_alloc(set->ctx, 1, 1 + dim); 490 if (!bounds) 491 return NULL; 492 493 isl_seq_clr(bounds->row[0], dim); 494 isl_int_set_si(bounds->row[0][1 + dim - 1], 1); 495 is_bound = uset_is_bound(set, bounds->row[0], 1 + dim); 496 if (is_bound < 0) 497 goto error; 498 isl_assert(set->ctx, is_bound, goto error); 499 isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim); 500 bounds->n_row = 1; 501 502 for (;;) { 503 slice = isl_set_copy(set); 504 slice = isl_set_add_basic_set_equality(slice, bounds->row[0]); 505 face = isl_set_affine_hull(slice); 506 if (!face) 507 goto error; 508 if (face->n_eq == 1) { 509 isl_basic_set_free(face); 510 break; 511 } 512 for (i = 0; i < face->n_eq; ++i) 513 if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) && 514 !isl_seq_is_neg(bounds->row[0], 515 face->eq[i], 1 + dim)) 516 break; 517 isl_assert(set->ctx, i < face->n_eq, goto error); 518 if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i])) 519 goto error; 520 isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col); 521 isl_basic_set_free(face); 522 } 523 524 return bounds; 525error: 526 isl_basic_set_free(face); 527 isl_mat_free(bounds); 528 return NULL; 529} 530 531/* Given the bounding constraint "c" of a facet of the convex hull of "set", 532 * compute a hyperplane description of the facet, i.e., compute the facets 533 * of the facet. 534 * 535 * We compute an affine transformation that transforms the constraint 536 * 537 * [ 1 ] 538 * c [ x ] = 0 539 * 540 * to the constraint 541 * 542 * z_1 = 0 543 * 544 * by computing the right inverse U of a matrix that starts with the rows 545 * 546 * [ 1 0 ] 547 * [ c ] 548 * 549 * Then 550 * [ 1 ] [ 1 ] 551 * [ x ] = U [ z ] 552 * and 553 * [ 1 ] [ 1 ] 554 * [ z ] = Q [ x ] 555 * 556 * with Q = U^{-1} 557 * Since z_1 is zero, we can drop this variable as well as the corresponding 558 * column of U to obtain 559 * 560 * [ 1 ] [ 1 ] 561 * [ x ] = U' [ z' ] 562 * and 563 * [ 1 ] [ 1 ] 564 * [ z' ] = Q' [ x ] 565 * 566 * with Q' equal to Q, but without the corresponding row. 567 * After computing the facets of the facet in the z' space, 568 * we convert them back to the x space through Q. 569 */ 570static struct isl_basic_set *compute_facet(struct isl_set *set, isl_int *c) 571{ 572 struct isl_mat *m, *U, *Q; 573 struct isl_basic_set *facet = NULL; 574 struct isl_ctx *ctx; 575 unsigned dim; 576 577 ctx = set->ctx; 578 set = isl_set_copy(set); 579 dim = isl_set_n_dim(set); 580 m = isl_mat_alloc(set->ctx, 2, 1 + dim); 581 if (!m) 582 goto error; 583 isl_int_set_si(m->row[0][0], 1); 584 isl_seq_clr(m->row[0]+1, dim); 585 isl_seq_cpy(m->row[1], c, 1+dim); 586 U = isl_mat_right_inverse(m); 587 Q = isl_mat_right_inverse(isl_mat_copy(U)); 588 U = isl_mat_drop_cols(U, 1, 1); 589 Q = isl_mat_drop_rows(Q, 1, 1); 590 set = isl_set_preimage(set, U); 591 facet = uset_convex_hull_wrap_bounded(set); 592 facet = isl_basic_set_preimage(facet, Q); 593 if (facet) 594 isl_assert(ctx, facet->n_eq == 0, goto error); 595 return facet; 596error: 597 isl_basic_set_free(facet); 598 isl_set_free(set); 599 return NULL; 600} 601 602/* Given an initial facet constraint, compute the remaining facets. 603 * We do this by running through all facets found so far and computing 604 * the adjacent facets through wrapping, adding those facets that we 605 * hadn't already found before. 606 * 607 * For each facet we have found so far, we first compute its facets 608 * in the resulting convex hull. That is, we compute the ridges 609 * of the resulting convex hull contained in the facet. 610 * We also compute the corresponding facet in the current approximation 611 * of the convex hull. There is no need to wrap around the ridges 612 * in this facet since that would result in a facet that is already 613 * present in the current approximation. 614 * 615 * This function can still be significantly optimized by checking which of 616 * the facets of the basic sets are also facets of the convex hull and 617 * using all the facets so far to help in constructing the facets of the 618 * facets 619 * and/or 620 * using the technique in section "3.1 Ridge Generation" of 621 * "Extended Convex Hull" by Fukuda et al. 622 */ 623static struct isl_basic_set *extend(struct isl_basic_set *hull, 624 struct isl_set *set) 625{ 626 int i, j, f; 627 int k; 628 struct isl_basic_set *facet = NULL; 629 struct isl_basic_set *hull_facet = NULL; 630 unsigned dim; 631 632 if (!hull) 633 return NULL; 634 635 isl_assert(set->ctx, set->n > 0, goto error); 636 637 dim = isl_set_n_dim(set); 638 639 for (i = 0; i < hull->n_ineq; ++i) { 640 facet = compute_facet(set, hull->ineq[i]); 641 facet = isl_basic_set_add_equality(facet, hull->ineq[i]); 642 facet = isl_basic_set_gauss(facet, NULL); 643 facet = isl_basic_set_normalize_constraints(facet); 644 hull_facet = isl_basic_set_copy(hull); 645 hull_facet = isl_basic_set_add_equality(hull_facet, hull->ineq[i]); 646 hull_facet = isl_basic_set_gauss(hull_facet, NULL); 647 hull_facet = isl_basic_set_normalize_constraints(hull_facet); 648 if (!facet || !hull_facet) 649 goto error; 650 hull = isl_basic_set_cow(hull); 651 hull = isl_basic_set_extend_space(hull, 652 isl_space_copy(hull->dim), 0, 0, facet->n_ineq); 653 if (!hull) 654 goto error; 655 for (j = 0; j < facet->n_ineq; ++j) { 656 for (f = 0; f < hull_facet->n_ineq; ++f) 657 if (isl_seq_eq(facet->ineq[j], 658 hull_facet->ineq[f], 1 + dim)) 659 break; 660 if (f < hull_facet->n_ineq) 661 continue; 662 k = isl_basic_set_alloc_inequality(hull); 663 if (k < 0) 664 goto error; 665 isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim); 666 if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j])) 667 goto error; 668 } 669 isl_basic_set_free(hull_facet); 670 isl_basic_set_free(facet); 671 } 672 hull = isl_basic_set_simplify(hull); 673 hull = isl_basic_set_finalize(hull); 674 return hull; 675error: 676 isl_basic_set_free(hull_facet); 677 isl_basic_set_free(facet); 678 isl_basic_set_free(hull); 679 return NULL; 680} 681 682/* Special case for computing the convex hull of a one dimensional set. 683 * We simply collect the lower and upper bounds of each basic set 684 * and the biggest of those. 685 */ 686static struct isl_basic_set *convex_hull_1d(struct isl_set *set) 687{ 688 struct isl_mat *c = NULL; 689 isl_int *lower = NULL; 690 isl_int *upper = NULL; 691 int i, j, k; 692 isl_int a, b; 693 struct isl_basic_set *hull; 694 695 for (i = 0; i < set->n; ++i) { 696 set->p[i] = isl_basic_set_simplify(set->p[i]); 697 if (!set->p[i]) 698 goto error; 699 } 700 set = isl_set_remove_empty_parts(set); 701 if (!set) 702 goto error; 703 isl_assert(set->ctx, set->n > 0, goto error); 704 c = isl_mat_alloc(set->ctx, 2, 2); 705 if (!c) 706 goto error; 707 708 if (set->p[0]->n_eq > 0) { 709 isl_assert(set->ctx, set->p[0]->n_eq == 1, goto error); 710 lower = c->row[0]; 711 upper = c->row[1]; 712 if (isl_int_is_pos(set->p[0]->eq[0][1])) { 713 isl_seq_cpy(lower, set->p[0]->eq[0], 2); 714 isl_seq_neg(upper, set->p[0]->eq[0], 2); 715 } else { 716 isl_seq_neg(lower, set->p[0]->eq[0], 2); 717 isl_seq_cpy(upper, set->p[0]->eq[0], 2); 718 } 719 } else { 720 for (j = 0; j < set->p[0]->n_ineq; ++j) { 721 if (isl_int_is_pos(set->p[0]->ineq[j][1])) { 722 lower = c->row[0]; 723 isl_seq_cpy(lower, set->p[0]->ineq[j], 2); 724 } else { 725 upper = c->row[1]; 726 isl_seq_cpy(upper, set->p[0]->ineq[j], 2); 727 } 728 } 729 } 730 731 isl_int_init(a); 732 isl_int_init(b); 733 for (i = 0; i < set->n; ++i) { 734 struct isl_basic_set *bset = set->p[i]; 735 int has_lower = 0; 736 int has_upper = 0; 737 738 for (j = 0; j < bset->n_eq; ++j) { 739 has_lower = 1; 740 has_upper = 1; 741 if (lower) { 742 isl_int_mul(a, lower[0], bset->eq[j][1]); 743 isl_int_mul(b, lower[1], bset->eq[j][0]); 744 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1])) 745 isl_seq_cpy(lower, bset->eq[j], 2); 746 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1])) 747 isl_seq_neg(lower, bset->eq[j], 2); 748 } 749 if (upper) { 750 isl_int_mul(a, upper[0], bset->eq[j][1]); 751 isl_int_mul(b, upper[1], bset->eq[j][0]); 752 if (isl_int_lt(a, b) && isl_int_is_pos(bset->eq[j][1])) 753 isl_seq_neg(upper, bset->eq[j], 2); 754 if (isl_int_gt(a, b) && isl_int_is_neg(bset->eq[j][1])) 755 isl_seq_cpy(upper, bset->eq[j], 2); 756 } 757 } 758 for (j = 0; j < bset->n_ineq; ++j) { 759 if (isl_int_is_pos(bset->ineq[j][1])) 760 has_lower = 1; 761 if (isl_int_is_neg(bset->ineq[j][1])) 762 has_upper = 1; 763 if (lower && isl_int_is_pos(bset->ineq[j][1])) { 764 isl_int_mul(a, lower[0], bset->ineq[j][1]); 765 isl_int_mul(b, lower[1], bset->ineq[j][0]); 766 if (isl_int_lt(a, b)) 767 isl_seq_cpy(lower, bset->ineq[j], 2); 768 } 769 if (upper && isl_int_is_neg(bset->ineq[j][1])) { 770 isl_int_mul(a, upper[0], bset->ineq[j][1]); 771 isl_int_mul(b, upper[1], bset->ineq[j][0]); 772 if (isl_int_gt(a, b)) 773 isl_seq_cpy(upper, bset->ineq[j], 2); 774 } 775 } 776 if (!has_lower) 777 lower = NULL; 778 if (!has_upper) 779 upper = NULL; 780 } 781 isl_int_clear(a); 782 isl_int_clear(b); 783 784 hull = isl_basic_set_alloc(set->ctx, 0, 1, 0, 0, 2); 785 hull = isl_basic_set_set_rational(hull); 786 if (!hull) 787 goto error; 788 if (lower) { 789 k = isl_basic_set_alloc_inequality(hull); 790 isl_seq_cpy(hull->ineq[k], lower, 2); 791 } 792 if (upper) { 793 k = isl_basic_set_alloc_inequality(hull); 794 isl_seq_cpy(hull->ineq[k], upper, 2); 795 } 796 hull = isl_basic_set_finalize(hull); 797 isl_set_free(set); 798 isl_mat_free(c); 799 return hull; 800error: 801 isl_set_free(set); 802 isl_mat_free(c); 803 return NULL; 804} 805 806static struct isl_basic_set *convex_hull_0d(struct isl_set *set) 807{ 808 struct isl_basic_set *convex_hull; 809 810 if (!set) 811 return NULL; 812 813 if (isl_set_is_empty(set)) 814 convex_hull = isl_basic_set_empty(isl_space_copy(set->dim)); 815 else 816 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim)); 817 isl_set_free(set); 818 return convex_hull; 819} 820 821/* Compute the convex hull of a pair of basic sets without any parameters or 822 * integer divisions using Fourier-Motzkin elimination. 823 * The convex hull is the set of all points that can be written as 824 * the sum of points from both basic sets (in homogeneous coordinates). 825 * We set up the constraints in a space with dimensions for each of 826 * the three sets and then project out the dimensions corresponding 827 * to the two original basic sets, retaining only those corresponding 828 * to the convex hull. 829 */ 830static struct isl_basic_set *convex_hull_pair_elim(struct isl_basic_set *bset1, 831 struct isl_basic_set *bset2) 832{ 833 int i, j, k; 834 struct isl_basic_set *bset[2]; 835 struct isl_basic_set *hull = NULL; 836 unsigned dim; 837 838 if (!bset1 || !bset2) 839 goto error; 840 841 dim = isl_basic_set_n_dim(bset1); 842 hull = isl_basic_set_alloc(bset1->ctx, 0, 2 + 3 * dim, 0, 843 1 + dim + bset1->n_eq + bset2->n_eq, 844 2 + bset1->n_ineq + bset2->n_ineq); 845 bset[0] = bset1; 846 bset[1] = bset2; 847 for (i = 0; i < 2; ++i) { 848 for (j = 0; j < bset[i]->n_eq; ++j) { 849 k = isl_basic_set_alloc_equality(hull); 850 if (k < 0) 851 goto error; 852 isl_seq_clr(hull->eq[k], (i+1) * (1+dim)); 853 isl_seq_clr(hull->eq[k]+(i+2)*(1+dim), (1-i)*(1+dim)); 854 isl_seq_cpy(hull->eq[k]+(i+1)*(1+dim), bset[i]->eq[j], 855 1+dim); 856 } 857 for (j = 0; j < bset[i]->n_ineq; ++j) { 858 k = isl_basic_set_alloc_inequality(hull); 859 if (k < 0) 860 goto error; 861 isl_seq_clr(hull->ineq[k], (i+1) * (1+dim)); 862 isl_seq_clr(hull->ineq[k]+(i+2)*(1+dim), (1-i)*(1+dim)); 863 isl_seq_cpy(hull->ineq[k]+(i+1)*(1+dim), 864 bset[i]->ineq[j], 1+dim); 865 } 866 k = isl_basic_set_alloc_inequality(hull); 867 if (k < 0) 868 goto error; 869 isl_seq_clr(hull->ineq[k], 1+2+3*dim); 870 isl_int_set_si(hull->ineq[k][(i+1)*(1+dim)], 1); 871 } 872 for (j = 0; j < 1+dim; ++j) { 873 k = isl_basic_set_alloc_equality(hull); 874 if (k < 0) 875 goto error; 876 isl_seq_clr(hull->eq[k], 1+2+3*dim); 877 isl_int_set_si(hull->eq[k][j], -1); 878 isl_int_set_si(hull->eq[k][1+dim+j], 1); 879 isl_int_set_si(hull->eq[k][2*(1+dim)+j], 1); 880 } 881 hull = isl_basic_set_set_rational(hull); 882 hull = isl_basic_set_remove_dims(hull, isl_dim_set, dim, 2*(1+dim)); 883 hull = isl_basic_set_remove_redundancies(hull); 884 isl_basic_set_free(bset1); 885 isl_basic_set_free(bset2); 886 return hull; 887error: 888 isl_basic_set_free(bset1); 889 isl_basic_set_free(bset2); 890 isl_basic_set_free(hull); 891 return NULL; 892} 893 894/* Is the set bounded for each value of the parameters? 895 */ 896int isl_basic_set_is_bounded(__isl_keep isl_basic_set *bset) 897{ 898 struct isl_tab *tab; 899 int bounded; 900 901 if (!bset) 902 return -1; 903 if (isl_basic_set_plain_is_empty(bset)) 904 return 1; 905 906 tab = isl_tab_from_recession_cone(bset, 1); 907 bounded = isl_tab_cone_is_bounded(tab); 908 isl_tab_free(tab); 909 return bounded; 910} 911 912/* Is the image bounded for each value of the parameters and 913 * the domain variables? 914 */ 915int isl_basic_map_image_is_bounded(__isl_keep isl_basic_map *bmap) 916{ 917 unsigned nparam = isl_basic_map_dim(bmap, isl_dim_param); 918 unsigned n_in = isl_basic_map_dim(bmap, isl_dim_in); 919 int bounded; 920 921 bmap = isl_basic_map_copy(bmap); 922 bmap = isl_basic_map_cow(bmap); 923 bmap = isl_basic_map_move_dims(bmap, isl_dim_param, nparam, 924 isl_dim_in, 0, n_in); 925 bounded = isl_basic_set_is_bounded((isl_basic_set *)bmap); 926 isl_basic_map_free(bmap); 927 928 return bounded; 929} 930 931/* Is the set bounded for each value of the parameters? 932 */ 933int isl_set_is_bounded(__isl_keep isl_set *set) 934{ 935 int i; 936 937 if (!set) 938 return -1; 939 940 for (i = 0; i < set->n; ++i) { 941 int bounded = isl_basic_set_is_bounded(set->p[i]); 942 if (!bounded || bounded < 0) 943 return bounded; 944 } 945 return 1; 946} 947 948/* Compute the lineality space of the convex hull of bset1 and bset2. 949 * 950 * We first compute the intersection of the recession cone of bset1 951 * with the negative of the recession cone of bset2 and then compute 952 * the linear hull of the resulting cone. 953 */ 954static struct isl_basic_set *induced_lineality_space( 955 struct isl_basic_set *bset1, struct isl_basic_set *bset2) 956{ 957 int i, k; 958 struct isl_basic_set *lin = NULL; 959 unsigned dim; 960 961 if (!bset1 || !bset2) 962 goto error; 963 964 dim = isl_basic_set_total_dim(bset1); 965 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset1), 0, 966 bset1->n_eq + bset2->n_eq, 967 bset1->n_ineq + bset2->n_ineq); 968 lin = isl_basic_set_set_rational(lin); 969 if (!lin) 970 goto error; 971 for (i = 0; i < bset1->n_eq; ++i) { 972 k = isl_basic_set_alloc_equality(lin); 973 if (k < 0) 974 goto error; 975 isl_int_set_si(lin->eq[k][0], 0); 976 isl_seq_cpy(lin->eq[k] + 1, bset1->eq[i] + 1, dim); 977 } 978 for (i = 0; i < bset1->n_ineq; ++i) { 979 k = isl_basic_set_alloc_inequality(lin); 980 if (k < 0) 981 goto error; 982 isl_int_set_si(lin->ineq[k][0], 0); 983 isl_seq_cpy(lin->ineq[k] + 1, bset1->ineq[i] + 1, dim); 984 } 985 for (i = 0; i < bset2->n_eq; ++i) { 986 k = isl_basic_set_alloc_equality(lin); 987 if (k < 0) 988 goto error; 989 isl_int_set_si(lin->eq[k][0], 0); 990 isl_seq_neg(lin->eq[k] + 1, bset2->eq[i] + 1, dim); 991 } 992 for (i = 0; i < bset2->n_ineq; ++i) { 993 k = isl_basic_set_alloc_inequality(lin); 994 if (k < 0) 995 goto error; 996 isl_int_set_si(lin->ineq[k][0], 0); 997 isl_seq_neg(lin->ineq[k] + 1, bset2->ineq[i] + 1, dim); 998 } 999 1000 isl_basic_set_free(bset1); 1001 isl_basic_set_free(bset2); 1002 return isl_basic_set_affine_hull(lin); 1003error: 1004 isl_basic_set_free(lin); 1005 isl_basic_set_free(bset1); 1006 isl_basic_set_free(bset2); 1007 return NULL; 1008} 1009 1010static struct isl_basic_set *uset_convex_hull(struct isl_set *set); 1011 1012/* Given a set and a linear space "lin" of dimension n > 0, 1013 * project the linear space from the set, compute the convex hull 1014 * and then map the set back to the original space. 1015 * 1016 * Let 1017 * 1018 * M x = 0 1019 * 1020 * describe the linear space. We first compute the Hermite normal 1021 * form H = M U of M = H Q, to obtain 1022 * 1023 * H Q x = 0 1024 * 1025 * The last n rows of H will be zero, so the last n variables of x' = Q x 1026 * are the one we want to project out. We do this by transforming each 1027 * basic set A x >= b to A U x' >= b and then removing the last n dimensions. 1028 * After computing the convex hull in x'_1, i.e., A' x'_1 >= b', 1029 * we transform the hull back to the original space as A' Q_1 x >= b', 1030 * with Q_1 all but the last n rows of Q. 1031 */ 1032static struct isl_basic_set *modulo_lineality(struct isl_set *set, 1033 struct isl_basic_set *lin) 1034{ 1035 unsigned total = isl_basic_set_total_dim(lin); 1036 unsigned lin_dim; 1037 struct isl_basic_set *hull; 1038 struct isl_mat *M, *U, *Q; 1039 1040 if (!set || !lin) 1041 goto error; 1042 lin_dim = total - lin->n_eq; 1043 M = isl_mat_sub_alloc6(set->ctx, lin->eq, 0, lin->n_eq, 1, total); 1044 M = isl_mat_left_hermite(M, 0, &U, &Q); 1045 if (!M) 1046 goto error; 1047 isl_mat_free(M); 1048 isl_basic_set_free(lin); 1049 1050 Q = isl_mat_drop_rows(Q, Q->n_row - lin_dim, lin_dim); 1051 1052 U = isl_mat_lin_to_aff(U); 1053 Q = isl_mat_lin_to_aff(Q); 1054 1055 set = isl_set_preimage(set, U); 1056 set = isl_set_remove_dims(set, isl_dim_set, total - lin_dim, lin_dim); 1057 hull = uset_convex_hull(set); 1058 hull = isl_basic_set_preimage(hull, Q); 1059 1060 return hull; 1061error: 1062 isl_basic_set_free(lin); 1063 isl_set_free(set); 1064 return NULL; 1065} 1066 1067/* Given two polyhedra with as constraints h_{ij} x >= 0 in homegeneous space, 1068 * set up an LP for solving 1069 * 1070 * \sum_j \alpha_{1j} h_{1j} = \sum_j \alpha_{2j} h_{2j} 1071 * 1072 * \alpha{i0} corresponds to the (implicit) positivity constraint 1 >= 0 1073 * The next \alpha{ij} correspond to the equalities and come in pairs. 1074 * The final \alpha{ij} correspond to the inequalities. 1075 */ 1076static struct isl_basic_set *valid_direction_lp( 1077 struct isl_basic_set *bset1, struct isl_basic_set *bset2) 1078{ 1079 isl_space *dim; 1080 struct isl_basic_set *lp; 1081 unsigned d; 1082 int n; 1083 int i, j, k; 1084 1085 if (!bset1 || !bset2) 1086 goto error; 1087 d = 1 + isl_basic_set_total_dim(bset1); 1088 n = 2 + 1089 2 * bset1->n_eq + bset1->n_ineq + 2 * bset2->n_eq + bset2->n_ineq; 1090 dim = isl_space_set_alloc(bset1->ctx, 0, n); 1091 lp = isl_basic_set_alloc_space(dim, 0, d, n); 1092 if (!lp) 1093 goto error; 1094 for (i = 0; i < n; ++i) { 1095 k = isl_basic_set_alloc_inequality(lp); 1096 if (k < 0) 1097 goto error; 1098 isl_seq_clr(lp->ineq[k] + 1, n); 1099 isl_int_set_si(lp->ineq[k][0], -1); 1100 isl_int_set_si(lp->ineq[k][1 + i], 1); 1101 } 1102 for (i = 0; i < d; ++i) { 1103 k = isl_basic_set_alloc_equality(lp); 1104 if (k < 0) 1105 goto error; 1106 n = 0; 1107 isl_int_set_si(lp->eq[k][n], 0); n++; 1108 /* positivity constraint 1 >= 0 */ 1109 isl_int_set_si(lp->eq[k][n], i == 0); n++; 1110 for (j = 0; j < bset1->n_eq; ++j) { 1111 isl_int_set(lp->eq[k][n], bset1->eq[j][i]); n++; 1112 isl_int_neg(lp->eq[k][n], bset1->eq[j][i]); n++; 1113 } 1114 for (j = 0; j < bset1->n_ineq; ++j) { 1115 isl_int_set(lp->eq[k][n], bset1->ineq[j][i]); n++; 1116 } 1117 /* positivity constraint 1 >= 0 */ 1118 isl_int_set_si(lp->eq[k][n], -(i == 0)); n++; 1119 for (j = 0; j < bset2->n_eq; ++j) { 1120 isl_int_neg(lp->eq[k][n], bset2->eq[j][i]); n++; 1121 isl_int_set(lp->eq[k][n], bset2->eq[j][i]); n++; 1122 } 1123 for (j = 0; j < bset2->n_ineq; ++j) { 1124 isl_int_neg(lp->eq[k][n], bset2->ineq[j][i]); n++; 1125 } 1126 } 1127 lp = isl_basic_set_gauss(lp, NULL); 1128 isl_basic_set_free(bset1); 1129 isl_basic_set_free(bset2); 1130 return lp; 1131error: 1132 isl_basic_set_free(bset1); 1133 isl_basic_set_free(bset2); 1134 return NULL; 1135} 1136 1137/* Compute a vector s in the homogeneous space such that <s, r> > 0 1138 * for all rays in the homogeneous space of the two cones that correspond 1139 * to the input polyhedra bset1 and bset2. 1140 * 1141 * We compute s as a vector that satisfies 1142 * 1143 * s = \sum_j \alpha_{ij} h_{ij} for i = 1,2 (*) 1144 * 1145 * with h_{ij} the normals of the facets of polyhedron i 1146 * (including the "positivity constraint" 1 >= 0) and \alpha_{ij} 1147 * strictly positive numbers. For simplicity we impose \alpha_{ij} >= 1. 1148 * We first set up an LP with as variables the \alpha{ij}. 1149 * In this formulation, for each polyhedron i, 1150 * the first constraint is the positivity constraint, followed by pairs 1151 * of variables for the equalities, followed by variables for the inequalities. 1152 * We then simply pick a feasible solution and compute s using (*). 1153 * 1154 * Note that we simply pick any valid direction and make no attempt 1155 * to pick a "good" or even the "best" valid direction. 1156 */ 1157static struct isl_vec *valid_direction( 1158 struct isl_basic_set *bset1, struct isl_basic_set *bset2) 1159{ 1160 struct isl_basic_set *lp; 1161 struct isl_tab *tab; 1162 struct isl_vec *sample = NULL; 1163 struct isl_vec *dir; 1164 unsigned d; 1165 int i; 1166 int n; 1167 1168 if (!bset1 || !bset2) 1169 goto error; 1170 lp = valid_direction_lp(isl_basic_set_copy(bset1), 1171 isl_basic_set_copy(bset2)); 1172 tab = isl_tab_from_basic_set(lp, 0); 1173 sample = isl_tab_get_sample_value(tab); 1174 isl_tab_free(tab); 1175 isl_basic_set_free(lp); 1176 if (!sample) 1177 goto error; 1178 d = isl_basic_set_total_dim(bset1); 1179 dir = isl_vec_alloc(bset1->ctx, 1 + d); 1180 if (!dir) 1181 goto error; 1182 isl_seq_clr(dir->block.data + 1, dir->size - 1); 1183 n = 1; 1184 /* positivity constraint 1 >= 0 */ 1185 isl_int_set(dir->block.data[0], sample->block.data[n]); n++; 1186 for (i = 0; i < bset1->n_eq; ++i) { 1187 isl_int_sub(sample->block.data[n], 1188 sample->block.data[n], sample->block.data[n+1]); 1189 isl_seq_combine(dir->block.data, 1190 bset1->ctx->one, dir->block.data, 1191 sample->block.data[n], bset1->eq[i], 1 + d); 1192 1193 n += 2; 1194 } 1195 for (i = 0; i < bset1->n_ineq; ++i) 1196 isl_seq_combine(dir->block.data, 1197 bset1->ctx->one, dir->block.data, 1198 sample->block.data[n++], bset1->ineq[i], 1 + d); 1199 isl_vec_free(sample); 1200 isl_seq_normalize(bset1->ctx, dir->el, dir->size); 1201 isl_basic_set_free(bset1); 1202 isl_basic_set_free(bset2); 1203 return dir; 1204error: 1205 isl_vec_free(sample); 1206 isl_basic_set_free(bset1); 1207 isl_basic_set_free(bset2); 1208 return NULL; 1209} 1210 1211/* Given a polyhedron b_i + A_i x >= 0 and a map T = S^{-1}, 1212 * compute b_i' + A_i' x' >= 0, with 1213 * 1214 * [ b_i A_i ] [ y' ] [ y' ] 1215 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0 1216 * 1217 * In particular, add the "positivity constraint" and then perform 1218 * the mapping. 1219 */ 1220static struct isl_basic_set *homogeneous_map(struct isl_basic_set *bset, 1221 struct isl_mat *T) 1222{ 1223 int k; 1224 1225 if (!bset) 1226 goto error; 1227 bset = isl_basic_set_extend_constraints(bset, 0, 1); 1228 k = isl_basic_set_alloc_inequality(bset); 1229 if (k < 0) 1230 goto error; 1231 isl_seq_clr(bset->ineq[k] + 1, isl_basic_set_total_dim(bset)); 1232 isl_int_set_si(bset->ineq[k][0], 1); 1233 bset = isl_basic_set_preimage(bset, T); 1234 return bset; 1235error: 1236 isl_mat_free(T); 1237 isl_basic_set_free(bset); 1238 return NULL; 1239} 1240 1241/* Compute the convex hull of a pair of basic sets without any parameters or 1242 * integer divisions, where the convex hull is known to be pointed, 1243 * but the basic sets may be unbounded. 1244 * 1245 * We turn this problem into the computation of a convex hull of a pair 1246 * _bounded_ polyhedra by "changing the direction of the homogeneous 1247 * dimension". This idea is due to Matthias Koeppe. 1248 * 1249 * Consider the cones in homogeneous space that correspond to the 1250 * input polyhedra. The rays of these cones are also rays of the 1251 * polyhedra if the coordinate that corresponds to the homogeneous 1252 * dimension is zero. That is, if the inner product of the rays 1253 * with the homogeneous direction is zero. 1254 * The cones in the homogeneous space can also be considered to 1255 * correspond to other pairs of polyhedra by chosing a different 1256 * homogeneous direction. To ensure that both of these polyhedra 1257 * are bounded, we need to make sure that all rays of the cones 1258 * correspond to vertices and not to rays. 1259 * Let s be a direction such that <s, r> > 0 for all rays r of both cones. 1260 * Then using s as a homogeneous direction, we obtain a pair of polytopes. 1261 * The vector s is computed in valid_direction. 1262 * 1263 * Note that we need to consider _all_ rays of the cones and not just 1264 * the rays that correspond to rays in the polyhedra. If we were to 1265 * only consider those rays and turn them into vertices, then we 1266 * may inadvertently turn some vertices into rays. 1267 * 1268 * The standard homogeneous direction is the unit vector in the 0th coordinate. 1269 * We therefore transform the two polyhedra such that the selected 1270 * direction is mapped onto this standard direction and then proceed 1271 * with the normal computation. 1272 * Let S be a non-singular square matrix with s as its first row, 1273 * then we want to map the polyhedra to the space 1274 * 1275 * [ y' ] [ y ] [ y ] [ y' ] 1276 * [ x' ] = S [ x ] i.e., [ x ] = S^{-1} [ x' ] 1277 * 1278 * We take S to be the unimodular completion of s to limit the growth 1279 * of the coefficients in the following computations. 1280 * 1281 * Let b_i + A_i x >= 0 be the constraints of polyhedron i. 1282 * We first move to the homogeneous dimension 1283 * 1284 * b_i y + A_i x >= 0 [ b_i A_i ] [ y ] [ 0 ] 1285 * y >= 0 or [ 1 0 ] [ x ] >= [ 0 ] 1286 * 1287 * Then we change directoin 1288 * 1289 * [ b_i A_i ] [ y' ] [ y' ] 1290 * [ 1 0 ] S^{-1} [ x' ] >= 0 or [ b_i' A_i' ] [ x' ] >= 0 1291 * 1292 * Then we compute the convex hull of the polytopes b_i' + A_i' x' >= 0 1293 * resulting in b' + A' x' >= 0, which we then convert back 1294 * 1295 * [ y ] [ y ] 1296 * [ b' A' ] S [ x ] >= 0 or [ b A ] [ x ] >= 0 1297 * 1298 * The polyhedron b + A x >= 0 is then the convex hull of the input polyhedra. 1299 */ 1300static struct isl_basic_set *convex_hull_pair_pointed( 1301 struct isl_basic_set *bset1, struct isl_basic_set *bset2) 1302{ 1303 struct isl_ctx *ctx = NULL; 1304 struct isl_vec *dir = NULL; 1305 struct isl_mat *T = NULL; 1306 struct isl_mat *T2 = NULL; 1307 struct isl_basic_set *hull; 1308 struct isl_set *set; 1309 1310 if (!bset1 || !bset2) 1311 goto error; 1312 ctx = bset1->ctx; 1313 dir = valid_direction(isl_basic_set_copy(bset1), 1314 isl_basic_set_copy(bset2)); 1315 if (!dir) 1316 goto error; 1317 T = isl_mat_alloc(bset1->ctx, dir->size, dir->size); 1318 if (!T) 1319 goto error; 1320 isl_seq_cpy(T->row[0], dir->block.data, dir->size); 1321 T = isl_mat_unimodular_complete(T, 1); 1322 T2 = isl_mat_right_inverse(isl_mat_copy(T)); 1323 1324 bset1 = homogeneous_map(bset1, isl_mat_copy(T2)); 1325 bset2 = homogeneous_map(bset2, T2); 1326 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0); 1327 set = isl_set_add_basic_set(set, bset1); 1328 set = isl_set_add_basic_set(set, bset2); 1329 hull = uset_convex_hull(set); 1330 hull = isl_basic_set_preimage(hull, T); 1331 1332 isl_vec_free(dir); 1333 1334 return hull; 1335error: 1336 isl_vec_free(dir); 1337 isl_basic_set_free(bset1); 1338 isl_basic_set_free(bset2); 1339 return NULL; 1340} 1341 1342static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set); 1343static struct isl_basic_set *modulo_affine_hull( 1344 struct isl_set *set, struct isl_basic_set *affine_hull); 1345 1346/* Compute the convex hull of a pair of basic sets without any parameters or 1347 * integer divisions. 1348 * 1349 * This function is called from uset_convex_hull_unbounded, which 1350 * means that the complete convex hull is unbounded. Some pairs 1351 * of basic sets may still be bounded, though. 1352 * They may even lie inside a lower dimensional space, in which 1353 * case they need to be handled inside their affine hull since 1354 * the main algorithm assumes that the result is full-dimensional. 1355 * 1356 * If the convex hull of the two basic sets would have a non-trivial 1357 * lineality space, we first project out this lineality space. 1358 */ 1359static struct isl_basic_set *convex_hull_pair(struct isl_basic_set *bset1, 1360 struct isl_basic_set *bset2) 1361{ 1362 isl_basic_set *lin, *aff; 1363 int bounded1, bounded2; 1364 1365 if (bset1->ctx->opt->convex == ISL_CONVEX_HULL_FM) 1366 return convex_hull_pair_elim(bset1, bset2); 1367 1368 aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1), 1369 isl_basic_set_copy(bset2))); 1370 if (!aff) 1371 goto error; 1372 if (aff->n_eq != 0) 1373 return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff); 1374 isl_basic_set_free(aff); 1375 1376 bounded1 = isl_basic_set_is_bounded(bset1); 1377 bounded2 = isl_basic_set_is_bounded(bset2); 1378 1379 if (bounded1 < 0 || bounded2 < 0) 1380 goto error; 1381 1382 if (bounded1 && bounded2) 1383 uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2)); 1384 1385 if (bounded1 || bounded2) 1386 return convex_hull_pair_pointed(bset1, bset2); 1387 1388 lin = induced_lineality_space(isl_basic_set_copy(bset1), 1389 isl_basic_set_copy(bset2)); 1390 if (!lin) 1391 goto error; 1392 if (isl_basic_set_is_universe(lin)) { 1393 isl_basic_set_free(bset1); 1394 isl_basic_set_free(bset2); 1395 return lin; 1396 } 1397 if (lin->n_eq < isl_basic_set_total_dim(lin)) { 1398 struct isl_set *set; 1399 set = isl_set_alloc_space(isl_basic_set_get_space(bset1), 2, 0); 1400 set = isl_set_add_basic_set(set, bset1); 1401 set = isl_set_add_basic_set(set, bset2); 1402 return modulo_lineality(set, lin); 1403 } 1404 isl_basic_set_free(lin); 1405 1406 return convex_hull_pair_pointed(bset1, bset2); 1407error: 1408 isl_basic_set_free(bset1); 1409 isl_basic_set_free(bset2); 1410 return NULL; 1411} 1412 1413/* Compute the lineality space of a basic set. 1414 * We currently do not allow the basic set to have any divs. 1415 * We basically just drop the constants and turn every inequality 1416 * into an equality. 1417 */ 1418struct isl_basic_set *isl_basic_set_lineality_space(struct isl_basic_set *bset) 1419{ 1420 int i, k; 1421 struct isl_basic_set *lin = NULL; 1422 unsigned dim; 1423 1424 if (!bset) 1425 goto error; 1426 isl_assert(bset->ctx, bset->n_div == 0, goto error); 1427 dim = isl_basic_set_total_dim(bset); 1428 1429 lin = isl_basic_set_alloc_space(isl_basic_set_get_space(bset), 0, dim, 0); 1430 if (!lin) 1431 goto error; 1432 for (i = 0; i < bset->n_eq; ++i) { 1433 k = isl_basic_set_alloc_equality(lin); 1434 if (k < 0) 1435 goto error; 1436 isl_int_set_si(lin->eq[k][0], 0); 1437 isl_seq_cpy(lin->eq[k] + 1, bset->eq[i] + 1, dim); 1438 } 1439 lin = isl_basic_set_gauss(lin, NULL); 1440 if (!lin) 1441 goto error; 1442 for (i = 0; i < bset->n_ineq && lin->n_eq < dim; ++i) { 1443 k = isl_basic_set_alloc_equality(lin); 1444 if (k < 0) 1445 goto error; 1446 isl_int_set_si(lin->eq[k][0], 0); 1447 isl_seq_cpy(lin->eq[k] + 1, bset->ineq[i] + 1, dim); 1448 lin = isl_basic_set_gauss(lin, NULL); 1449 if (!lin) 1450 goto error; 1451 } 1452 isl_basic_set_free(bset); 1453 return lin; 1454error: 1455 isl_basic_set_free(lin); 1456 isl_basic_set_free(bset); 1457 return NULL; 1458} 1459 1460/* Compute the (linear) hull of the lineality spaces of the basic sets in the 1461 * "underlying" set "set". 1462 */ 1463static struct isl_basic_set *uset_combined_lineality_space(struct isl_set *set) 1464{ 1465 int i; 1466 struct isl_set *lin = NULL; 1467 1468 if (!set) 1469 return NULL; 1470 if (set->n == 0) { 1471 isl_space *dim = isl_set_get_space(set); 1472 isl_set_free(set); 1473 return isl_basic_set_empty(dim); 1474 } 1475 1476 lin = isl_set_alloc_space(isl_set_get_space(set), set->n, 0); 1477 for (i = 0; i < set->n; ++i) 1478 lin = isl_set_add_basic_set(lin, 1479 isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i]))); 1480 isl_set_free(set); 1481 return isl_set_affine_hull(lin); 1482} 1483 1484/* Compute the convex hull of a set without any parameters or 1485 * integer divisions. 1486 * In each step, we combined two basic sets until only one 1487 * basic set is left. 1488 * The input basic sets are assumed not to have a non-trivial 1489 * lineality space. If any of the intermediate results has 1490 * a non-trivial lineality space, it is projected out. 1491 */ 1492static struct isl_basic_set *uset_convex_hull_unbounded(struct isl_set *set) 1493{ 1494 struct isl_basic_set *convex_hull = NULL; 1495 1496 convex_hull = isl_set_copy_basic_set(set); 1497 set = isl_set_drop_basic_set(set, convex_hull); 1498 if (!set) 1499 goto error; 1500 while (set->n > 0) { 1501 struct isl_basic_set *t; 1502 t = isl_set_copy_basic_set(set); 1503 if (!t) 1504 goto error; 1505 set = isl_set_drop_basic_set(set, t); 1506 if (!set) 1507 goto error; 1508 convex_hull = convex_hull_pair(convex_hull, t); 1509 if (set->n == 0) 1510 break; 1511 t = isl_basic_set_lineality_space(isl_basic_set_copy(convex_hull)); 1512 if (!t) 1513 goto error; 1514 if (isl_basic_set_is_universe(t)) { 1515 isl_basic_set_free(convex_hull); 1516 convex_hull = t; 1517 break; 1518 } 1519 if (t->n_eq < isl_basic_set_total_dim(t)) { 1520 set = isl_set_add_basic_set(set, convex_hull); 1521 return modulo_lineality(set, t); 1522 } 1523 isl_basic_set_free(t); 1524 } 1525 isl_set_free(set); 1526 return convex_hull; 1527error: 1528 isl_set_free(set); 1529 isl_basic_set_free(convex_hull); 1530 return NULL; 1531} 1532 1533/* Compute an initial hull for wrapping containing a single initial 1534 * facet. 1535 * This function assumes that the given set is bounded. 1536 */ 1537static struct isl_basic_set *initial_hull(struct isl_basic_set *hull, 1538 struct isl_set *set) 1539{ 1540 struct isl_mat *bounds = NULL; 1541 unsigned dim; 1542 int k; 1543 1544 if (!hull) 1545 goto error; 1546 bounds = initial_facet_constraint(set); 1547 if (!bounds) 1548 goto error; 1549 k = isl_basic_set_alloc_inequality(hull); 1550 if (k < 0) 1551 goto error; 1552 dim = isl_set_n_dim(set); 1553 isl_assert(set->ctx, 1 + dim == bounds->n_col, goto error); 1554 isl_seq_cpy(hull->ineq[k], bounds->row[0], bounds->n_col); 1555 isl_mat_free(bounds); 1556 1557 return hull; 1558error: 1559 isl_basic_set_free(hull); 1560 isl_mat_free(bounds); 1561 return NULL; 1562} 1563 1564struct max_constraint { 1565 struct isl_mat *c; 1566 int count; 1567 int ineq; 1568}; 1569 1570static int max_constraint_equal(const void *entry, const void *val) 1571{ 1572 struct max_constraint *a = (struct max_constraint *)entry; 1573 isl_int *b = (isl_int *)val; 1574 1575 return isl_seq_eq(a->c->row[0] + 1, b, a->c->n_col - 1); 1576} 1577 1578static void update_constraint(struct isl_ctx *ctx, struct isl_hash_table *table, 1579 isl_int *con, unsigned len, int n, int ineq) 1580{ 1581 struct isl_hash_table_entry *entry; 1582 struct max_constraint *c; 1583 uint32_t c_hash; 1584 1585 c_hash = isl_seq_get_hash(con + 1, len); 1586 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal, 1587 con + 1, 0); 1588 if (!entry) 1589 return; 1590 c = entry->data; 1591 if (c->count < n) { 1592 isl_hash_table_remove(ctx, table, entry); 1593 return; 1594 } 1595 c->count++; 1596 if (isl_int_gt(c->c->row[0][0], con[0])) 1597 return; 1598 if (isl_int_eq(c->c->row[0][0], con[0])) { 1599 if (ineq) 1600 c->ineq = ineq; 1601 return; 1602 } 1603 c->c = isl_mat_cow(c->c); 1604 isl_int_set(c->c->row[0][0], con[0]); 1605 c->ineq = ineq; 1606} 1607 1608/* Check whether the constraint hash table "table" constains the constraint 1609 * "con". 1610 */ 1611static int has_constraint(struct isl_ctx *ctx, struct isl_hash_table *table, 1612 isl_int *con, unsigned len, int n) 1613{ 1614 struct isl_hash_table_entry *entry; 1615 struct max_constraint *c; 1616 uint32_t c_hash; 1617 1618 c_hash = isl_seq_get_hash(con + 1, len); 1619 entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal, 1620 con + 1, 0); 1621 if (!entry) 1622 return 0; 1623 c = entry->data; 1624 if (c->count < n) 1625 return 0; 1626 return isl_int_eq(c->c->row[0][0], con[0]); 1627} 1628 1629/* Check for inequality constraints of a basic set without equalities 1630 * such that the same or more stringent copies of the constraint appear 1631 * in all of the basic sets. Such constraints are necessarily facet 1632 * constraints of the convex hull. 1633 * 1634 * If the resulting basic set is by chance identical to one of 1635 * the basic sets in "set", then we know that this basic set contains 1636 * all other basic sets and is therefore the convex hull of set. 1637 * In this case we set *is_hull to 1. 1638 */ 1639static struct isl_basic_set *common_constraints(struct isl_basic_set *hull, 1640 struct isl_set *set, int *is_hull) 1641{ 1642 int i, j, s, n; 1643 int min_constraints; 1644 int best; 1645 struct max_constraint *constraints = NULL; 1646 struct isl_hash_table *table = NULL; 1647 unsigned total; 1648 1649 *is_hull = 0; 1650 1651 for (i = 0; i < set->n; ++i) 1652 if (set->p[i]->n_eq == 0) 1653 break; 1654 if (i >= set->n) 1655 return hull; 1656 min_constraints = set->p[i]->n_ineq; 1657 best = i; 1658 for (i = best + 1; i < set->n; ++i) { 1659 if (set->p[i]->n_eq != 0) 1660 continue; 1661 if (set->p[i]->n_ineq >= min_constraints) 1662 continue; 1663 min_constraints = set->p[i]->n_ineq; 1664 best = i; 1665 } 1666 constraints = isl_calloc_array(hull->ctx, struct max_constraint, 1667 min_constraints); 1668 if (!constraints) 1669 return hull; 1670 table = isl_alloc_type(hull->ctx, struct isl_hash_table); 1671 if (isl_hash_table_init(hull->ctx, table, min_constraints)) 1672 goto error; 1673 1674 total = isl_space_dim(set->dim, isl_dim_all); 1675 for (i = 0; i < set->p[best]->n_ineq; ++i) { 1676 constraints[i].c = isl_mat_sub_alloc6(hull->ctx, 1677 set->p[best]->ineq + i, 0, 1, 0, 1 + total); 1678 if (!constraints[i].c) 1679 goto error; 1680 constraints[i].ineq = 1; 1681 } 1682 for (i = 0; i < min_constraints; ++i) { 1683 struct isl_hash_table_entry *entry; 1684 uint32_t c_hash; 1685 c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total); 1686 entry = isl_hash_table_find(hull->ctx, table, c_hash, 1687 max_constraint_equal, constraints[i].c->row[0] + 1, 1); 1688 if (!entry) 1689 goto error; 1690 isl_assert(hull->ctx, !entry->data, goto error); 1691 entry->data = &constraints[i]; 1692 } 1693 1694 n = 0; 1695 for (s = 0; s < set->n; ++s) { 1696 if (s == best) 1697 continue; 1698 1699 for (i = 0; i < set->p[s]->n_eq; ++i) { 1700 isl_int *eq = set->p[s]->eq[i]; 1701 for (j = 0; j < 2; ++j) { 1702 isl_seq_neg(eq, eq, 1 + total); 1703 update_constraint(hull->ctx, table, 1704 eq, total, n, 0); 1705 } 1706 } 1707 for (i = 0; i < set->p[s]->n_ineq; ++i) { 1708 isl_int *ineq = set->p[s]->ineq[i]; 1709 update_constraint(hull->ctx, table, ineq, total, n, 1710 set->p[s]->n_eq == 0); 1711 } 1712 ++n; 1713 } 1714 1715 for (i = 0; i < min_constraints; ++i) { 1716 if (constraints[i].count < n) 1717 continue; 1718 if (!constraints[i].ineq) 1719 continue; 1720 j = isl_basic_set_alloc_inequality(hull); 1721 if (j < 0) 1722 goto error; 1723 isl_seq_cpy(hull->ineq[j], constraints[i].c->row[0], 1 + total); 1724 } 1725 1726 for (s = 0; s < set->n; ++s) { 1727 if (set->p[s]->n_eq) 1728 continue; 1729 if (set->p[s]->n_ineq != hull->n_ineq) 1730 continue; 1731 for (i = 0; i < set->p[s]->n_ineq; ++i) { 1732 isl_int *ineq = set->p[s]->ineq[i]; 1733 if (!has_constraint(hull->ctx, table, ineq, total, n)) 1734 break; 1735 } 1736 if (i == set->p[s]->n_ineq) 1737 *is_hull = 1; 1738 } 1739 1740 isl_hash_table_clear(table); 1741 for (i = 0; i < min_constraints; ++i) 1742 isl_mat_free(constraints[i].c); 1743 free(constraints); 1744 free(table); 1745 return hull; 1746error: 1747 isl_hash_table_clear(table); 1748 free(table); 1749 if (constraints) 1750 for (i = 0; i < min_constraints; ++i) 1751 isl_mat_free(constraints[i].c); 1752 free(constraints); 1753 return hull; 1754} 1755 1756/* Create a template for the convex hull of "set" and fill it up 1757 * obvious facet constraints, if any. If the result happens to 1758 * be the convex hull of "set" then *is_hull is set to 1. 1759 */ 1760static struct isl_basic_set *proto_hull(struct isl_set *set, int *is_hull) 1761{ 1762 struct isl_basic_set *hull; 1763 unsigned n_ineq; 1764 int i; 1765 1766 n_ineq = 1; 1767 for (i = 0; i < set->n; ++i) { 1768 n_ineq += set->p[i]->n_eq; 1769 n_ineq += set->p[i]->n_ineq; 1770 } 1771 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq); 1772 hull = isl_basic_set_set_rational(hull); 1773 if (!hull) 1774 return NULL; 1775 return common_constraints(hull, set, is_hull); 1776} 1777 1778static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set) 1779{ 1780 struct isl_basic_set *hull; 1781 int is_hull; 1782 1783 hull = proto_hull(set, &is_hull); 1784 if (hull && !is_hull) { 1785 if (hull->n_ineq == 0) 1786 hull = initial_hull(hull, set); 1787 hull = extend(hull, set); 1788 } 1789 isl_set_free(set); 1790 1791 return hull; 1792} 1793 1794/* Compute the convex hull of a set without any parameters or 1795 * integer divisions. Depending on whether the set is bounded, 1796 * we pass control to the wrapping based convex hull or 1797 * the Fourier-Motzkin elimination based convex hull. 1798 * We also handle a few special cases before checking the boundedness. 1799 */ 1800static struct isl_basic_set *uset_convex_hull(struct isl_set *set) 1801{ 1802 struct isl_basic_set *convex_hull = NULL; 1803 struct isl_basic_set *lin; 1804 1805 if (isl_set_n_dim(set) == 0) 1806 return convex_hull_0d(set); 1807 1808 set = isl_set_coalesce(set); 1809 set = isl_set_set_rational(set); 1810 1811 if (!set) 1812 goto error; 1813 if (!set) 1814 return NULL; 1815 if (set->n == 1) { 1816 convex_hull = isl_basic_set_copy(set->p[0]); 1817 isl_set_free(set); 1818 return convex_hull; 1819 } 1820 if (isl_set_n_dim(set) == 1) 1821 return convex_hull_1d(set); 1822 1823 if (isl_set_is_bounded(set) && 1824 set->ctx->opt->convex == ISL_CONVEX_HULL_WRAP) 1825 return uset_convex_hull_wrap(set); 1826 1827 lin = uset_combined_lineality_space(isl_set_copy(set)); 1828 if (!lin) 1829 goto error; 1830 if (isl_basic_set_is_universe(lin)) { 1831 isl_set_free(set); 1832 return lin; 1833 } 1834 if (lin->n_eq < isl_basic_set_total_dim(lin)) 1835 return modulo_lineality(set, lin); 1836 isl_basic_set_free(lin); 1837 1838 return uset_convex_hull_unbounded(set); 1839error: 1840 isl_set_free(set); 1841 isl_basic_set_free(convex_hull); 1842 return NULL; 1843} 1844 1845/* This is the core procedure, where "set" is a "pure" set, i.e., 1846 * without parameters or divs and where the convex hull of set is 1847 * known to be full-dimensional. 1848 */ 1849static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set) 1850{ 1851 struct isl_basic_set *convex_hull = NULL; 1852 1853 if (!set) 1854 goto error; 1855 1856 if (isl_set_n_dim(set) == 0) { 1857 convex_hull = isl_basic_set_universe(isl_space_copy(set->dim)); 1858 isl_set_free(set); 1859 convex_hull = isl_basic_set_set_rational(convex_hull); 1860 return convex_hull; 1861 } 1862 1863 set = isl_set_set_rational(set); 1864 set = isl_set_coalesce(set); 1865 if (!set) 1866 goto error; 1867 if (set->n == 1) { 1868 convex_hull = isl_basic_set_copy(set->p[0]); 1869 isl_set_free(set); 1870 convex_hull = isl_basic_map_remove_redundancies(convex_hull); 1871 return convex_hull; 1872 } 1873 if (isl_set_n_dim(set) == 1) 1874 return convex_hull_1d(set); 1875 1876 return uset_convex_hull_wrap(set); 1877error: 1878 isl_set_free(set); 1879 return NULL; 1880} 1881 1882/* Compute the convex hull of set "set" with affine hull "affine_hull", 1883 * We first remove the equalities (transforming the set), compute the 1884 * convex hull of the transformed set and then add the equalities back 1885 * (after performing the inverse transformation. 1886 */ 1887static struct isl_basic_set *modulo_affine_hull( 1888 struct isl_set *set, struct isl_basic_set *affine_hull) 1889{ 1890 struct isl_mat *T; 1891 struct isl_mat *T2; 1892 struct isl_basic_set *dummy; 1893 struct isl_basic_set *convex_hull; 1894 1895 dummy = isl_basic_set_remove_equalities( 1896 isl_basic_set_copy(affine_hull), &T, &T2); 1897 if (!dummy) 1898 goto error; 1899 isl_basic_set_free(dummy); 1900 set = isl_set_preimage(set, T); 1901 convex_hull = uset_convex_hull(set); 1902 convex_hull = isl_basic_set_preimage(convex_hull, T2); 1903 convex_hull = isl_basic_set_intersect(convex_hull, affine_hull); 1904 return convex_hull; 1905error: 1906 isl_basic_set_free(affine_hull); 1907 isl_set_free(set); 1908 return NULL; 1909} 1910 1911/* Compute the convex hull of a map. 1912 * 1913 * The implementation was inspired by "Extended Convex Hull" by Fukuda et al., 1914 * specifically, the wrapping of facets to obtain new facets. 1915 */ 1916struct isl_basic_map *isl_map_convex_hull(struct isl_map *map) 1917{ 1918 struct isl_basic_set *bset; 1919 struct isl_basic_map *model = NULL; 1920 struct isl_basic_set *affine_hull = NULL; 1921 struct isl_basic_map *convex_hull = NULL; 1922 struct isl_set *set = NULL; 1923 struct isl_ctx *ctx; 1924 1925 if (!map) 1926 goto error; 1927 1928 ctx = map->ctx; 1929 if (map->n == 0) { 1930 convex_hull = isl_basic_map_empty_like_map(map); 1931 isl_map_free(map); 1932 return convex_hull; 1933 } 1934 1935 map = isl_map_detect_equalities(map); 1936 map = isl_map_align_divs(map); 1937 if (!map) 1938 goto error; 1939 model = isl_basic_map_copy(map->p[0]); 1940 set = isl_map_underlying_set(map); 1941 if (!set) 1942 goto error; 1943 1944 affine_hull = isl_set_affine_hull(isl_set_copy(set)); 1945 if (!affine_hull) 1946 goto error; 1947 if (affine_hull->n_eq != 0) 1948 bset = modulo_affine_hull(set, affine_hull); 1949 else { 1950 isl_basic_set_free(affine_hull); 1951 bset = uset_convex_hull(set); 1952 } 1953 1954 convex_hull = isl_basic_map_overlying_set(bset, model); 1955 if (!convex_hull) 1956 return NULL; 1957 1958 ISL_F_SET(convex_hull, ISL_BASIC_MAP_NO_IMPLICIT); 1959 ISL_F_SET(convex_hull, ISL_BASIC_MAP_ALL_EQUALITIES); 1960 ISL_F_CLR(convex_hull, ISL_BASIC_MAP_RATIONAL); 1961 return convex_hull; 1962error: 1963 isl_set_free(set); 1964 isl_basic_map_free(model); 1965 return NULL; 1966} 1967 1968struct isl_basic_set *isl_set_convex_hull(struct isl_set *set) 1969{ 1970 return (struct isl_basic_set *) 1971 isl_map_convex_hull((struct isl_map *)set); 1972} 1973 1974__isl_give isl_basic_map *isl_map_polyhedral_hull(__isl_take isl_map *map) 1975{ 1976 isl_basic_map *hull; 1977 1978 hull = isl_map_convex_hull(map); 1979 return isl_basic_map_remove_divs(hull); 1980} 1981 1982__isl_give isl_basic_set *isl_set_polyhedral_hull(__isl_take isl_set *set) 1983{ 1984 return (isl_basic_set *)isl_map_polyhedral_hull((isl_map *)set); 1985} 1986 1987struct sh_data_entry { 1988 struct isl_hash_table *table; 1989 struct isl_tab *tab; 1990}; 1991 1992/* Holds the data needed during the simple hull computation. 1993 * In particular, 1994 * n the number of basic sets in the original set 1995 * hull_table a hash table of already computed constraints 1996 * in the simple hull 1997 * p for each basic set, 1998 * table a hash table of the constraints 1999 * tab the tableau corresponding to the basic set 2000 */ 2001struct sh_data { 2002 struct isl_ctx *ctx; 2003 unsigned n; 2004 struct isl_hash_table *hull_table; 2005 struct sh_data_entry p[1]; 2006}; 2007 2008static void sh_data_free(struct sh_data *data) 2009{ 2010 int i; 2011 2012 if (!data) 2013 return; 2014 isl_hash_table_free(data->ctx, data->hull_table); 2015 for (i = 0; i < data->n; ++i) { 2016 isl_hash_table_free(data->ctx, data->p[i].table); 2017 isl_tab_free(data->p[i].tab); 2018 } 2019 free(data); 2020} 2021 2022struct ineq_cmp_data { 2023 unsigned len; 2024 isl_int *p; 2025}; 2026 2027static int has_ineq(const void *entry, const void *val) 2028{ 2029 isl_int *row = (isl_int *)entry; 2030 struct ineq_cmp_data *v = (struct ineq_cmp_data *)val; 2031 2032 return isl_seq_eq(row + 1, v->p + 1, v->len) || 2033 isl_seq_is_neg(row + 1, v->p + 1, v->len); 2034} 2035 2036static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table, 2037 isl_int *ineq, unsigned len) 2038{ 2039 uint32_t c_hash; 2040 struct ineq_cmp_data v; 2041 struct isl_hash_table_entry *entry; 2042 2043 v.len = len; 2044 v.p = ineq; 2045 c_hash = isl_seq_get_hash(ineq + 1, len); 2046 entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1); 2047 if (!entry) 2048 return - 1; 2049 entry->data = ineq; 2050 return 0; 2051} 2052 2053/* Fill hash table "table" with the constraints of "bset". 2054 * Equalities are added as two inequalities. 2055 * The value in the hash table is a pointer to the (in)equality of "bset". 2056 */ 2057static int hash_basic_set(struct isl_hash_table *table, 2058 struct isl_basic_set *bset) 2059{ 2060 int i, j; 2061 unsigned dim = isl_basic_set_total_dim(bset); 2062 2063 for (i = 0; i < bset->n_eq; ++i) { 2064 for (j = 0; j < 2; ++j) { 2065 isl_seq_neg(bset->eq[i], bset->eq[i], 1 + dim); 2066 if (hash_ineq(bset->ctx, table, bset->eq[i], dim) < 0) 2067 return -1; 2068 } 2069 } 2070 for (i = 0; i < bset->n_ineq; ++i) { 2071 if (hash_ineq(bset->ctx, table, bset->ineq[i], dim) < 0) 2072 return -1; 2073 } 2074 return 0; 2075} 2076 2077static struct sh_data *sh_data_alloc(struct isl_set *set, unsigned n_ineq) 2078{ 2079 struct sh_data *data; 2080 int i; 2081 2082 data = isl_calloc(set->ctx, struct sh_data, 2083 sizeof(struct sh_data) + 2084 (set->n - 1) * sizeof(struct sh_data_entry)); 2085 if (!data) 2086 return NULL; 2087 data->ctx = set->ctx; 2088 data->n = set->n; 2089 data->hull_table = isl_hash_table_alloc(set->ctx, n_ineq); 2090 if (!data->hull_table) 2091 goto error; 2092 for (i = 0; i < set->n; ++i) { 2093 data->p[i].table = isl_hash_table_alloc(set->ctx, 2094 2 * set->p[i]->n_eq + set->p[i]->n_ineq); 2095 if (!data->p[i].table) 2096 goto error; 2097 if (hash_basic_set(data->p[i].table, set->p[i]) < 0) 2098 goto error; 2099 } 2100 return data; 2101error: 2102 sh_data_free(data); 2103 return NULL; 2104} 2105 2106/* Check if inequality "ineq" is a bound for basic set "j" or if 2107 * it can be relaxed (by increasing the constant term) to become 2108 * a bound for that basic set. In the latter case, the constant 2109 * term is updated. 2110 * Relaxation of the constant term is only allowed if "shift" is set. 2111 * 2112 * Return 1 if "ineq" is a bound 2113 * 0 if "ineq" may attain arbitrarily small values on basic set "j" 2114 * -1 if some error occurred 2115 */ 2116static int is_bound(struct sh_data *data, struct isl_set *set, int j, 2117 isl_int *ineq, int shift) 2118{ 2119 enum isl_lp_result res; 2120 isl_int opt; 2121 2122 if (!data->p[j].tab) { 2123 data->p[j].tab = isl_tab_from_basic_set(set->p[j], 0); 2124 if (!data->p[j].tab) 2125 return -1; 2126 } 2127 2128 isl_int_init(opt); 2129 2130 res = isl_tab_min(data->p[j].tab, ineq, data->ctx->one, 2131 &opt, NULL, 0); 2132 if (res == isl_lp_ok && isl_int_is_neg(opt)) { 2133 if (shift) 2134 isl_int_sub(ineq[0], ineq[0], opt); 2135 else 2136 res = isl_lp_unbounded; 2137 } 2138 2139 isl_int_clear(opt); 2140 2141 return (res == isl_lp_ok || res == isl_lp_empty) ? 1 : 2142 res == isl_lp_unbounded ? 0 : -1; 2143} 2144 2145/* Check if inequality "ineq" from basic set "i" is or can be relaxed to 2146 * become a bound on the whole set. If so, add the (relaxed) inequality 2147 * to "hull". Relaxation is only allowed if "shift" is set. 2148 * 2149 * We first check if "hull" already contains a translate of the inequality. 2150 * If so, we are done. 2151 * Then, we check if any of the previous basic sets contains a translate 2152 * of the inequality. If so, then we have already considered this 2153 * inequality and we are done. 2154 * Otherwise, for each basic set other than "i", we check if the inequality 2155 * is a bound on the basic set. 2156 * For previous basic sets, we know that they do not contain a translate 2157 * of the inequality, so we directly call is_bound. 2158 * For following basic sets, we first check if a translate of the 2159 * inequality appears in its description and if so directly update 2160 * the inequality accordingly. 2161 */ 2162static struct isl_basic_set *add_bound(struct isl_basic_set *hull, 2163 struct sh_data *data, struct isl_set *set, int i, isl_int *ineq, 2164 int shift) 2165{ 2166 uint32_t c_hash; 2167 struct ineq_cmp_data v; 2168 struct isl_hash_table_entry *entry; 2169 int j, k; 2170 2171 if (!hull) 2172 return NULL; 2173 2174 v.len = isl_basic_set_total_dim(hull); 2175 v.p = ineq; 2176 c_hash = isl_seq_get_hash(ineq + 1, v.len); 2177 2178 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash, 2179 has_ineq, &v, 0); 2180 if (entry) 2181 return hull; 2182 2183 for (j = 0; j < i; ++j) { 2184 entry = isl_hash_table_find(hull->ctx, data->p[j].table, 2185 c_hash, has_ineq, &v, 0); 2186 if (entry) 2187 break; 2188 } 2189 if (j < i) 2190 return hull; 2191 2192 k = isl_basic_set_alloc_inequality(hull); 2193 isl_seq_cpy(hull->ineq[k], ineq, 1 + v.len); 2194 if (k < 0) 2195 goto error; 2196 2197 for (j = 0; j < i; ++j) { 2198 int bound; 2199 bound = is_bound(data, set, j, hull->ineq[k], shift); 2200 if (bound < 0) 2201 goto error; 2202 if (!bound) 2203 break; 2204 } 2205 if (j < i) { 2206 isl_basic_set_free_inequality(hull, 1); 2207 return hull; 2208 } 2209 2210 for (j = i + 1; j < set->n; ++j) { 2211 int bound, neg; 2212 isl_int *ineq_j; 2213 entry = isl_hash_table_find(hull->ctx, data->p[j].table, 2214 c_hash, has_ineq, &v, 0); 2215 if (entry) { 2216 ineq_j = entry->data; 2217 neg = isl_seq_is_neg(ineq_j + 1, 2218 hull->ineq[k] + 1, v.len); 2219 if (neg) 2220 isl_int_neg(ineq_j[0], ineq_j[0]); 2221 if (isl_int_gt(ineq_j[0], hull->ineq[k][0])) 2222 isl_int_set(hull->ineq[k][0], ineq_j[0]); 2223 if (neg) 2224 isl_int_neg(ineq_j[0], ineq_j[0]); 2225 continue; 2226 } 2227 bound = is_bound(data, set, j, hull->ineq[k], shift); 2228 if (bound < 0) 2229 goto error; 2230 if (!bound) 2231 break; 2232 } 2233 if (j < set->n) { 2234 isl_basic_set_free_inequality(hull, 1); 2235 return hull; 2236 } 2237 2238 entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash, 2239 has_ineq, &v, 1); 2240 if (!entry) 2241 goto error; 2242 entry->data = hull->ineq[k]; 2243 2244 return hull; 2245error: 2246 isl_basic_set_free(hull); 2247 return NULL; 2248} 2249 2250/* Check if any inequality from basic set "i" is or can be relaxed to 2251 * become a bound on the whole set. If so, add the (relaxed) inequality 2252 * to "hull". Relaxation is only allowed if "shift" is set. 2253 */ 2254static struct isl_basic_set *add_bounds(struct isl_basic_set *bset, 2255 struct sh_data *data, struct isl_set *set, int i, int shift) 2256{ 2257 int j, k; 2258 unsigned dim = isl_basic_set_total_dim(bset); 2259 2260 for (j = 0; j < set->p[i]->n_eq; ++j) { 2261 for (k = 0; k < 2; ++k) { 2262 isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim); 2263 bset = add_bound(bset, data, set, i, set->p[i]->eq[j], 2264 shift); 2265 } 2266 } 2267 for (j = 0; j < set->p[i]->n_ineq; ++j) 2268 bset = add_bound(bset, data, set, i, set->p[i]->ineq[j], shift); 2269 return bset; 2270} 2271 2272/* Compute a superset of the convex hull of set that is described 2273 * by only (translates of) the constraints in the constituents of set. 2274 * Translation is only allowed if "shift" is set. 2275 */ 2276static __isl_give isl_basic_set *uset_simple_hull(__isl_take isl_set *set, 2277 int shift) 2278{ 2279 struct sh_data *data = NULL; 2280 struct isl_basic_set *hull = NULL; 2281 unsigned n_ineq; 2282 int i; 2283 2284 if (!set) 2285 return NULL; 2286 2287 n_ineq = 0; 2288 for (i = 0; i < set->n; ++i) { 2289 if (!set->p[i]) 2290 goto error; 2291 n_ineq += 2 * set->p[i]->n_eq + set->p[i]->n_ineq; 2292 } 2293 2294 hull = isl_basic_set_alloc_space(isl_space_copy(set->dim), 0, 0, n_ineq); 2295 if (!hull) 2296 goto error; 2297 2298 data = sh_data_alloc(set, n_ineq); 2299 if (!data) 2300 goto error; 2301 2302 for (i = 0; i < set->n; ++i) 2303 hull = add_bounds(hull, data, set, i, shift); 2304 2305 sh_data_free(data); 2306 isl_set_free(set); 2307 2308 return hull; 2309error: 2310 sh_data_free(data); 2311 isl_basic_set_free(hull); 2312 isl_set_free(set); 2313 return NULL; 2314} 2315 2316/* Compute a superset of the convex hull of map that is described 2317 * by only (translates of) the constraints in the constituents of map. 2318 * Translation is only allowed if "shift" is set. 2319 */ 2320static __isl_give isl_basic_map *map_simple_hull(__isl_take isl_map *map, 2321 int shift) 2322{ 2323 struct isl_set *set = NULL; 2324 struct isl_basic_map *model = NULL; 2325 struct isl_basic_map *hull; 2326 struct isl_basic_map *affine_hull; 2327 struct isl_basic_set *bset = NULL; 2328 2329 if (!map) 2330 return NULL; 2331 if (map->n == 0) { 2332 hull = isl_basic_map_empty_like_map(map); 2333 isl_map_free(map); 2334 return hull; 2335 } 2336 if (map->n == 1) { 2337 hull = isl_basic_map_copy(map->p[0]); 2338 isl_map_free(map); 2339 return hull; 2340 } 2341 2342 map = isl_map_detect_equalities(map); 2343 affine_hull = isl_map_affine_hull(isl_map_copy(map)); 2344 map = isl_map_align_divs(map); 2345 model = map ? isl_basic_map_copy(map->p[0]) : NULL; 2346 2347 set = isl_map_underlying_set(map); 2348 2349 bset = uset_simple_hull(set, shift); 2350 2351 hull = isl_basic_map_overlying_set(bset, model); 2352 2353 hull = isl_basic_map_intersect(hull, affine_hull); 2354 hull = isl_basic_map_remove_redundancies(hull); 2355 2356 if (!hull) 2357 return NULL; 2358 ISL_F_SET(hull, ISL_BASIC_MAP_NO_IMPLICIT); 2359 ISL_F_SET(hull, ISL_BASIC_MAP_ALL_EQUALITIES); 2360 2361 return hull; 2362} 2363 2364/* Compute a superset of the convex hull of map that is described 2365 * by only translates of the constraints in the constituents of map. 2366 */ 2367__isl_give isl_basic_map *isl_map_simple_hull(__isl_take isl_map *map) 2368{ 2369 return map_simple_hull(map, 1); 2370} 2371 2372struct isl_basic_set *isl_set_simple_hull(struct isl_set *set) 2373{ 2374 return (struct isl_basic_set *) 2375 isl_map_simple_hull((struct isl_map *)set); 2376} 2377 2378/* Compute a superset of the convex hull of map that is described 2379 * by only the constraints in the constituents of map. 2380 */ 2381__isl_give isl_basic_map *isl_map_unshifted_simple_hull( 2382 __isl_take isl_map *map) 2383{ 2384 return map_simple_hull(map, 0); 2385} 2386 2387__isl_give isl_basic_set *isl_set_unshifted_simple_hull( 2388 __isl_take isl_set *set) 2389{ 2390 return isl_map_unshifted_simple_hull(set); 2391} 2392 2393/* Given a set "set", return parametric bounds on the dimension "dim". 2394 */ 2395static struct isl_basic_set *set_bounds(struct isl_set *set, int dim) 2396{ 2397 unsigned set_dim = isl_set_dim(set, isl_dim_set); 2398 set = isl_set_copy(set); 2399 set = isl_set_eliminate_dims(set, dim + 1, set_dim - (dim + 1)); 2400 set = isl_set_eliminate_dims(set, 0, dim); 2401 return isl_set_convex_hull(set); 2402} 2403 2404/* Computes a "simple hull" and then check if each dimension in the 2405 * resulting hull is bounded by a symbolic constant. If not, the 2406 * hull is intersected with the corresponding bounds on the whole set. 2407 */ 2408struct isl_basic_set *isl_set_bounded_simple_hull(struct isl_set *set) 2409{ 2410 int i, j; 2411 struct isl_basic_set *hull; 2412 unsigned nparam, left; 2413 int removed_divs = 0; 2414 2415 hull = isl_set_simple_hull(isl_set_copy(set)); 2416 if (!hull) 2417 goto error; 2418 2419 nparam = isl_basic_set_dim(hull, isl_dim_param); 2420 for (i = 0; i < isl_basic_set_dim(hull, isl_dim_set); ++i) { 2421 int lower = 0, upper = 0; 2422 struct isl_basic_set *bounds; 2423 2424 left = isl_basic_set_total_dim(hull) - nparam - i - 1; 2425 for (j = 0; j < hull->n_eq; ++j) { 2426 if (isl_int_is_zero(hull->eq[j][1 + nparam + i])) 2427 continue; 2428 if (isl_seq_first_non_zero(hull->eq[j]+1+nparam+i+1, 2429 left) == -1) 2430 break; 2431 } 2432 if (j < hull->n_eq) 2433 continue; 2434 2435 for (j = 0; j < hull->n_ineq; ++j) { 2436 if (isl_int_is_zero(hull->ineq[j][1 + nparam + i])) 2437 continue; 2438 if (isl_seq_first_non_zero(hull->ineq[j]+1+nparam+i+1, 2439 left) != -1 || 2440 isl_seq_first_non_zero(hull->ineq[j]+1+nparam, 2441 i) != -1) 2442 continue; 2443 if (isl_int_is_pos(hull->ineq[j][1 + nparam + i])) 2444 lower = 1; 2445 else 2446 upper = 1; 2447 if (lower && upper) 2448 break; 2449 } 2450 2451 if (lower && upper) 2452 continue; 2453 2454 if (!removed_divs) { 2455 set = isl_set_remove_divs(set); 2456 if (!set) 2457 goto error; 2458 removed_divs = 1; 2459 } 2460 bounds = set_bounds(set, i); 2461 hull = isl_basic_set_intersect(hull, bounds); 2462 if (!hull) 2463 goto error; 2464 } 2465 2466 isl_set_free(set); 2467 return hull; 2468error: 2469 isl_set_free(set); 2470 return NULL; 2471} 2472