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