basis_reduction_glpk.c

   1 #include <assert.h>
   2 #include <math.h>
   3 #include <glpk.h>
   4 #include <barvinok/basis_reduction.h>
   5
   6 static glp_prob *init_lp(Polyhedron *P);
   7 static void set_lp_obj(glp_prob *lp, Value *row, int dim);
   8 static int solve_lp(glp_prob *lp);
   9 static void get_obj_val(glp_prob* lp, double *F);
  10 static int add_lp_row(glp_prob *lp, Value *row, int dim);
  11 static void get_alpha(glp_prob* lp, int row, double *alpha);
  12 static void del_lp_row(glp_prob *lp);
  13
  14 #define GBR_LP                              glp_prob
  15 #define GBR_type                            double
  16 #define GBR_init(v)                         do { } while(0)
  17 #define GBR_clear(v)                        do { } while(0)
  18 #define GBR_set(a,b)                        a = b
  19 #define GBR_set_ui(a,b)                     a = b
  20 #define GBR_mul(a,b,c)                      a = b * c
  21 #define GBR_lt(a,b)                         (a < b)
  22 #define GBR_floor(a,b)                      value_set_si(a,(int)floor(b+1e-10))
  23 #define GBR_ceil(a,b)                       value_set_si(a,(int)ceil(b-1e-10))
  24 #define GBR_lp_init(P)                      init_lp(P)
  25 #define GBR_lp_set_obj(lp, obj, dim)        set_lp_obj(lp, obj, dim)
  26 #define GBR_lp_solve(lp)                    solve_lp(lp)
  27 #define GBR_lp_get_obj_val(lp, F)           get_obj_val(lp, F)
  28 #define GBR_lp_delete(lp)                   glp_delete_prob(lp)
  29 #define GBR_lp_next_row(lp)                 glp_get_num_rows(lp)+1
  30 #define GBR_lp_add_row(lp, row, dim)        add_lp_row(lp, row, dim)
  31 #define GBR_lp_get_alpha(lp, row, alpha)    get_alpha(lp, row, alpha)
  32 #define GBR_lp_del_row(lp)                  del_lp_row(lp);
  33 #define Polyhedron_Reduced_Basis            glpk_Polyhedron_Reduced_Basis
  34 #include "basis_reduction_templ.c"
  35
  36 static glp_prob *init_lp(Polyhedron *P)
  37 {
  38     glp_prob *lp;
  39     int *ind;
  40     double *val;
  41     int i, j, k, l;
  42
  43     ind = ALLOCN(int, 1+P->Dimension);
  44     val = ALLOCN(double, 1+P->Dimension);
  45     lp = glp_create_prob();
  46     glp_set_obj_dir(lp, GLP_MAX);
  47     glp_add_rows(lp, 2*P->NbConstraints);
  48     glp_add_cols(lp, 2*P->Dimension);
  49     for (i = 0; i < 2; ++i) {
  50         for (j = 0; j < P->NbConstraints; ++j) {
  51             int type = j < P->NbEq ? GLP_FX : GLP_LO;
  52             for (k = 0, l = 0; k < P->Dimension; ++k) {
  53                 if (value_zero_p(P->Constraint[j][1+k]))
  54                     continue;
  55                 ind[1+l] = 1+i*P->Dimension+k;
  56                 val[1+l] = VALUE_TO_DOUBLE(P->Constraint[j][1+k]);
  57                 ++l;
  58             }
  59             glp_set_mat_row(lp, 1+i*P->NbConstraints+j, l, ind, val);
  60             glp_set_row_bnds(lp, 1+i*P->NbConstraints+j, type,
  61                              -VALUE_TO_DOUBLE(P->Constraint[j][1+P->Dimension]), 0);
  62         }
  63         for (k = 0, l = 0; k < P->Dimension; ++k) {
  64             glp_set_col_bnds(lp, 1+i*P->Dimension+k, GLP_FR, 0, 0);
  65         }
  66     }
  67     free(ind);
  68     free(val);
  69     return lp;
  70 }
  71
  72 static void set_lp_obj(glp_prob *lp, Value *row, int dim)
  73 {
  74     int j;
  75     for (j = 0; j < dim; ++j) {
  76         glp_set_obj_coef(lp, 1+j, VALUE_TO_DOUBLE(row[j]));
  77         glp_set_obj_coef(lp, 1+dim+j, -VALUE_TO_DOUBLE(row[j]));
  78     }
  79 }
  80
  81 static int solve_lp(glp_prob *lp)
  82 {
  83     glp_smcp parm;
  84     glp_adv_basis(lp, 0);
  85     glp_init_smcp(&parm);
  86     parm.msg_lev = GLP_MSG_OFF;
  87     glp_simplex(lp, &parm);
  88     return glp_get_status(lp) == GLP_UNBND;
  89 }
  90
  91 static void get_obj_val(glp_prob* lp, double *F)
  92 {
  93     assert(glp_get_status(lp) == GLP_OPT);
  94     *F = glp_get_obj_val(lp);
  95     assert(*F > -1e-10);
  96     if (*F < 0)
  97         *F = 0;
  98 }
  99
 100 static int add_lp_row(glp_prob *lp, Value *row, int dim)
 101 {
 102     int j, l;
 103     int nr = glp_add_rows(lp, 1);
 104     int *ind;
 105     double *val;
 106
 107     ind = ALLOCN(int, 1+2*dim);
 108     val = ALLOCN(double, 1+2*dim);
 109     for (j = 0, l = 0; j < dim; ++j) {
 110         if (value_zero_p(row[j]))
 111             continue;
 112         ind[1+l] = 1+j;
 113         val[1+l] = VALUE_TO_DOUBLE(row[j]);
 114         ind[1+l+1] = 1+dim+j;
 115         val[1+l+1] = -VALUE_TO_DOUBLE(row[j]);
 116         l += 2;
 117     }
 118     glp_set_mat_row(lp, nr, l, ind, val);
 119     glp_set_row_bnds(lp, nr, GLP_FX, 0, 0);
 120     free(ind);
 121     free(val);
 122
 123     return nr;
 124 }
 125
 126 static void get_alpha(glp_prob* lp, int row, double *alpha)
 127 {
 128     *alpha = -glp_get_row_dual(lp, row);
 129 }
 130
 131 static void del_lp_row(glp_prob *lp)
 132 {
 133     int rows[2];
 134     rows[1] = glp_get_num_rows(lp);
 135     glp_del_rows(lp, 1, rows);
 136 }
 137