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