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/ilp.h> 13#include "isl_sample.h" 14#include <isl/seq.h> 15#include "isl_equalities.h" 16#include <isl_aff_private.h> 17#include <isl_local_space_private.h> 18#include <isl_mat_private.h> 19#include <isl_val_private.h> 20 21/* Given a basic set "bset", construct a basic set U such that for 22 * each element x in U, the whole unit box positioned at x is inside 23 * the given basic set. 24 * Note that U may not contain all points that satisfy this property. 25 * 26 * We simply add the sum of all negative coefficients to the constant 27 * term. This ensures that if x satisfies the resulting constraints, 28 * then x plus any sum of unit vectors satisfies the original constraints. 29 */ 30static struct isl_basic_set *unit_box_base_points(struct isl_basic_set *bset) 31{ 32 int i, j, k; 33 struct isl_basic_set *unit_box = NULL; 34 unsigned total; 35 36 if (!bset) 37 goto error; 38 39 if (bset->n_eq != 0) { 40 unit_box = isl_basic_set_empty_like(bset); 41 isl_basic_set_free(bset); 42 return unit_box; 43 } 44 45 total = isl_basic_set_total_dim(bset); 46 unit_box = isl_basic_set_alloc_space(isl_basic_set_get_space(bset), 47 0, 0, bset->n_ineq); 48 49 for (i = 0; i < bset->n_ineq; ++i) { 50 k = isl_basic_set_alloc_inequality(unit_box); 51 if (k < 0) 52 goto error; 53 isl_seq_cpy(unit_box->ineq[k], bset->ineq[i], 1 + total); 54 for (j = 0; j < total; ++j) { 55 if (isl_int_is_nonneg(unit_box->ineq[k][1 + j])) 56 continue; 57 isl_int_add(unit_box->ineq[k][0], 58 unit_box->ineq[k][0], unit_box->ineq[k][1 + j]); 59 } 60 } 61 62 isl_basic_set_free(bset); 63 return unit_box; 64error: 65 isl_basic_set_free(bset); 66 isl_basic_set_free(unit_box); 67 return NULL; 68} 69 70/* Find an integer point in "bset", preferably one that is 71 * close to minimizing "f". 72 * 73 * We first check if we can easily put unit boxes inside bset. 74 * If so, we take the best base point of any of the unit boxes we can find 75 * and round it up to the nearest integer. 76 * If not, we simply pick any integer point in "bset". 77 */ 78static struct isl_vec *initial_solution(struct isl_basic_set *bset, isl_int *f) 79{ 80 enum isl_lp_result res; 81 struct isl_basic_set *unit_box; 82 struct isl_vec *sol; 83 84 unit_box = unit_box_base_points(isl_basic_set_copy(bset)); 85 86 res = isl_basic_set_solve_lp(unit_box, 0, f, bset->ctx->one, 87 NULL, NULL, &sol); 88 if (res == isl_lp_ok) { 89 isl_basic_set_free(unit_box); 90 return isl_vec_ceil(sol); 91 } 92 93 isl_basic_set_free(unit_box); 94 95 return isl_basic_set_sample_vec(isl_basic_set_copy(bset)); 96} 97 98/* Restrict "bset" to those points with values for f in the interval [l, u]. 99 */ 100static struct isl_basic_set *add_bounds(struct isl_basic_set *bset, 101 isl_int *f, isl_int l, isl_int u) 102{ 103 int k; 104 unsigned total; 105 106 total = isl_basic_set_total_dim(bset); 107 bset = isl_basic_set_extend_constraints(bset, 0, 2); 108 109 k = isl_basic_set_alloc_inequality(bset); 110 if (k < 0) 111 goto error; 112 isl_seq_cpy(bset->ineq[k], f, 1 + total); 113 isl_int_sub(bset->ineq[k][0], bset->ineq[k][0], l); 114 115 k = isl_basic_set_alloc_inequality(bset); 116 if (k < 0) 117 goto error; 118 isl_seq_neg(bset->ineq[k], f, 1 + total); 119 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], u); 120 121 return bset; 122error: 123 isl_basic_set_free(bset); 124 return NULL; 125} 126 127/* Find an integer point in "bset" that minimizes f (in any) such that 128 * the value of f lies inside the interval [l, u]. 129 * Return this integer point if it can be found. 130 * Otherwise, return sol. 131 * 132 * We perform a number of steps until l > u. 133 * In each step, we look for an integer point with value in either 134 * the whole interval [l, u] or half of the interval [l, l+floor(u-l-1/2)]. 135 * The choice depends on whether we have found an integer point in the 136 * previous step. If so, we look for the next point in half of the remaining 137 * interval. 138 * If we find a point, the current solution is updated and u is set 139 * to its value minus 1. 140 * If no point can be found, we update l to the upper bound of the interval 141 * we checked (u or l+floor(u-l-1/2)) plus 1. 142 */ 143static struct isl_vec *solve_ilp_search(struct isl_basic_set *bset, 144 isl_int *f, isl_int *opt, struct isl_vec *sol, isl_int l, isl_int u) 145{ 146 isl_int tmp; 147 int divide = 1; 148 149 isl_int_init(tmp); 150 151 while (isl_int_le(l, u)) { 152 struct isl_basic_set *slice; 153 struct isl_vec *sample; 154 155 if (!divide) 156 isl_int_set(tmp, u); 157 else { 158 isl_int_sub(tmp, u, l); 159 isl_int_fdiv_q_ui(tmp, tmp, 2); 160 isl_int_add(tmp, tmp, l); 161 } 162 slice = add_bounds(isl_basic_set_copy(bset), f, l, tmp); 163 sample = isl_basic_set_sample_vec(slice); 164 if (!sample) { 165 isl_vec_free(sol); 166 sol = NULL; 167 break; 168 } 169 if (sample->size > 0) { 170 isl_vec_free(sol); 171 sol = sample; 172 isl_seq_inner_product(f, sol->el, sol->size, opt); 173 isl_int_sub_ui(u, *opt, 1); 174 divide = 1; 175 } else { 176 isl_vec_free(sample); 177 if (!divide) 178 break; 179 isl_int_add_ui(l, tmp, 1); 180 divide = 0; 181 } 182 } 183 184 isl_int_clear(tmp); 185 186 return sol; 187} 188 189/* Find an integer point in "bset" that minimizes f (if any). 190 * If sol_p is not NULL then the integer point is returned in *sol_p. 191 * The optimal value of f is returned in *opt. 192 * 193 * The algorithm maintains a currently best solution and an interval [l, u] 194 * of values of f for which integer solutions could potentially still be found. 195 * The initial value of the best solution so far is any solution. 196 * The initial value of l is minimal value of f over the rationals 197 * (rounded up to the nearest integer). 198 * The initial value of u is the value of f at the initial solution minus 1. 199 * 200 * We then call solve_ilp_search to perform a binary search on the interval. 201 */ 202static enum isl_lp_result solve_ilp(struct isl_basic_set *bset, 203 isl_int *f, isl_int *opt, 204 struct isl_vec **sol_p) 205{ 206 enum isl_lp_result res; 207 isl_int l, u; 208 struct isl_vec *sol; 209 210 res = isl_basic_set_solve_lp(bset, 0, f, bset->ctx->one, 211 opt, NULL, &sol); 212 if (res == isl_lp_ok && isl_int_is_one(sol->el[0])) { 213 if (sol_p) 214 *sol_p = sol; 215 else 216 isl_vec_free(sol); 217 return isl_lp_ok; 218 } 219 isl_vec_free(sol); 220 if (res == isl_lp_error || res == isl_lp_empty) 221 return res; 222 223 sol = initial_solution(bset, f); 224 if (!sol) 225 return isl_lp_error; 226 if (sol->size == 0) { 227 isl_vec_free(sol); 228 return isl_lp_empty; 229 } 230 if (res == isl_lp_unbounded) { 231 isl_vec_free(sol); 232 return isl_lp_unbounded; 233 } 234 235 isl_int_init(l); 236 isl_int_init(u); 237 238 isl_int_set(l, *opt); 239 240 isl_seq_inner_product(f, sol->el, sol->size, opt); 241 isl_int_sub_ui(u, *opt, 1); 242 243 sol = solve_ilp_search(bset, f, opt, sol, l, u); 244 if (!sol) 245 res = isl_lp_error; 246 247 isl_int_clear(l); 248 isl_int_clear(u); 249 250 if (sol_p) 251 *sol_p = sol; 252 else 253 isl_vec_free(sol); 254 255 return res; 256} 257 258static enum isl_lp_result solve_ilp_with_eq(struct isl_basic_set *bset, int max, 259 isl_int *f, isl_int *opt, 260 struct isl_vec **sol_p) 261{ 262 unsigned dim; 263 enum isl_lp_result res; 264 struct isl_mat *T = NULL; 265 struct isl_vec *v; 266 267 bset = isl_basic_set_copy(bset); 268 dim = isl_basic_set_total_dim(bset); 269 v = isl_vec_alloc(bset->ctx, 1 + dim); 270 if (!v) 271 goto error; 272 isl_seq_cpy(v->el, f, 1 + dim); 273 bset = isl_basic_set_remove_equalities(bset, &T, NULL); 274 v = isl_vec_mat_product(v, isl_mat_copy(T)); 275 if (!v) 276 goto error; 277 res = isl_basic_set_solve_ilp(bset, max, v->el, opt, sol_p); 278 isl_vec_free(v); 279 if (res == isl_lp_ok && sol_p) { 280 *sol_p = isl_mat_vec_product(T, *sol_p); 281 if (!*sol_p) 282 res = isl_lp_error; 283 } else 284 isl_mat_free(T); 285 isl_basic_set_free(bset); 286 return res; 287error: 288 isl_mat_free(T); 289 isl_basic_set_free(bset); 290 return isl_lp_error; 291} 292 293/* Find an integer point in "bset" that minimizes (or maximizes if max is set) 294 * f (if any). 295 * If sol_p is not NULL then the integer point is returned in *sol_p. 296 * The optimal value of f is returned in *opt. 297 * 298 * If there is any equality among the points in "bset", then we first 299 * project it out. Otherwise, we continue with solve_ilp above. 300 */ 301enum isl_lp_result isl_basic_set_solve_ilp(struct isl_basic_set *bset, int max, 302 isl_int *f, isl_int *opt, 303 struct isl_vec **sol_p) 304{ 305 unsigned dim; 306 enum isl_lp_result res; 307 308 if (!bset) 309 return isl_lp_error; 310 if (sol_p) 311 *sol_p = NULL; 312 313 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error); 314 315 if (isl_basic_set_plain_is_empty(bset)) 316 return isl_lp_empty; 317 318 if (bset->n_eq) 319 return solve_ilp_with_eq(bset, max, f, opt, sol_p); 320 321 dim = isl_basic_set_total_dim(bset); 322 323 if (max) 324 isl_seq_neg(f, f, 1 + dim); 325 326 res = solve_ilp(bset, f, opt, sol_p); 327 328 if (max) { 329 isl_seq_neg(f, f, 1 + dim); 330 isl_int_neg(*opt, *opt); 331 } 332 333 return res; 334error: 335 isl_basic_set_free(bset); 336 return isl_lp_error; 337} 338 339static enum isl_lp_result basic_set_opt(__isl_keep isl_basic_set *bset, int max, 340 __isl_keep isl_aff *obj, isl_int *opt) 341{ 342 enum isl_lp_result res; 343 344 if (!obj) 345 return isl_lp_error; 346 bset = isl_basic_set_copy(bset); 347 bset = isl_basic_set_underlying_set(bset); 348 res = isl_basic_set_solve_ilp(bset, max, obj->v->el + 1, opt, NULL); 349 isl_basic_set_free(bset); 350 return res; 351} 352 353static __isl_give isl_mat *extract_divs(__isl_keep isl_basic_set *bset) 354{ 355 int i; 356 isl_ctx *ctx = isl_basic_set_get_ctx(bset); 357 isl_mat *div; 358 359 div = isl_mat_alloc(ctx, bset->n_div, 360 1 + 1 + isl_basic_set_total_dim(bset)); 361 if (!div) 362 return NULL; 363 364 for (i = 0; i < bset->n_div; ++i) 365 isl_seq_cpy(div->row[i], bset->div[i], div->n_col); 366 367 return div; 368} 369 370enum isl_lp_result isl_basic_set_opt(__isl_keep isl_basic_set *bset, int max, 371 __isl_keep isl_aff *obj, isl_int *opt) 372{ 373 int *exp1 = NULL; 374 int *exp2 = NULL; 375 isl_ctx *ctx; 376 isl_mat *bset_div = NULL; 377 isl_mat *div = NULL; 378 enum isl_lp_result res; 379 int bset_n_div, obj_n_div; 380 381 if (!bset || !obj) 382 return isl_lp_error; 383 384 ctx = isl_aff_get_ctx(obj); 385 if (!isl_space_is_equal(bset->dim, obj->ls->dim)) 386 isl_die(ctx, isl_error_invalid, 387 "spaces don't match", return isl_lp_error); 388 if (!isl_int_is_one(obj->v->el[0])) 389 isl_die(ctx, isl_error_unsupported, 390 "expecting integer affine expression", 391 return isl_lp_error); 392 393 bset_n_div = isl_basic_set_dim(bset, isl_dim_div); 394 obj_n_div = isl_aff_dim(obj, isl_dim_div); 395 if (bset_n_div == 0 && obj_n_div == 0) 396 return basic_set_opt(bset, max, obj, opt); 397 398 bset = isl_basic_set_copy(bset); 399 obj = isl_aff_copy(obj); 400 401 bset_div = extract_divs(bset); 402 exp1 = isl_alloc_array(ctx, int, bset_n_div); 403 exp2 = isl_alloc_array(ctx, int, obj_n_div); 404 if (!bset_div || (bset_n_div && !exp1) || (obj_n_div && !exp2)) 405 goto error; 406 407 div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2); 408 409 bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1); 410 obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2); 411 412 res = basic_set_opt(bset, max, obj, opt); 413 414 isl_mat_free(bset_div); 415 isl_mat_free(div); 416 free(exp1); 417 free(exp2); 418 isl_basic_set_free(bset); 419 isl_aff_free(obj); 420 421 return res; 422error: 423 isl_mat_free(div); 424 isl_mat_free(bset_div); 425 free(exp1); 426 free(exp2); 427 isl_basic_set_free(bset); 428 isl_aff_free(obj); 429 return isl_lp_error; 430} 431 432/* Compute the minimum (maximum if max is set) of the integer affine 433 * expression obj over the points in set and put the result in *opt. 434 * 435 * The parameters are assumed to have been aligned. 436 */ 437static enum isl_lp_result isl_set_opt_aligned(__isl_keep isl_set *set, int max, 438 __isl_keep isl_aff *obj, isl_int *opt) 439{ 440 int i; 441 enum isl_lp_result res; 442 int empty = 1; 443 isl_int opt_i; 444 445 if (!set || !obj) 446 return isl_lp_error; 447 if (set->n == 0) 448 return isl_lp_empty; 449 450 res = isl_basic_set_opt(set->p[0], max, obj, opt); 451 if (res == isl_lp_error || res == isl_lp_unbounded) 452 return res; 453 if (set->n == 1) 454 return res; 455 if (res == isl_lp_ok) 456 empty = 0; 457 458 isl_int_init(opt_i); 459 for (i = 1; i < set->n; ++i) { 460 res = isl_basic_set_opt(set->p[i], max, obj, &opt_i); 461 if (res == isl_lp_error || res == isl_lp_unbounded) { 462 isl_int_clear(opt_i); 463 return res; 464 } 465 if (res == isl_lp_ok) 466 empty = 0; 467 if (max ? isl_int_gt(opt_i, *opt) : isl_int_lt(opt_i, *opt)) 468 isl_int_set(*opt, opt_i); 469 } 470 isl_int_clear(opt_i); 471 472 return empty ? isl_lp_empty : isl_lp_ok; 473} 474 475/* Compute the minimum (maximum if max is set) of the integer affine 476 * expression obj over the points in set and put the result in *opt. 477 */ 478enum isl_lp_result isl_set_opt(__isl_keep isl_set *set, int max, 479 __isl_keep isl_aff *obj, isl_int *opt) 480{ 481 enum isl_lp_result res; 482 483 if (!set || !obj) 484 return isl_lp_error; 485 486 if (isl_space_match(set->dim, isl_dim_param, 487 obj->ls->dim, isl_dim_param)) 488 return isl_set_opt_aligned(set, max, obj, opt); 489 490 set = isl_set_copy(set); 491 obj = isl_aff_copy(obj); 492 set = isl_set_align_params(set, isl_aff_get_domain_space(obj)); 493 obj = isl_aff_align_params(obj, isl_set_get_space(set)); 494 495 res = isl_set_opt_aligned(set, max, obj, opt); 496 497 isl_set_free(set); 498 isl_aff_free(obj); 499 500 return res; 501} 502 503enum isl_lp_result isl_basic_set_max(__isl_keep isl_basic_set *bset, 504 __isl_keep isl_aff *obj, isl_int *opt) 505{ 506 return isl_basic_set_opt(bset, 1, obj, opt); 507} 508 509enum isl_lp_result isl_set_max(__isl_keep isl_set *set, 510 __isl_keep isl_aff *obj, isl_int *opt) 511{ 512 return isl_set_opt(set, 1, obj, opt); 513} 514 515enum isl_lp_result isl_set_min(__isl_keep isl_set *set, 516 __isl_keep isl_aff *obj, isl_int *opt) 517{ 518 return isl_set_opt(set, 0, obj, opt); 519} 520 521/* Convert the result of a function that returns an isl_lp_result 522 * to an isl_val. The numerator of "v" is set to the optimal value 523 * if lp_res is isl_lp_ok. "max" is set if a maximum was computed. 524 * 525 * Return "v" with denominator set to 1 if lp_res is isl_lp_ok. 526 * Return NULL on error. 527 * Return a NaN if lp_res is isl_lp_empty. 528 * Return infinity or negative infinity if lp_res is isl_lp_unbounded, 529 * depending on "max". 530 */ 531static __isl_give isl_val *convert_lp_result(enum isl_lp_result lp_res, 532 __isl_take isl_val *v, int max) 533{ 534 isl_ctx *ctx; 535 536 if (lp_res == isl_lp_ok) { 537 isl_int_set_si(v->d, 1); 538 return isl_val_normalize(v); 539 } 540 ctx = isl_val_get_ctx(v); 541 isl_val_free(v); 542 if (lp_res == isl_lp_error) 543 return NULL; 544 if (lp_res == isl_lp_empty) 545 return isl_val_nan(ctx); 546 if (max) 547 return isl_val_infty(ctx); 548 else 549 return isl_val_neginfty(ctx); 550} 551 552/* Return the minimum (maximum if max is set) of the integer affine 553 * expression "obj" over the points in "bset". 554 * 555 * Return infinity or negative infinity if the optimal value is unbounded and 556 * NaN if "bset" is empty. 557 * 558 * Call isl_basic_set_opt and translate the results. 559 */ 560__isl_give isl_val *isl_basic_set_opt_val(__isl_keep isl_basic_set *bset, 561 int max, __isl_keep isl_aff *obj) 562{ 563 isl_ctx *ctx; 564 isl_val *res; 565 enum isl_lp_result lp_res; 566 567 if (!bset || !obj) 568 return NULL; 569 570 ctx = isl_aff_get_ctx(obj); 571 res = isl_val_alloc(ctx); 572 if (!res) 573 return NULL; 574 lp_res = isl_basic_set_opt(bset, max, obj, &res->n); 575 return convert_lp_result(lp_res, res, max); 576} 577 578/* Return the maximum of the integer affine 579 * expression "obj" over the points in "bset". 580 * 581 * Return infinity or negative infinity if the optimal value is unbounded and 582 * NaN if "bset" is empty. 583 */ 584__isl_give isl_val *isl_basic_set_max_val(__isl_keep isl_basic_set *bset, 585 __isl_keep isl_aff *obj) 586{ 587 return isl_basic_set_opt_val(bset, 1, obj); 588} 589 590/* Return the minimum (maximum if max is set) of the integer affine 591 * expression "obj" over the points in "set". 592 * 593 * Return infinity or negative infinity if the optimal value is unbounded and 594 * NaN if "bset" is empty. 595 * 596 * Call isl_set_opt and translate the results. 597 */ 598__isl_give isl_val *isl_set_opt_val(__isl_keep isl_set *set, int max, 599 __isl_keep isl_aff *obj) 600{ 601 isl_ctx *ctx; 602 isl_val *res; 603 enum isl_lp_result lp_res; 604 605 if (!set || !obj) 606 return NULL; 607 608 ctx = isl_aff_get_ctx(obj); 609 res = isl_val_alloc(ctx); 610 if (!res) 611 return NULL; 612 lp_res = isl_set_opt(set, max, obj, &res->n); 613 return convert_lp_result(lp_res, res, max); 614} 615 616/* Return the minimum of the integer affine 617 * expression "obj" over the points in "set". 618 * 619 * Return infinity or negative infinity if the optimal value is unbounded and 620 * NaN if "bset" is empty. 621 */ 622__isl_give isl_val *isl_set_min_val(__isl_keep isl_set *set, 623 __isl_keep isl_aff *obj) 624{ 625 return isl_set_opt_val(set, 0, obj); 626} 627 628/* Return the maximum of the integer affine 629 * expression "obj" over the points in "set". 630 * 631 * Return infinity or negative infinity if the optimal value is unbounded and 632 * NaN if "bset" is empty. 633 */ 634__isl_give isl_val *isl_set_max_val(__isl_keep isl_set *set, 635 __isl_keep isl_aff *obj) 636{ 637 return isl_set_opt_val(set, 1, obj); 638} 639