basis_reduction_glpk.c

   1 #include <math.h>
   2 #include <glpk.h>
   3 #include <barvinok/basis_reduction.h>
   4
   5 static LPX *init_lp(Polyhedron *P);
   6 static void set_lp_obj(LPX *lp, Value *row, int dim);
   7 static void solve_lp(LPX *lp);
   8 static void get_obj_val(LPX* lp, double *F);
   9 static int add_lp_row(LPX *lp, Value *row, int dim);
  10 static void get_alpha(LPX* lp, int row, double *alpha);
  11 static void del_lp_row(LPX *lp);
  12
  13 #define GBR_LP                              LPX
  14 #define GBR_type                            double
  15 #define GBR_init(v)                         do { } while(0)
  16 #define GBR_clear(v)                        do { } while(0)
  17 #define GBR_set(a,b)                        a = b
  18 #define GBR_set_ui(a,b)                     a = b
  19 #define GBR_mul(a,b,c)                      a = b * c
  20 #define GBR_lt(a,b)                         (a < b)
  21 #define GBR_floor(a,b)                      value_set_si(a,(int)floor(b+1e-10))
  22 #define GBR_ceil(a,b)                       value_set_si(a,(int)ceil(b-1e-10))
  23 #define GBR_lp_init(P)                      init_lp(P)
  24 #define GBR_lp_set_obj(lp, obj, dim)        set_lp_obj(lp, obj, dim)
  25 #define GBR_lp_solve(lp)                    solve_lp(lp)
  26 #define GBR_lp_get_obj_val(lp, F)           get_obj_val(lp, F)
  27 #define GBR_lp_delete(lp)                   lpx_delete_prob(lp)
  28 #define GBR_lp_next_row(lp)                 lpx_get_num_rows(lp)+1
  29 #define GBR_lp_add_row(lp, row, dim)        add_lp_row(lp, row, dim)
  30 #define GBR_lp_get_alpha(lp, row, alpha)    get_alpha(lp, row, alpha)
  31 #define GBR_lp_del_row(lp)                  del_lp_row(lp);
  32 #define Polyhedron_Reduced_Basis            glpk_Polyhedron_Reduced_Basis
  33 #include "basis_reduction_templ.c"
  34
  35 static LPX *init_lp(Polyhedron *P)
  36 {
  37     LPX *lp;
  38     int *ind;
  39     double *val;
  40     int i, j, k, l;
  41
  42     ind = ALLOCN(int, 1+P->Dimension);
  43     val = ALLOCN(double, 1+P->Dimension);
  44     lp = lpx_create_prob();
  45     lpx_set_obj_dir(lp, LPX_MAX);
  46     lpx_add_rows(lp, 2*P->NbConstraints);
  47     lpx_add_cols(lp, 2*P->Dimension);
  48     for (i = 0; i < 2; ++i) {
  49         for (j = 0; j < P->NbConstraints; ++j) {
  50             for (k = 0, l = 0; k < P->Dimension; ++k) {
  51                 if (value_zero_p(P->Constraint[j][1+k]))
  52                     continue;
  53                 ind[1+l] = 1+i*P->Dimension+k;
  54                 val[1+l] = VALUE_TO_DOUBLE(P->Constraint[j][1+k]);
  55                 ++l;
  56             }
  57             lpx_set_mat_row(lp, 1+i*P->NbConstraints+j, l, ind, val);
  58             lpx_set_row_bnds(lp, 1+i*P->NbConstraints+j, LPX_LO,
  59                              -VALUE_TO_DOUBLE(P->Constraint[j][1+P->Dimension]), 0);
  60         }
  61         for (k = 0, l = 0; k < P->Dimension; ++k) {
  62             lpx_set_col_bnds(lp, 1+i*P->Dimension+k, LPX_FR, 0, 0);
  63         }
  64     }
  65     free(ind);
  66     free(val);
  67     lpx_set_int_parm(lp, LPX_K_MSGLEV, 0);
  68     return lp;
  69 }
  70
  71 static void set_lp_obj(LPX *lp, Value *row, int dim)
  72 {
  73     int j;
  74     for (j = 0; j < dim; ++j) {
  75         lpx_set_obj_coef(lp, 1+j, VALUE_TO_DOUBLE(row[j]));
  76         lpx_set_obj_coef(lp, 1+dim+j, -VALUE_TO_DOUBLE(row[j]));
  77     }
  78 }
  79
  80 static void solve_lp(LPX *lp)
  81 {
  82     lpx_adv_basis(lp);
  83     lpx_simplex(lp);
  84 }
  85
  86 static void get_obj_val(LPX* lp, double *F)
  87 {
  88     assert(lpx_get_status(lp) == LPX_OPT);
  89     *F = lpx_get_obj_val(lp);
  90     assert(*F > -1e-10);
  91     if (*F < 0)
  92         *F = 0;
  93 }
  94
  95 static int add_lp_row(LPX *lp, Value *row, int dim)
  96 {
  97     int j, l;
  98     int nr = lpx_add_rows(lp, 1);
  99     int *ind;
 100     double *val;
 101
 102     ind = ALLOCN(int, 1+2*dim);
 103     val = ALLOCN(double, 1+2*dim);
 104     for (j = 0, l = 0; j < dim; ++j) {
 105         if (value_zero_p(row[j]))
 106             continue;
 107         ind[1+l] = 1+j;
 108         val[1+l] = VALUE_TO_DOUBLE(row[j]);
 109         ind[1+l+1] = 1+dim+j;
 110         val[1+l+1] = -VALUE_TO_DOUBLE(row[j]);
 111         l += 2;
 112     }
 113     lpx_set_mat_row(lp, nr, l, ind, val);
 114     lpx_set_row_bnds(lp, nr, LPX_FX, 0, 0);
 115     free(ind);
 116     free(val);
 117
 118     return nr;
 119 }
 120
 121 static void get_alpha(LPX* lp, int row, double *alpha)
 122 {
 123     *alpha = -lpx_get_row_dual(lp, row);
 124 }
 125
 126 static void del_lp_row(LPX *lp)
 127 {
 128     int rows[2];
 129     rows[1] = lpx_get_num_rows(lp);
 130     lpx_del_rows(lp, 1, rows);
 131 }
 132