polysign_glpk.c

   1 #include <assert.h>
   2 #include <math.h>
   3 #include <glpk.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/util.h>
   6 #include "polysign.h"
   7
   8 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
   9
  10 #define EMPTY_DOMAIN    -2
  11
  12 static LPX *solve_lp(int dir, Polyhedron *P, Value *f, Value denom)
  13 {
  14     LPX *lp;
  15     int *ind;
  16     double *val;
  17     int j, k, l;
  18
  19     ind = ALLOCN(int, 1+P->Dimension);
  20     val = ALLOCN(double, 1+P->Dimension);
  21     lp = lpx_create_prob();
  22     lpx_set_obj_dir(lp, dir);
  23     lpx_add_rows(lp, P->NbConstraints);
  24     lpx_add_cols(lp, P->Dimension);
  25
  26     for (j = 0; j < P->NbConstraints; ++j) {
  27         int type = j < P->NbEq ? LPX_FX : LPX_LO;
  28         for (k = 0, l = 0; k < P->Dimension; ++k) {
  29             if (value_zero_p(P->Constraint[j][1+k]))
  30                 continue;
  31             ind[1+l] = 1+k;
  32             val[1+l] = VALUE_TO_DOUBLE(P->Constraint[j][1+k]);
  33             ++l;
  34         }
  35         lpx_set_mat_row(lp, 1+j, l, ind, val);
  36         lpx_set_row_bnds(lp, 1+j, type,
  37                          -VALUE_TO_DOUBLE(P->Constraint[j][1+P->Dimension]), 0);
  38     }
  39     for (k = 0, l = 0; k < P->Dimension; ++k) {
  40         lpx_set_col_bnds(lp, 1+k, LPX_FR, 0, 0);
  41     }
  42     free(ind);
  43     free(val);
  44     lpx_set_int_parm(lp, LPX_K_MSGLEV, 0);
  45
  46     /* objective function */
  47     for (j = 0; j < P->Dimension; ++j)
  48         lpx_set_obj_coef(lp, 1+j, VALUE_TO_DOUBLE(f[j]) /
  49                                     VALUE_TO_DOUBLE(denom));
  50     lpx_set_obj_coef(lp, 0, VALUE_TO_DOUBLE(f[P->Dimension]) /
  51                                 VALUE_TO_DOUBLE(denom));
  52
  53     lpx_adv_basis(lp);
  54     lpx_simplex(lp);
  55
  56     return lp;
  57 }
  58
  59 static enum lp_result polyhedron_affine_minmax(int dir, Polyhedron *P,
  60                                               Value *f, Value denom, Value *opt)
  61 {
  62     enum lp_result res = lp_ok;
  63     LPX *lp = solve_lp(dir, P, f, denom);
  64
  65     switch (lpx_get_status(lp)) {
  66     case LPX_OPT:
  67         if (dir == LPX_MIN)
  68             value_set_si(*opt, (int)ceil(lpx_get_obj_val(lp)-1e-10));
  69         else
  70             value_set_si(*opt, (int)floor(lpx_get_obj_val(lp)+1e-10));
  71         break;
  72     case LPX_UNBND:
  73         res = lp_unbounded;
  74         break;
  75     case LPX_NOFEAS:
  76         res = lp_empty;
  77         break;
  78     default:
  79         assert(0);
  80     }
  81     lpx_delete_prob(lp);
  82     return res;
  83 }
  84
  85 static int polyhedron_affine_minmax_sign(int dir, Polyhedron *P, Matrix *T,
  86                                          int rational)
  87 {
  88     LPX *lp;
  89     int sign;
  90     double opt;
  91     assert(P->Dimension == T->NbColumns-1);
  92     assert(T->NbRows == 2);
  93
  94     lp = solve_lp(dir, P, T->p[0], T->p[1][P->Dimension]);
  95     switch (lpx_get_status(lp)) {
  96     case LPX_OPT:
  97         opt = lpx_get_obj_val(lp);
  98         if (rational) {
  99             sign = opt < 0 ? -1 : opt > 0 ? 1 : 0;
 100         } else {
 101             if (opt < -0.5/VALUE_TO_DOUBLE(T->p[1][P->Dimension]))
 102                 sign = -1;
 103             else if (opt > 0.5/VALUE_TO_DOUBLE(T->p[1][P->Dimension]))
 104                 sign = 1;
 105             else
 106                 sign = 0;
 107         }
 108         break;
 109     case LPX_UNBND:
 110         if (dir == LPX_MIN)
 111             sign = -1;
 112         else
 113             sign = 1;
 114         break;
 115     case LPX_NOFEAS:
 116         sign = EMPTY_DOMAIN;
 117         break;
 118     default:
 119         assert(0);
 120     }
 121     lpx_delete_prob(lp);
 122     return sign;
 123 }
 124
 125 enum order_sign glpk_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
 126                                             struct barvinok_options *options)
 127 {
 128     int rational = !POL_ISSET(options->MaxRays, POL_INTEGER);
 129
 130     if (emptyQ2(D))
 131         return order_undefined;
 132
 133     int min = polyhedron_affine_minmax_sign(LPX_MIN, D, T, rational);
 134     if (min == EMPTY_DOMAIN)
 135         return order_undefined;
 136     if (min > 0)
 137         return order_gt;
 138     int max = polyhedron_affine_minmax_sign(LPX_MAX, D, T, rational);
 139     assert(max != EMPTY_DOMAIN);
 140     if (max < 0)
 141         return order_lt;
 142     if (min == max)
 143         return order_eq;
 144     if (max == 0)
 145         return order_le;
 146     if (min == 0)
 147         return order_ge;
 148     return order_unknown;
 149 }
 150
 151 enum lp_result glpk_polyhedron_range(Polyhedron *D, Value *obj, Value denom,
 152                                 Value *min, Value *max,
 153                                 struct barvinok_options *options)
 154 {
 155     enum lp_result res;
 156
 157     if (emptyQ2(D))
 158         return lp_empty;
 159
 160     res = polyhedron_affine_minmax(LPX_MIN, D, obj, denom, min);
 161     if (res != lp_ok)
 162         return res;
 163     res = polyhedron_affine_minmax(LPX_MAX, D, obj, denom, max);
 164     return res;
 165 }