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