isl_ilp.c

   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  */
  30 static 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;
  64 error:
  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  */
  78 static 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  */
 100 static 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;
 122 error:
 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  */
 143 static 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  */
 202 static 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
 258 static 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;
 287 error:
 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  */
 301 enum 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;
 334 error:
 335         isl_basic_set_free(bset);
 336         return isl_lp_error;
 337 }
 338
 339 static 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
 353 static __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
 370 enum 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;
 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         if (bset_n_div == 0 && obj->ls->div->n_row == 0)
 395                 return basic_set_opt(bset, max, obj, opt);
 396
 397         bset = isl_basic_set_copy(bset);
 398         obj = isl_aff_copy(obj);
 399
 400         bset_div = extract_divs(bset);
 401         exp1 = isl_alloc_array(ctx, int, bset_n_div);
 402         exp2 = isl_alloc_array(ctx, int, obj->ls->div->n_row);
 403         if (!bset_div || !exp1 || !exp2)
 404                 goto error;
 405
 406         div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2);
 407
 408         bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1);
 409         obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2);
 410
 411         res = basic_set_opt(bset, max, obj, opt);
 412
 413         isl_mat_free(bset_div);
 414         isl_mat_free(div);
 415         free(exp1);
 416         free(exp2);
 417         isl_basic_set_free(bset);
 418         isl_aff_free(obj);
 419
 420         return res;
 421 error:
 422         isl_mat_free(div);
 423         isl_mat_free(bset_div);
 424         free(exp1);
 425         free(exp2);
 426         isl_basic_set_free(bset);
 427         isl_aff_free(obj);
 428         return isl_lp_error;
 429 }
 430
 431 /* Compute the minimum (maximum if max is set) of the integer affine
 432  * expression obj over the points in set and put the result in *opt.
 433  *
 434  * The parameters are assumed to have been aligned.
 435  */
 436 static enum isl_lp_result isl_set_opt_aligned(__isl_keep isl_set *set, int max,
 437         __isl_keep isl_aff *obj, isl_int *opt)
 438 {
 439         int i;
 440         enum isl_lp_result res;
 441         int empty = 1;
 442         isl_int opt_i;
 443
 444         if (!set || !obj)
 445                 return isl_lp_error;
 446         if (set->n == 0)
 447                 return isl_lp_empty;
 448
 449         res = isl_basic_set_opt(set->p[0], max, obj, opt);
 450         if (res == isl_lp_error || res == isl_lp_unbounded)
 451                 return res;
 452         if (set->n == 1)
 453                 return res;
 454         if (res == isl_lp_ok)
 455                 empty = 0;
 456
 457         isl_int_init(opt_i);
 458         for (i = 1; i < set->n; ++i) {
 459                 res = isl_basic_set_opt(set->p[i], max, obj, &opt_i);
 460                 if (res == isl_lp_error || res == isl_lp_unbounded) {
 461                         isl_int_clear(opt_i);
 462                         return res;
 463                 }
 464                 if (res == isl_lp_ok)
 465                         empty = 0;
 466                 if (isl_int_gt(opt_i, *opt))
 467                         isl_int_set(*opt, opt_i);
 468         }
 469         isl_int_clear(opt_i);
 470
 471         return empty ? isl_lp_empty : isl_lp_ok;
 472 }
 473
 474 /* Compute the minimum (maximum if max is set) of the integer affine
 475  * expression obj over the points in set and put the result in *opt.
 476  */
 477 enum isl_lp_result isl_set_opt(__isl_keep isl_set *set, int max,
 478         __isl_keep isl_aff *obj, isl_int *opt)
 479 {
 480         enum isl_lp_result res;
 481
 482         if (!set || !obj)
 483                 return isl_lp_error;
 484
 485         if (isl_space_match(set->dim, isl_dim_param,
 486                             obj->ls->dim, isl_dim_param))
 487                 return isl_set_opt_aligned(set, max, obj, opt);
 488
 489         set = isl_set_copy(set);
 490         obj = isl_aff_copy(obj);
 491         set = isl_set_align_params(set, isl_aff_get_domain_space(obj));
 492         obj = isl_aff_align_params(obj, isl_set_get_space(set));
 493
 494         res = isl_set_opt_aligned(set, max, obj, opt);
 495
 496         isl_set_free(set);
 497         isl_aff_free(obj);
 498
 499         return res;
 500 }
 501
 502 enum isl_lp_result isl_basic_set_max(__isl_keep isl_basic_set *bset,
 503         __isl_keep isl_aff *obj, isl_int *opt)
 504 {
 505         return isl_basic_set_opt(bset, 1, obj, opt);
 506 }
 507
 508 enum isl_lp_result isl_set_max(__isl_keep isl_set *set,
 509         __isl_keep isl_aff *obj, isl_int *opt)
 510 {
 511         return isl_set_opt(set, 1, obj, opt);
 512 }
 513
 514 enum isl_lp_result isl_set_min(__isl_keep isl_set *set,
 515         __isl_keep isl_aff *obj, isl_int *opt)
 516 {
 517         return isl_set_opt(set, 0, obj, opt);
 518 }
 519
 520 /* Convert the result of a function that returns an isl_lp_result
 521  * to an isl_val.  The numerator of "v" is set to the optimal value
 522  * if lp_res is isl_lp_ok.  "max" is set if a maximum was computed.
 523  *
 524  * Return "v" with denominator set to 1 if lp_res is isl_lp_ok.
 525  * Return NULL on error.
 526  * Return a NaN if lp_res is isl_lp_empty.
 527  * Return infinity or negative infinity if lp_res is isl_lp_unbounded,
 528  * depending on "max".
 529  */
 530 static __isl_give isl_val *convert_lp_result(enum isl_lp_result lp_res,
 531         __isl_take isl_val *v, int max)
 532 {
 533         isl_ctx *ctx;
 534
 535         if (lp_res == isl_lp_ok) {
 536                 isl_int_set_si(v->d, 1);
 537                 return isl_val_normalize(v);
 538         }
 539         ctx = isl_val_get_ctx(v);
 540         isl_val_free(v);
 541         if (lp_res == isl_lp_error)
 542                 return NULL;
 543         if (lp_res == isl_lp_empty)
 544                 return isl_val_nan(ctx);
 545         if (max)
 546                 return isl_val_infty(ctx);
 547         else
 548                 return isl_val_neginfty(ctx);
 549 }
 550
 551 /* Return the minimum (maximum if max is set) of the integer affine
 552  * expression "obj" over the points in "bset".
 553  *
 554  * Return infinity or negative infinity if the optimal value is unbounded and
 555  * NaN if "bset" is empty.
 556  *
 557  * Call isl_basic_set_opt and translate the results.
 558  */
 559 __isl_give isl_val *isl_basic_set_opt_val(__isl_keep isl_basic_set *bset,
 560         int max, __isl_keep isl_aff *obj)
 561 {
 562         isl_ctx *ctx;
 563         isl_val *res;
 564         enum isl_lp_result lp_res;
 565
 566         if (!bset || !obj)
 567                 return NULL;
 568
 569         ctx = isl_aff_get_ctx(obj);
 570         res = isl_val_alloc(ctx);
 571         if (!res)
 572                 return NULL;
 573         lp_res = isl_basic_set_opt(bset, max, obj, &res->n);
 574         return convert_lp_result(lp_res, res, max);
 575 }
 576
 577 /* Return the maximum of the integer affine
 578  * expression "obj" over the points in "bset".
 579  *
 580  * Return infinity or negative infinity if the optimal value is unbounded and
 581  * NaN if "bset" is empty.
 582  */
 583 __isl_give isl_val *isl_basic_set_max_val(__isl_keep isl_basic_set *bset,
 584         __isl_keep isl_aff *obj)
 585 {
 586         return isl_basic_set_opt_val(bset, 1, obj);
 587 }
 588
 589 /* Return the minimum (maximum if max is set) of the integer affine
 590  * expression "obj" over the points in "set".
 591  *
 592  * Return infinity or negative infinity if the optimal value is unbounded and
 593  * NaN if "bset" is empty.
 594  *
 595  * Call isl_set_opt and translate the results.
 596  */
 597 __isl_give isl_val *isl_set_opt_val(__isl_keep isl_set *set, int max,
 598         __isl_keep isl_aff *obj)
 599 {
 600         isl_ctx *ctx;
 601         isl_val *res;
 602         enum isl_lp_result lp_res;
 603
 604         if (!set || !obj)
 605                 return NULL;
 606
 607         ctx = isl_aff_get_ctx(obj);
 608         res = isl_val_alloc(ctx);
 609         if (!res)
 610                 return NULL;
 611         lp_res = isl_set_opt(set, max, obj, &res->n);
 612         return convert_lp_result(lp_res, res, max);
 613 }
 614
 615 /* Return the minimum of the integer affine
 616  * expression "obj" over the points in "set".
 617  *
 618  * Return infinity or negative infinity if the optimal value is unbounded and
 619  * NaN if "bset" is empty.
 620  */
 621 __isl_give isl_val *isl_set_min_val(__isl_keep isl_set *set,
 622         __isl_keep isl_aff *obj)
 623 {
 624         return isl_set_opt_val(set, 0, obj);
 625 }
 626
 627 /* Return the maximum of the integer affine
 628  * expression "obj" over the points in "set".
 629  *
 630  * Return infinity or negative infinity if the optimal value is unbounded and
 631  * NaN if "bset" is empty.
 632  */
 633 __isl_give isl_val *isl_set_max_val(__isl_keep isl_set *set,
 634         __isl_keep isl_aff *obj)
 635 {
 636         return isl_set_opt_val(set, 1, obj);
 637 }