polysign_glpk.c

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