README: update note on shared libraries and NTL
[barvinok.git] / polysign_glpk.c
blob734cca2d7bb1260c93487d8a420dc5fea43961e8
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"
8 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
10 #define EMPTY_DOMAIN -2
12 static glp_prob *solve_lp(int dir, Matrix *C, Value *f, Value denom)
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;
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);
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;
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);
41 for (k = 0, l = 0; k < dim; ++k) {
42 glp_set_col_bnds(lp, 1+k, GLP_FR, 0, 0);
44 free(ind);
45 free(val);
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));
54 glp_adv_basis(lp, 0);
55 glp_init_smcp(&parm);
56 parm.msg_lev = GLP_MSG_OFF;
57 glp_simplex(lp, &parm);
59 return lp;
62 static enum lp_result constraints_affine_minmax(int dir, Matrix *C,
63 Value *f, Value denom, Value *opt)
65 enum lp_result res = lp_ok;
66 glp_prob *lp = solve_lp(dir, C, f, denom);
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);
84 glp_delete_prob(lp);
85 return res;
88 static int constraints_affine_minmax_sign(int dir, Matrix *C, Matrix *T,
89 int rational)
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);
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;
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);
125 glp_delete_prob(lp);
126 return sign;
129 enum order_sign glpk_polyhedron_affine_sign_0D(Polyhedron *D, Matrix *T)
132 int sgn;
134 POL_ENSURE_VERTICES(D);
136 if (emptyQ(D))
137 return order_undefined;
139 sgn = value_sign(T->p[0][0]);
141 return sgn < 0 ? order_lt : sgn == 0 ? order_eq : order_gt;
144 enum order_sign glpk_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
145 struct barvinok_options *options)
147 int rational = !POL_ISSET(options->MaxRays, POL_INTEGER);
148 Matrix M;
150 if (emptyQ2(D))
151 return order_undefined;
153 if (D->Dimension == 0)
154 return glpk_polyhedron_affine_sign_0D(D, T);
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;
175 enum lp_result glpk_constraints_opt(Matrix *C, Value *obj, Value denom,
176 enum lp_dir dir, Value *opt)
178 int glpk_dir = dir == lp_min ? GLP_MIN : GLP_MAX;
179 return constraints_affine_minmax(glpk_dir, C, obj, denom, opt);