barvinok 0.41.8
[barvinok.git] / basis_reduction_glpk.c
blob5806dbe5a3bc46641ac3e53b2364f99f11911052
1 #include <assert.h>
2 #include <math.h>
3 #include <glpk.h>
4 #include <barvinok/basis_reduction.h>
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);
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"
36 static glp_prob *init_lp(Polyhedron *P)
38 glp_prob *lp;
39 int *ind;
40 double *val;
41 int i, j, k, l;
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;
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);
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);
67 free(ind);
68 free(val);
69 return lp;
72 static void set_lp_obj(glp_prob *lp, Value *row, int dim)
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]));
81 static int solve_lp(glp_prob *lp)
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;
91 static void get_obj_val(glp_prob* lp, double *F)
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;
100 static int add_lp_row(glp_prob *lp, Value *row, int dim)
102 int j, l;
103 int nr = glp_add_rows(lp, 1);
104 int *ind;
105 double *val;
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;
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);
123 return nr;
126 static void get_alpha(glp_prob* lp, int row, double *alpha)
128 *alpha = -glp_get_row_dual(lp, row);
131 static void del_lp_row(glp_prob *lp)
133 int rows[2];
134 rows[1] = glp_get_num_rows(lp);
135 glp_del_rows(lp, 1, rows);