isl_lp.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/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
  18 enum isl_lp_result isl_tab_solve_lp(struct isl_basic_map *bmap, int maximize,
  19                                       isl_int *f, isl_int denom, isl_int *opt,
  20                                       isl_int *opt_denom,
  21                                       struct isl_vec **sol)
  22 {
  23         struct isl_tab *tab;
  24         enum isl_lp_result res;
  25         unsigned dim = isl_basic_map_total_dim(bmap);
  26
  27         if (maximize)
  28                 isl_seq_neg(f, f, 1 + dim);
  29
  30         bmap = isl_basic_map_gauss(bmap, NULL);
  31         tab = isl_tab_from_basic_map(bmap, 0);
  32         res = isl_tab_min(tab, f, denom, opt, opt_denom, 0);
  33         if (res == isl_lp_ok && sol) {
  34                 *sol = isl_tab_get_sample_value(tab);
  35                 if (!*sol)
  36                         res = isl_lp_error;
  37         }
  38         isl_tab_free(tab);
  39
  40         if (maximize)
  41                 isl_seq_neg(f, f, 1 + dim);
  42         if (maximize && opt)
  43                 isl_int_neg(*opt, *opt);
  44
  45         return res;
  46 }
  47
  48 /* Given a basic map "bmap" and an affine combination of the variables "f"
  49  * with denominator "denom", set *opt / *opt_denom to the minimal
  50  * (or maximal if "maximize" is true) value attained by f/d over "bmap",
  51  * assuming the basic map is not empty and the expression cannot attain
  52  * arbitrarily small (or large) values.
  53  * If opt_denom is NULL, then *opt is rounded up (or down)
  54  * to the nearest integer.
  55  * The return value reflects the nature of the result (empty, unbounded,
  56  * minmimal or maximal value returned in *opt).
  57  */
  58 enum isl_lp_result isl_basic_map_solve_lp(struct isl_basic_map *bmap, int max,
  59                                       isl_int *f, isl_int d, isl_int *opt,
  60                                       isl_int *opt_denom,
  61                                       struct isl_vec **sol)
  62 {
  63         if (sol)
  64                 *sol = NULL;
  65
  66         if (!bmap)
  67                 return isl_lp_error;
  68
  69         switch (bmap->ctx->opt->lp_solver) {
  70         case ISL_LP_PIP:
  71                 return isl_pip_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
  72         case ISL_LP_TAB:
  73                 return isl_tab_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
  74         default:
  75                 return isl_lp_error;
  76         }
  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((struct isl_basic_map *)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((struct isl_map *)set, max,
 203                                         f, d, opt, opt_denom, sol);
 204 }