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