isl_ilp.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 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
  20 /* Given a basic set "bset", construct a basic set U such that for
  21  * each element x in U, the whole unit box positioned at x is inside
  22  * the given basic set.
  23  * Note that U may not contain all points that satisfy this property.
  24  *
  25  * We simply add the sum of all negative coefficients to the constant
  26  * term.  This ensures that if x satisfies the resulting constraints,
  27  * then x plus any sum of unit vectors satisfies the original constraints.
  28  */
  29 static struct isl_basic_set *unit_box_base_points(struct isl_basic_set *bset)
  30 {
  31         int i, j, k;
  32         struct isl_basic_set *unit_box = NULL;
  33         unsigned total;
  34
  35         if (!bset)
  36                 goto error;
  37
  38         if (bset->n_eq != 0) {
  39                 unit_box = isl_basic_set_empty_like(bset);
  40                 isl_basic_set_free(bset);
  41                 return unit_box;
  42         }
  43
  44         total = isl_basic_set_total_dim(bset);
  45         unit_box = isl_basic_set_alloc_dim(isl_basic_set_get_dim(bset),
  46                                         0, 0, bset->n_ineq);
  47
  48         for (i = 0; i < bset->n_ineq; ++i) {
  49                 k = isl_basic_set_alloc_inequality(unit_box);
  50                 if (k < 0)
  51                         goto error;
  52                 isl_seq_cpy(unit_box->ineq[k], bset->ineq[i], 1 + total);
  53                 for (j = 0; j < total; ++j) {
  54                         if (isl_int_is_nonneg(unit_box->ineq[k][1 + j]))
  55                                 continue;
  56                         isl_int_add(unit_box->ineq[k][0],
  57                                 unit_box->ineq[k][0], unit_box->ineq[k][1 + j]);
  58                 }
  59         }
  60
  61         isl_basic_set_free(bset);
  62         return unit_box;
  63 error:
  64         isl_basic_set_free(bset);
  65         isl_basic_set_free(unit_box);
  66         return NULL;
  67 }
  68
  69 /* Find an integer point in "bset", preferably one that is
  70  * close to minimizing "f".
  71  *
  72  * We first check if we can easily put unit boxes inside bset.
  73  * If so, we take the best base point of any of the unit boxes we can find
  74  * and round it up to the nearest integer.
  75  * If not, we simply pick any integer point in "bset".
  76  */
  77 static struct isl_vec *initial_solution(struct isl_basic_set *bset, isl_int *f)
  78 {
  79         enum isl_lp_result res;
  80         struct isl_basic_set *unit_box;
  81         struct isl_vec *sol;
  82
  83         unit_box = unit_box_base_points(isl_basic_set_copy(bset));
  84
  85         res = isl_basic_set_solve_lp(unit_box, 0, f, bset->ctx->one,
  86                                         NULL, NULL, &sol);
  87         if (res == isl_lp_ok) {
  88                 isl_basic_set_free(unit_box);
  89                 return isl_vec_ceil(sol);
  90         }
  91
  92         isl_basic_set_free(unit_box);
  93
  94         return isl_basic_set_sample_vec(isl_basic_set_copy(bset));
  95 }
  96
  97 /* Restrict "bset" to those points with values for f in the interval [l, u].
  98  */
  99 static struct isl_basic_set *add_bounds(struct isl_basic_set *bset,
 100         isl_int *f, isl_int l, isl_int u)
 101 {
 102         int k;
 103         unsigned total;
 104
 105         total = isl_basic_set_total_dim(bset);
 106         bset = isl_basic_set_extend_constraints(bset, 0, 2);
 107
 108         k = isl_basic_set_alloc_inequality(bset);
 109         if (k < 0)
 110                 goto error;
 111         isl_seq_cpy(bset->ineq[k], f, 1 + total);
 112         isl_int_sub(bset->ineq[k][0], bset->ineq[k][0], l);
 113
 114         k = isl_basic_set_alloc_inequality(bset);
 115         if (k < 0)
 116                 goto error;
 117         isl_seq_neg(bset->ineq[k], f, 1 + total);
 118         isl_int_add(bset->ineq[k][0], bset->ineq[k][0], u);
 119
 120         return bset;
 121 error:
 122         isl_basic_set_free(bset);
 123         return NULL;
 124 }
 125
 126 /* Find an integer point in "bset" that minimizes f (in any) such that
 127  * the value of f lies inside the interval [l, u].
 128  * Return this integer point if it can be found.
 129  * Otherwise, return sol.
 130  *
 131  * We perform a number of steps until l > u.
 132  * In each step, we look for an integer point with value in either
 133  * the whole interval [l, u] or half of the interval [l, l+floor(u-l-1/2)].
 134  * The choice depends on whether we have found an integer point in the
 135  * previous step.  If so, we look for the next point in half of the remaining
 136  * interval.
 137  * If we find a point, the current solution is updated and u is set
 138  * to its value minus 1.
 139  * If no point can be found, we update l to the upper bound of the interval
 140  * we checked (u or l+floor(u-l-1/2)) plus 1.
 141  */
 142 static struct isl_vec *solve_ilp_search(struct isl_basic_set *bset,
 143         isl_int *f, isl_int *opt, struct isl_vec *sol, isl_int l, isl_int u)
 144 {
 145         isl_int tmp;
 146         int divide = 1;
 147
 148         isl_int_init(tmp);
 149
 150         while (isl_int_le(l, u)) {
 151                 struct isl_basic_set *slice;
 152                 struct isl_vec *sample;
 153
 154                 if (!divide)
 155                         isl_int_set(tmp, u);
 156                 else {
 157                         isl_int_sub(tmp, u, l);
 158                         isl_int_fdiv_q_ui(tmp, tmp, 2);
 159                         isl_int_add(tmp, tmp, l);
 160                 }
 161                 slice = add_bounds(isl_basic_set_copy(bset), f, l, tmp);
 162                 sample = isl_basic_set_sample_vec(slice);
 163                 if (!sample) {
 164                         isl_vec_free(sol);
 165                         sol = NULL;
 166                         break;
 167                 }
 168                 if (sample->size > 0) {
 169                         isl_vec_free(sol);
 170                         sol = sample;
 171                         isl_seq_inner_product(f, sol->el, sol->size, opt);
 172                         isl_int_sub_ui(u, *opt, 1);
 173                         divide = 1;
 174                 } else {
 175                         isl_vec_free(sample);
 176                         if (!divide)
 177                                 break;
 178                         isl_int_add_ui(l, tmp, 1);
 179                         divide = 0;
 180                 }
 181         }
 182
 183         isl_int_clear(tmp);
 184
 185         return sol;
 186 }
 187
 188 /* Find an integer point in "bset" that minimizes f (if any).
 189  * If sol_p is not NULL then the integer point is returned in *sol_p.
 190  * The optimal value of f is returned in *opt.
 191  *
 192  * The algorithm maintains a currently best solution and an interval [l, u]
 193  * of values of f for which integer solutions could potentially still be found.
 194  * The initial value of the best solution so far is any solution.
 195  * The initial value of l is minimal value of f over the rationals
 196  * (rounded up to the nearest integer).
 197  * The initial value of u is the value of f at the initial solution minus 1.
 198  *
 199  * We then call solve_ilp_search to perform a binary search on the interval.
 200  */
 201 static enum isl_lp_result solve_ilp(struct isl_basic_set *bset,
 202                                       isl_int *f, isl_int *opt,
 203                                       struct isl_vec **sol_p)
 204 {
 205         enum isl_lp_result res;
 206         isl_int l, u;
 207         struct isl_vec *sol;
 208
 209         res = isl_basic_set_solve_lp(bset, 0, f, bset->ctx->one,
 210                                         opt, NULL, &sol);
 211         if (res == isl_lp_ok && isl_int_is_one(sol->el[0])) {
 212                 if (sol_p)
 213                         *sol_p = sol;
 214                 else
 215                         isl_vec_free(sol);
 216                 return isl_lp_ok;
 217         }
 218         isl_vec_free(sol);
 219         if (res == isl_lp_error || res == isl_lp_empty)
 220                 return res;
 221
 222         sol = initial_solution(bset, f);
 223         if (!sol)
 224                 return isl_lp_error;
 225         if (sol->size == 0) {
 226                 isl_vec_free(sol);
 227                 return isl_lp_empty;
 228         }
 229         if (res == isl_lp_unbounded) {
 230                 isl_vec_free(sol);
 231                 return isl_lp_unbounded;
 232         }
 233
 234         isl_int_init(l);
 235         isl_int_init(u);
 236
 237         isl_int_set(l, *opt);
 238
 239         isl_seq_inner_product(f, sol->el, sol->size, opt);
 240         isl_int_sub_ui(u, *opt, 1);
 241
 242         sol = solve_ilp_search(bset, f, opt, sol, l, u);
 243         if (!sol)
 244                 res = isl_lp_error;
 245
 246         isl_int_clear(l);
 247         isl_int_clear(u);
 248
 249         if (sol_p)
 250                 *sol_p = sol;
 251         else
 252                 isl_vec_free(sol);
 253
 254         return res;
 255 }
 256
 257 static enum isl_lp_result solve_ilp_with_eq(struct isl_basic_set *bset, int max,
 258                                       isl_int *f, isl_int *opt,
 259                                       struct isl_vec **sol_p)
 260 {
 261         unsigned dim;
 262         enum isl_lp_result res;
 263         struct isl_mat *T = NULL;
 264         struct isl_vec *v;
 265
 266         bset = isl_basic_set_copy(bset);
 267         dim = isl_basic_set_total_dim(bset);
 268         v = isl_vec_alloc(bset->ctx, 1 + dim);
 269         if (!v)
 270                 goto error;
 271         isl_seq_cpy(v->el, f, 1 + dim);
 272         bset = isl_basic_set_remove_equalities(bset, &T, NULL);
 273         v = isl_vec_mat_product(v, isl_mat_copy(T));
 274         if (!v)
 275                 goto error;
 276         res = isl_basic_set_solve_ilp(bset, max, v->el, opt, sol_p);
 277         isl_vec_free(v);
 278         if (res == isl_lp_ok && sol_p) {
 279                 *sol_p = isl_mat_vec_product(T, *sol_p);
 280                 if (!*sol_p)
 281                         res = isl_lp_error;
 282         } else
 283                 isl_mat_free(T);
 284         isl_basic_set_free(bset);
 285         return res;
 286 error:
 287         isl_mat_free(T);
 288         isl_basic_set_free(bset);
 289         return isl_lp_error;
 290 }
 291
 292 /* Find an integer point in "bset" that minimizes (or maximizes if max is set)
 293  * f (if any).
 294  * If sol_p is not NULL then the integer point is returned in *sol_p.
 295  * The optimal value of f is returned in *opt.
 296  *
 297  * If there is any equality among the points in "bset", then we first
 298  * project it out.  Otherwise, we continue with solve_ilp above.
 299  */
 300 enum isl_lp_result isl_basic_set_solve_ilp(struct isl_basic_set *bset, int max,
 301                                       isl_int *f, isl_int *opt,
 302                                       struct isl_vec **sol_p)
 303 {
 304         unsigned dim;
 305         enum isl_lp_result res;
 306
 307         if (!bset)
 308                 return isl_lp_error;
 309         if (sol_p)
 310                 *sol_p = NULL;
 311
 312         isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
 313
 314         if (isl_basic_set_plain_is_empty(bset))
 315                 return isl_lp_empty;
 316
 317         if (bset->n_eq)
 318                 return solve_ilp_with_eq(bset, max, f, opt, sol_p);
 319
 320         dim = isl_basic_set_total_dim(bset);
 321
 322         if (max)
 323                 isl_seq_neg(f, f, 1 + dim);
 324
 325         res = solve_ilp(bset, f, opt, sol_p);
 326
 327         if (max) {
 328                 isl_seq_neg(f, f, 1 + dim);
 329                 isl_int_neg(*opt, *opt);
 330         }
 331
 332         return res;
 333 error:
 334         isl_basic_set_free(bset);
 335         return isl_lp_error;
 336 }
 337
 338 static enum isl_lp_result basic_set_opt(__isl_keep isl_basic_set *bset, int max,
 339         __isl_keep isl_aff *obj, isl_int *opt)
 340 {
 341         enum isl_lp_result res;
 342
 343         if (!obj)
 344                 return isl_lp_error;
 345         bset = isl_basic_set_copy(bset);
 346         bset = isl_basic_set_underlying_set(bset);
 347         res = isl_basic_set_solve_ilp(bset, max, obj->v->el + 1, opt, NULL);
 348         isl_basic_set_free(bset);
 349         return res;
 350 }
 351
 352 static __isl_give isl_mat *extract_divs(__isl_keep isl_basic_set *bset)
 353 {
 354         int i;
 355         isl_ctx *ctx = isl_basic_set_get_ctx(bset);
 356         isl_mat *div;
 357
 358         div = isl_mat_alloc(ctx, bset->n_div,
 359                             1 + 1 + isl_basic_set_total_dim(bset));
 360         if (!div)
 361                 return NULL;
 362
 363         for (i = 0; i < bset->n_div; ++i)
 364                 isl_seq_cpy(div->row[i], bset->div[i], div->n_col);
 365
 366         return div;
 367 }
 368
 369 enum isl_lp_result isl_basic_set_opt(__isl_keep isl_basic_set *bset, int max,
 370         __isl_keep isl_aff *obj, isl_int *opt)
 371 {
 372         int *exp1 = NULL;
 373         int *exp2 = NULL;
 374         isl_ctx *ctx;
 375         isl_mat *bset_div = NULL;
 376         isl_mat *div = NULL;
 377         enum isl_lp_result res;
 378
 379         if (!bset || !obj)
 380                 return isl_lp_error;
 381
 382         ctx = isl_aff_get_ctx(obj);
 383         if (!isl_dim_equal(bset->dim, obj->ls->dim))
 384                 isl_die(ctx, isl_error_invalid,
 385                         "spaces don't match", return isl_lp_error);
 386         if (!isl_int_is_one(obj->v->el[0]))
 387                 isl_die(ctx, isl_error_unsupported,
 388                         "expecting integer affine expression",
 389                         return isl_lp_error);
 390
 391         if (bset->n_div == 0 && obj->ls->div->n_row == 0)
 392                 return basic_set_opt(bset, max, obj, opt);
 393
 394         bset = isl_basic_set_copy(bset);
 395         obj = isl_aff_copy(obj);
 396
 397         bset_div = extract_divs(bset);
 398         exp1 = isl_alloc_array(ctx, int, bset_div->n_row);
 399         exp2 = isl_alloc_array(ctx, int, obj->ls->div->n_row);
 400         if (!bset_div || !exp1 || !exp2)
 401                 goto error;
 402
 403         div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2);
 404
 405         bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1);
 406         obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2);
 407
 408         res = basic_set_opt(bset, max, obj, opt);
 409
 410         isl_mat_free(bset_div);
 411         isl_mat_free(div);
 412         free(exp1);
 413         free(exp2);
 414         isl_basic_set_free(bset);
 415         isl_aff_free(obj);
 416
 417         return res;
 418 error:
 419         isl_mat_free(div);
 420         isl_mat_free(bset_div);
 421         free(exp1);
 422         free(exp2);
 423         isl_basic_set_free(bset);
 424         isl_aff_free(obj);
 425         return isl_lp_error;
 426 }
 427
 428 /* Compute the minimum (maximum if max is set) of the integer affine
 429  * expression obj over the points in set and put the result in *opt.
 430  */
 431 enum isl_lp_result isl_set_opt(__isl_keep isl_set *set, int max,
 432         __isl_keep isl_aff *obj, isl_int *opt)
 433 {
 434         int i;
 435         enum isl_lp_result res;
 436         int empty = 1;
 437         isl_int opt_i;
 438
 439         if (!set || !obj)
 440                 return isl_lp_error;
 441         if (set->n == 0)
 442                 return isl_lp_empty;
 443
 444         res = isl_basic_set_opt(set->p[0], max, obj, opt);
 445         if (res == isl_lp_error || res == isl_lp_unbounded)
 446                 return res;
 447         if (set->n == 1)
 448                 return res;
 449         if (res == isl_lp_ok)
 450                 empty = 0;
 451
 452         isl_int_init(opt_i);
 453         for (i = 1; i < set->n; ++i) {
 454                 res = isl_basic_set_opt(set->p[i], max, obj, &opt_i);
 455                 if (res == isl_lp_error || res == isl_lp_unbounded) {
 456                         isl_int_clear(opt_i);
 457                         return res;
 458                 }
 459                 if (res == isl_lp_ok)
 460                         empty = 0;
 461                 if (isl_int_gt(opt_i, *opt))
 462                         isl_int_set(*opt, opt_i);
 463         }
 464         isl_int_clear(opt_i);
 465
 466         return empty ? isl_lp_empty : isl_lp_ok;
 467 }
 468
 469 enum isl_lp_result isl_basic_set_max(__isl_keep isl_basic_set *bset,
 470         __isl_keep isl_aff *obj, isl_int *opt)
 471 {
 472         return isl_basic_set_opt(bset, 1, obj, opt);
 473 }
 474
 475 enum isl_lp_result isl_set_max(__isl_keep isl_set *set,
 476         __isl_keep isl_aff *obj, isl_int *opt)
 477 {
 478         return isl_set_opt(set, 1, obj, opt);
 479 }