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
  17 /* Given a basic set "bset", construct a basic set U such that for
  18  * each element x in U, the whole unit box positioned at x is inside
  19  * the given basic set.
  20  * Note that U may not contain all points that satisfy this property.
  21  *
  22  * We simply add the sum of all negative coefficients to the constant
  23  * term.  This ensures that if x satisfies the resulting constraints,
  24  * then x plus any sum of unit vectors satisfies the original constraints.
  25  */
  26 static struct isl_basic_set *unit_box_base_points(struct isl_basic_set *bset)
  27 {
  28         int i, j, k;
  29         struct isl_basic_set *unit_box = NULL;
  30         unsigned total;
  31
  32         if (!bset)
  33                 goto error;
  34
  35         if (bset->n_eq != 0) {
  36                 unit_box = isl_basic_set_empty_like(bset);
  37                 isl_basic_set_free(bset);
  38                 return unit_box;
  39         }
  40
  41         total = isl_basic_set_total_dim(bset);
  42         unit_box = isl_basic_set_alloc_dim(isl_basic_set_get_dim(bset),
  43                                         0, 0, bset->n_ineq);
  44
  45         for (i = 0; i < bset->n_ineq; ++i) {
  46                 k = isl_basic_set_alloc_inequality(unit_box);
  47                 if (k < 0)
  48                         goto error;
  49                 isl_seq_cpy(unit_box->ineq[k], bset->ineq[i], 1 + total);
  50                 for (j = 0; j < total; ++j) {
  51                         if (isl_int_is_nonneg(unit_box->ineq[k][1 + j]))
  52                                 continue;
  53                         isl_int_add(unit_box->ineq[k][0],
  54                                 unit_box->ineq[k][0], unit_box->ineq[k][1 + j]);
  55                 }
  56         }
  57
  58         isl_basic_set_free(bset);
  59         return unit_box;
  60 error:
  61         isl_basic_set_free(bset);
  62         isl_basic_set_free(unit_box);
  63         return NULL;
  64 }
  65
  66 /* Find an integer point in "bset", preferably one that is
  67  * close to minimizing "f".
  68  *
  69  * We first check if we can easily put unit boxes inside bset.
  70  * If so, we take the best base point of any of the unit boxes we can find
  71  * and round it up to the nearest integer.
  72  * If not, we simply pick any integer point in "bset".
  73  */
  74 static struct isl_vec *initial_solution(struct isl_basic_set *bset, isl_int *f)
  75 {
  76         enum isl_lp_result res;
  77         struct isl_basic_set *unit_box;
  78         struct isl_vec *sol;
  79
  80         unit_box = unit_box_base_points(isl_basic_set_copy(bset));
  81
  82         res = isl_basic_set_solve_lp(unit_box, 0, f, bset->ctx->one,
  83                                         NULL, NULL, &sol);
  84         if (res == isl_lp_ok) {
  85                 isl_basic_set_free(unit_box);
  86                 return isl_vec_ceil(sol);
  87         }
  88
  89         isl_basic_set_free(unit_box);
  90
  91         return isl_basic_set_sample_vec(isl_basic_set_copy(bset));
  92 }
  93
  94 /* Restrict "bset" to those points with values for f in the interval [l, u].
  95  */
  96 static struct isl_basic_set *add_bounds(struct isl_basic_set *bset,
  97         isl_int *f, isl_int l, isl_int u)
  98 {
  99         int k;
 100         unsigned total;
 101
 102         total = isl_basic_set_total_dim(bset);
 103         bset = isl_basic_set_extend_constraints(bset, 0, 2);
 104
 105         k = isl_basic_set_alloc_inequality(bset);
 106         if (k < 0)
 107                 goto error;
 108         isl_seq_cpy(bset->ineq[k], f, 1 + total);
 109         isl_int_sub(bset->ineq[k][0], bset->ineq[k][0], l);
 110
 111         k = isl_basic_set_alloc_inequality(bset);
 112         if (k < 0)
 113                 goto error;
 114         isl_seq_neg(bset->ineq[k], f, 1 + total);
 115         isl_int_add(bset->ineq[k][0], bset->ineq[k][0], u);
 116
 117         return bset;
 118 error:
 119         isl_basic_set_free(bset);
 120         return NULL;
 121 }
 122
 123 /* Find an integer point in "bset" that minimizes f (in any) such that
 124  * the value of f lies inside the interval [l, u].
 125  * Return this integer point if it can be found.
 126  * Otherwise, return sol.
 127  *
 128  * We perform a number of steps until l > u.
 129  * In each step, we look for an integer point with value in either
 130  * the whole interval [l, u] or half of the interval [l, l+floor(u-l-1/2)].
 131  * The choice depends on whether we have found an integer point in the
 132  * previous step.  If so, we look for the next point in half of the remaining
 133  * interval.
 134  * If we find a point, the current solution is updated and u is set
 135  * to its value minus 1.
 136  * If no point can be found, we update l to the upper bound of the interval
 137  * we checked (u or l+floor(u-l-1/2)) plus 1.
 138  */
 139 static struct isl_vec *solve_ilp_search(struct isl_basic_set *bset,
 140         isl_int *f, isl_int *opt, struct isl_vec *sol, isl_int l, isl_int u)
 141 {
 142         isl_int tmp;
 143         int divide = 1;
 144
 145         isl_int_init(tmp);
 146
 147         while (isl_int_le(l, u)) {
 148                 struct isl_basic_set *slice;
 149                 struct isl_vec *sample;
 150
 151                 if (!divide)
 152                         isl_int_set(tmp, u);
 153                 else {
 154                         isl_int_sub(tmp, u, l);
 155                         isl_int_fdiv_q_ui(tmp, tmp, 2);
 156                         isl_int_add(tmp, tmp, l);
 157                 }
 158                 slice = add_bounds(isl_basic_set_copy(bset), f, l, tmp);
 159                 sample = isl_basic_set_sample_vec(slice);
 160                 if (!sample) {
 161                         isl_vec_free(sol);
 162                         sol = NULL;
 163                         break;
 164                 }
 165                 if (sample->size > 0) {
 166                         isl_vec_free(sol);
 167                         sol = sample;
 168                         isl_seq_inner_product(f, sol->el, sol->size, opt);
 169                         isl_int_sub_ui(u, *opt, 1);
 170                         divide = 1;
 171                 } else {
 172                         isl_vec_free(sample);
 173                         if (!divide)
 174                                 break;
 175                         isl_int_add_ui(l, tmp, 1);
 176                         divide = 0;
 177                 }
 178         }
 179
 180         isl_int_clear(tmp);
 181
 182         return sol;
 183 }
 184
 185 /* Find an integer point in "bset" that minimizes f (if any).
 186  * If sol_p is not NULL then the integer point is returned in *sol_p.
 187  * The optimal value of f is returned in *opt.
 188  *
 189  * The algorithm maintains a currently best solution and an interval [l, u]
 190  * of values of f for which integer solutions could potentially still be found.
 191  * The initial value of the best solution so far is any solution.
 192  * The initial value of l is minimal value of f over the rationals
 193  * (rounded up to the nearest integer).
 194  * The initial value of u is the value of f at the initial solution minus 1.
 195  *
 196  * We then call solve_ilp_search to perform a binary search on the interval.
 197  */
 198 static enum isl_lp_result solve_ilp(struct isl_basic_set *bset,
 199                                       isl_int *f, isl_int *opt,
 200                                       struct isl_vec **sol_p)
 201 {
 202         enum isl_lp_result res;
 203         isl_int l, u;
 204         struct isl_vec *sol;
 205
 206         res = isl_basic_set_solve_lp(bset, 0, f, bset->ctx->one,
 207                                         opt, NULL, &sol);
 208         if (res == isl_lp_ok && isl_int_is_one(sol->el[0])) {
 209                 if (sol_p)
 210                         *sol_p = sol;
 211                 else
 212                         isl_vec_free(sol);
 213                 return isl_lp_ok;
 214         }
 215         isl_vec_free(sol);
 216         if (res == isl_lp_error || res == isl_lp_empty)
 217                 return res;
 218
 219         sol = initial_solution(bset, f);
 220         if (!sol)
 221                 return isl_lp_error;
 222         if (sol->size == 0) {
 223                 isl_vec_free(sol);
 224                 return isl_lp_empty;
 225         }
 226         if (res == isl_lp_unbounded) {
 227                 isl_vec_free(sol);
 228                 return isl_lp_unbounded;
 229         }
 230
 231         isl_int_init(l);
 232         isl_int_init(u);
 233
 234         isl_int_set(l, *opt);
 235
 236         isl_seq_inner_product(f, sol->el, sol->size, opt);
 237         isl_int_sub_ui(u, *opt, 1);
 238
 239         sol = solve_ilp_search(bset, f, opt, sol, l, u);
 240         if (!sol)
 241                 res = isl_lp_error;
 242
 243         isl_int_clear(l);
 244         isl_int_clear(u);
 245
 246         if (sol_p)
 247                 *sol_p = sol;
 248         else
 249                 isl_vec_free(sol);
 250
 251         return res;
 252 }
 253
 254 static enum isl_lp_result solve_ilp_with_eq(struct isl_basic_set *bset, int max,
 255                                       isl_int *f, isl_int *opt,
 256                                       struct isl_vec **sol_p)
 257 {
 258         unsigned dim;
 259         enum isl_lp_result res;
 260         struct isl_mat *T = NULL;
 261         struct isl_vec *v;
 262
 263         bset = isl_basic_set_copy(bset);
 264         dim = isl_basic_set_total_dim(bset);
 265         v = isl_vec_alloc(bset->ctx, 1 + dim);
 266         if (!v)
 267                 goto error;
 268         isl_seq_cpy(v->el, f, 1 + dim);
 269         bset = isl_basic_set_remove_equalities(bset, &T, NULL);
 270         v = isl_vec_mat_product(v, isl_mat_copy(T));
 271         if (!v)
 272                 goto error;
 273         res = isl_basic_set_solve_ilp(bset, max, v->el, opt, sol_p);
 274         isl_vec_free(v);
 275         if (res == isl_lp_ok && sol_p) {
 276                 *sol_p = isl_mat_vec_product(T, *sol_p);
 277                 if (!*sol_p)
 278                         res = isl_lp_error;
 279         } else
 280                 isl_mat_free(T);
 281         isl_basic_set_free(bset);
 282         return res;
 283 error:
 284         isl_mat_free(T);
 285         isl_basic_set_free(bset);
 286         return isl_lp_error;
 287 }
 288
 289 /* Find an integer point in "bset" that minimizes (or maximizes if max is set)
 290  * f (if any).
 291  * If sol_p is not NULL then the integer point is returned in *sol_p.
 292  * The optimal value of f is returned in *opt.
 293  *
 294  * If there is any equality among the points in "bset", then we first
 295  * project it out.  Otherwise, we continue with solve_ilp above.
 296  */
 297 enum isl_lp_result isl_basic_set_solve_ilp(struct isl_basic_set *bset, int max,
 298                                       isl_int *f, isl_int *opt,
 299                                       struct isl_vec **sol_p)
 300 {
 301         unsigned dim;
 302         enum isl_lp_result res;
 303
 304         if (!bset)
 305                 return isl_lp_error;
 306         if (sol_p)
 307                 *sol_p = NULL;
 308
 309         isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
 310
 311         if (isl_basic_set_fast_is_empty(bset))
 312                 return isl_lp_empty;
 313
 314         if (bset->n_eq)
 315                 return solve_ilp_with_eq(bset, max, f, opt, sol_p);
 316
 317         dim = isl_basic_set_total_dim(bset);
 318
 319         if (max)
 320                 isl_seq_neg(f, f, 1 + dim);
 321
 322         res = solve_ilp(bset, f, opt, sol_p);
 323
 324         if (max) {
 325                 isl_seq_neg(f, f, 1 + dim);
 326                 isl_int_neg(*opt, *opt);
 327         }
 328
 329         return res;
 330 error:
 331         isl_basic_set_free(bset);
 332         return isl_lp_error;
 333 }