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