isl_lp.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/lp.h>
  13 #include "isl_lp_piplib.h"
  14 #include <isl_seq.h>
  15 #include "isl_tab.h"
  16 #include <isl_options_private.h>
  17 #include <isl_local_space_private.h>
  18 #include <isl_aff_private.h>
  19 #include <isl_mat_private.h>
  20 #include <isl_val_private.h>
  21 #include <isl_vec_private.h>
  22
  23 enum isl_lp_result isl_tab_solve_lp(struct isl_basic_map *bmap, int maximize,
  24                                       isl_int *f, isl_int denom, isl_int *opt,
  25                                       isl_int *opt_denom,
  26                                       struct isl_vec **sol)
  27 {
  28         struct isl_tab *tab;
  29         enum isl_lp_result res;
  30         unsigned dim = isl_basic_map_total_dim(bmap);
  31
  32         if (maximize)
  33                 isl_seq_neg(f, f, 1 + dim);
  34
  35         bmap = isl_basic_map_gauss(bmap, NULL);
  36         tab = isl_tab_from_basic_map(bmap, 0);
  37         res = isl_tab_min(tab, f, denom, opt, opt_denom, 0);
  38         if (res == isl_lp_ok && sol) {
  39                 *sol = isl_tab_get_sample_value(tab);
  40                 if (!*sol)
  41                         res = isl_lp_error;
  42         }
  43         isl_tab_free(tab);
  44
  45         if (maximize)
  46                 isl_seq_neg(f, f, 1 + dim);
  47         if (maximize && opt)
  48                 isl_int_neg(*opt, *opt);
  49
  50         return res;
  51 }
  52
  53 /* Given a basic map "bmap" and an affine combination of the variables "f"
  54  * with denominator "denom", set *opt / *opt_denom to the minimal
  55  * (or maximal if "maximize" is true) value attained by f/d over "bmap",
  56  * assuming the basic map is not empty and the expression cannot attain
  57  * arbitrarily small (or large) values.
  58  * If opt_denom is NULL, then *opt is rounded up (or down)
  59  * to the nearest integer.
  60  * The return value reflects the nature of the result (empty, unbounded,
  61  * minmimal or maximal value returned in *opt).
  62  */
  63 enum isl_lp_result isl_basic_map_solve_lp(struct isl_basic_map *bmap, int max,
  64                                       isl_int *f, isl_int d, isl_int *opt,
  65                                       isl_int *opt_denom,
  66                                       struct isl_vec **sol)
  67 {
  68         if (sol)
  69                 *sol = NULL;
  70
  71         if (!bmap)
  72                 return isl_lp_error;
  73
  74         switch (bmap->ctx->opt->lp_solver) {
  75         case ISL_LP_PIP:
  76                 return isl_pip_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
  77         case ISL_LP_TAB:
  78                 return isl_tab_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
  79         default:
  80                 return isl_lp_error;
  81         }
  82 }
  83
  84 enum isl_lp_result isl_basic_set_solve_lp(struct isl_basic_set *bset, int max,
  85                                       isl_int *f, isl_int d, isl_int *opt,
  86                                       isl_int *opt_denom,
  87                                       struct isl_vec **sol)
  88 {
  89         return isl_basic_map_solve_lp((struct isl_basic_map *)bset, max,
  90                                         f, d, opt, opt_denom, sol);
  91 }
  92
  93 enum isl_lp_result isl_map_solve_lp(__isl_keep isl_map *map, int max,
  94                                       isl_int *f, isl_int d, isl_int *opt,
  95                                       isl_int *opt_denom,
  96                                       struct isl_vec **sol)
  97 {
  98         int i;
  99         isl_int o;
 100         isl_int t;
 101         isl_int opt_i;
 102         isl_int opt_denom_i;
 103         enum isl_lp_result res;
 104         int max_div;
 105         isl_vec *v = NULL;
 106
 107         if (!map)
 108                 return isl_lp_error;
 109         if (map->n == 0)
 110                 return isl_lp_empty;
 111
 112         max_div = 0;
 113         for (i = 0; i < map->n; ++i)
 114                 if (map->p[i]->n_div > max_div)
 115                         max_div = map->p[i]->n_div;
 116         if (max_div > 0) {
 117                 unsigned total = isl_space_dim(map->dim, isl_dim_all);
 118                 v = isl_vec_alloc(map->ctx, 1 + total + max_div);
 119                 if (!v)
 120                         return isl_lp_error;
 121                 isl_seq_cpy(v->el, f, 1 + total);
 122                 isl_seq_clr(v->el + 1 + total, max_div);
 123                 f = v->el;
 124         }
 125
 126         if (!opt && map->n > 1 && sol) {
 127                 isl_int_init(o);
 128                 opt = &o;
 129         }
 130         if (map->n > 0)
 131                 isl_int_init(opt_i);
 132         if (map->n > 0 && opt_denom) {
 133                 isl_int_init(opt_denom_i);
 134                 isl_int_init(t);
 135         }
 136
 137         res = isl_basic_map_solve_lp(map->p[0], max, f, d,
 138                                         opt, opt_denom, sol);
 139         if (res == isl_lp_error || res == isl_lp_unbounded)
 140                 goto done;
 141
 142         if (sol)
 143                 *sol = NULL;
 144
 145         for (i = 1; i < map->n; ++i) {
 146                 isl_vec *sol_i = NULL;
 147                 enum isl_lp_result res_i;
 148                 int better;
 149
 150                 res_i = isl_basic_map_solve_lp(map->p[i], max, f, d,
 151                                             &opt_i,
 152                                             opt_denom ? &opt_denom_i : NULL,
 153                                             sol ? &sol_i : NULL);
 154                 if (res_i == isl_lp_error || res_i == isl_lp_unbounded) {
 155                         res = res_i;
 156                         goto done;
 157                 }
 158                 if (res_i == isl_lp_empty)
 159                         continue;
 160                 if (res == isl_lp_empty) {
 161                         better = 1;
 162                 } else if (!opt_denom) {
 163                         if (max)
 164                                 better = isl_int_gt(opt_i, *opt);
 165                         else
 166                                 better = isl_int_lt(opt_i, *opt);
 167                 } else {
 168                         isl_int_mul(t, opt_i, *opt_denom);
 169                         isl_int_submul(t, *opt, opt_denom_i);
 170                         if (max)
 171                                 better = isl_int_is_pos(t);
 172                         else
 173                                 better = isl_int_is_neg(t);
 174                 }
 175                 if (better) {
 176                         res = res_i;
 177                         if (opt)
 178                                 isl_int_set(*opt, opt_i);
 179                         if (opt_denom)
 180                                 isl_int_set(*opt_denom, opt_denom_i);
 181                         if (sol) {
 182                                 isl_vec_free(*sol);
 183                                 *sol = sol_i;
 184                         }
 185                 } else
 186                         isl_vec_free(sol_i);
 187         }
 188
 189 done:
 190         isl_vec_free(v);
 191         if (map->n > 0 && opt_denom) {
 192                 isl_int_clear(opt_denom_i);
 193                 isl_int_clear(t);
 194         }
 195         if (map->n > 0)
 196                 isl_int_clear(opt_i);
 197         if (opt == &o)
 198                 isl_int_clear(o);
 199         return res;
 200 }
 201
 202 enum isl_lp_result isl_set_solve_lp(__isl_keep isl_set *set, int max,
 203                                       isl_int *f, isl_int d, isl_int *opt,
 204                                       isl_int *opt_denom,
 205                                       struct isl_vec **sol)
 206 {
 207         return isl_map_solve_lp((struct isl_map *)set, max,
 208                                         f, d, opt, opt_denom, sol);
 209 }
 210
 211 /* Return the optimal (rational) value of "obj" over "bset", assuming
 212  * that "obj" and "bset" have aligned parameters and divs.
 213  * If "max" is set, then the maximal value is computed.
 214  * Otherwise, the minimal value is computed.
 215  *
 216  * Return infinity or negative infinity if the optimal value is unbounded and
 217  * NaN if "bset" is empty.
 218  *
 219  * Call isl_basic_set_solve_lp and translate the results.
 220  */
 221 static __isl_give isl_val *basic_set_opt_lp(
 222         __isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
 223 {
 224         isl_ctx *ctx;
 225         isl_val *res;
 226         enum isl_lp_result lp_res;
 227
 228         if (!bset || !obj)
 229                 return NULL;
 230
 231         ctx = isl_aff_get_ctx(obj);
 232         res = isl_val_alloc(ctx);
 233         if (!res)
 234                 return NULL;
 235         lp_res = isl_basic_set_solve_lp(bset, max, obj->v->el + 1,
 236                                         obj->v->el[0], &res->n, &res->d, NULL);
 237         if (lp_res == isl_lp_ok)
 238                 return isl_val_normalize(res);
 239         isl_val_free(res);
 240         if (lp_res == isl_lp_error)
 241                 return NULL;
 242         if (lp_res == isl_lp_empty)
 243                 return isl_val_nan(ctx);
 244         if (max)
 245                 return isl_val_infty(ctx);
 246         else
 247                 return isl_val_neginfty(ctx);
 248 }
 249
 250 /* Return the optimal (rational) value of "obj" over "bset", assuming
 251  * that "obj" and "bset" have aligned parameters.
 252  * If "max" is set, then the maximal value is computed.
 253  * Otherwise, the minimal value is computed.
 254  *
 255  * Return infinity or negative infinity if the optimal value is unbounded and
 256  * NaN if "bset" is empty.
 257  *
 258  * Align the divs of "bset" and "obj" and call basic_set_opt_lp.
 259  */
 260 static __isl_give isl_val *isl_basic_set_opt_lp_val_aligned(
 261         __isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
 262 {
 263         int *exp1 = NULL;
 264         int *exp2 = NULL;
 265         isl_ctx *ctx;
 266         isl_mat *bset_div = NULL;
 267         isl_mat *div = NULL;
 268         isl_val *res;
 269
 270         if (!bset || !obj)
 271                 return NULL;
 272
 273         ctx = isl_aff_get_ctx(obj);
 274         if (!isl_space_is_equal(bset->dim, obj->ls->dim))
 275                 isl_die(ctx, isl_error_invalid,
 276                         "spaces don't match", return NULL);
 277
 278         if (bset->n_div == 0 && obj->ls->div->n_row == 0)
 279                 return basic_set_opt_lp(bset, max, obj);
 280
 281         bset = isl_basic_set_copy(bset);
 282         obj = isl_aff_copy(obj);
 283
 284         bset_div = isl_basic_set_get_divs(bset);
 285         exp1 = isl_alloc_array(ctx, int, bset_div->n_row);
 286         exp2 = isl_alloc_array(ctx, int, obj->ls->div->n_row);
 287         if (!bset_div || !exp1 || !exp2)
 288                 goto error;
 289
 290         div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2);
 291
 292         bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1);
 293         obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2);
 294
 295         res = basic_set_opt_lp(bset, max, obj);
 296
 297         isl_mat_free(bset_div);
 298         isl_mat_free(div);
 299         free(exp1);
 300         free(exp2);
 301         isl_basic_set_free(bset);
 302         isl_aff_free(obj);
 303
 304         return res;
 305 error:
 306         isl_mat_free(div);
 307         isl_mat_free(bset_div);
 308         free(exp1);
 309         free(exp2);
 310         isl_basic_set_free(bset);
 311         isl_aff_free(obj);
 312         return NULL;
 313 }
 314
 315 /* Return the optimal (rational) value of "obj" over "bset".
 316  * If "max" is set, then the maximal value is computed.
 317  * Otherwise, the minimal value is computed.
 318  *
 319  * Return infinity or negative infinity if the optimal value is unbounded and
 320  * NaN if "bset" is empty.
 321  */
 322 static __isl_give isl_val *isl_basic_set_opt_lp_val(
 323         __isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
 324 {
 325         isl_val *res;
 326
 327         if (!bset || !obj)
 328                 return NULL;
 329
 330         if (isl_space_match(bset->dim, isl_dim_param,
 331                             obj->ls->dim, isl_dim_param))
 332                 return isl_basic_set_opt_lp_val_aligned(bset, max, obj);
 333
 334         bset = isl_basic_set_copy(bset);
 335         obj = isl_aff_copy(obj);
 336         bset = isl_basic_set_align_params(bset, isl_aff_get_domain_space(obj));
 337         obj = isl_aff_align_params(obj, isl_basic_set_get_space(bset));
 338
 339         res = isl_basic_set_opt_lp_val_aligned(bset, max, obj);
 340
 341         isl_basic_set_free(bset);
 342         isl_aff_free(obj);
 343
 344         return res;
 345 }
 346
 347 /* Return the minimal (rational) value of "obj" over "bset".
 348  *
 349  * Return negative infinity if the minimal value is unbounded and
 350  * NaN if "bset" is empty.
 351  */
 352 __isl_give isl_val *isl_basic_set_min_lp_val(__isl_keep isl_basic_set *bset,
 353         __isl_keep isl_aff *obj)
 354 {
 355         return isl_basic_set_opt_lp_val(bset, 0, obj);
 356 }
 357
 358 /* Return the maximal (rational) value of "obj" over "bset".
 359  *
 360  * Return infinity if the maximal value is unbounded and
 361  * NaN if "bset" is empty.
 362  */
 363 __isl_give isl_val *isl_basic_set_max_lp_val(__isl_keep isl_basic_set *bset,
 364         __isl_keep isl_aff *obj)
 365 {
 366         return isl_basic_set_opt_lp_val(bset, 1, obj);
 367 }