barvinok_bound: rewrite in C
[barvinok.git] / basis_reduction_pip.c
bloba0871a77ef74574aad0bcc7ca833f7b63aa0caac
1 #include <assert.h>
2 #include <piplib/piplibMP.h>
3 #include <barvinok/basis_reduction.h>
5 struct pip_lp {
6 Polyhedron *P;
7 Value *obj;
8 Matrix *eq;
9 int neq;
10 PipQuast *sol;
13 static struct pip_lp *init_lp(Polyhedron *P);
14 static void set_lp_obj(struct pip_lp *lp, Value *row, int dim);
15 static int solve_lp(struct pip_lp *lp);
16 static void get_obj_val(struct pip_lp* lp, mpq_t *F);
17 static void delete_lp(struct pip_lp *lp);
18 static int add_lp_row(struct pip_lp *lp, Value *row, int dim);
19 static void get_alpha(struct pip_lp* lp, int row, mpq_t *alpha);
20 static void del_lp_row(struct pip_lp *lp);
22 #define GBR_LP struct pip_lp
23 #define GBR_type mpq_t
24 #define GBR_init(v) mpq_init(v)
25 #define GBR_clear(v) mpq_clear(v)
26 #define GBR_set(a,b) mpq_set(a,b)
27 #define GBR_set_ui(a,b) mpq_set_ui(a,b,1)
28 #define GBR_mul(a,b,c) mpq_mul(a,b,c)
29 #define GBR_lt(a,b) (mpq_cmp(a,b) < 0)
30 #define GBR_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b))
31 #define GBR_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b))
32 #define GBR_lp_init(P) init_lp(P)
33 #define GBR_lp_set_obj(lp, obj, dim) set_lp_obj(lp, obj, dim)
34 #define GBR_lp_solve(lp) solve_lp(lp)
35 #define GBR_lp_get_obj_val(lp, F) get_obj_val(lp, F)
36 #define GBR_lp_delete(lp) delete_lp(lp)
37 #define GBR_lp_next_row(lp) lp->neq
38 #define GBR_lp_add_row(lp, row, dim) add_lp_row(lp, row, dim)
39 #define GBR_lp_get_alpha(lp, row, alpha) get_alpha(lp, row, alpha)
40 #define GBR_lp_del_row(lp) del_lp_row(lp);
41 #define Polyhedron_Reduced_Basis pip_Polyhedron_Reduced_Basis
42 #include "basis_reduction_templ.c"
44 #define ALLOC(type) (type*)malloc(sizeof(type))
46 static struct pip_lp *init_lp(Polyhedron *P)
48 struct pip_lp *lp = ALLOC(struct pip_lp);
49 lp->P = P;
50 lp->eq = Matrix_Alloc(P->Dimension, P->Dimension);
51 lp->neq = 0;
52 lp->sol = NULL;
54 return lp;
57 static void set_lp_obj(struct pip_lp *lp, Value *row, int dim)
59 assert(lp->P->Dimension == dim);
60 lp->obj = row;
63 static int solve_lp(struct pip_lp *lp)
65 PipOptions *options;
66 PipMatrix *domain;
67 unsigned dim = lp->P->Dimension;
68 int rows = 1 + 2*lp->P->NbConstraints + lp->neq;
69 int cols = 2*dim + 3;
70 int i, j;
72 if (lp->sol)
73 pip_quast_free(lp->sol);
75 domain = pip_matrix_alloc(rows, cols);
76 for (i = 0; i < lp->neq; ++i)
77 for (j = 0; j < dim; ++j) {
78 value_assign(domain->p[i][1+1+j], lp->eq->p[i][j]);
79 value_oppose(domain->p[i][1+1+dim+j], lp->eq->p[i][j]);
81 value_set_si(domain->p[lp->neq][1], -1);
82 for (j = 0; j < dim; ++j) {
83 value_assign(domain->p[lp->neq][1+1+j], lp->obj[j]);
84 value_oppose(domain->p[lp->neq][1+1+dim+j], lp->obj[j]);
86 for (i = 0; i < lp->P->NbConstraints; ++i) {
87 value_assign(domain->p[lp->neq+1+i][0], lp->P->Constraint[i][0]);
88 value_assign(domain->p[lp->neq+1+lp->P->NbConstraints+i][0],
89 lp->P->Constraint[i][0]);
90 for (j = 0; j < dim; ++j) {
91 value_assign(domain->p[lp->neq+1+i][2+j], lp->P->Constraint[i][1+j]);
92 value_assign(domain->p[lp->neq+1+lp->P->NbConstraints+i][2+dim+j],
93 lp->P->Constraint[i][1+j]);
95 value_assign(domain->p[lp->neq+1+i][cols-1], lp->P->Constraint[i][1+dim]);
96 value_assign(domain->p[lp->neq+1+lp->P->NbConstraints+i][cols-1],
97 lp->P->Constraint[i][1+dim]);
100 options = pip_options_init();
101 options->Urs_unknowns = -1;
102 options->Maximize = 1;
103 options->Nq = 0;
104 if (lp->neq)
105 options->Compute_dual = 1;
106 lp->sol = pip_solve(domain, NULL, -1, options);
107 pip_matrix_free(domain);
109 pip_options_free(options);
111 /* We only call this function on a polytope that is known
112 * to be (rationally) non-empty.
113 * However, if the test for rational emptiness was inaccurate
114 * (e.g., when using GLPK with big coefficients), the input
115 * polytope may be empty. It's better to just abort here
116 * so the user can select a different linear solver.
118 assert(lp->sol->list);
119 return value_zero_p(lp->sol->list->vector->the_deno[0]);
122 static void get_obj_val(struct pip_lp* lp, mpq_t *F)
124 value_assign(mpq_numref(*F), lp->sol->list->vector->the_vector[0]);
125 value_assign(mpq_denref(*F), lp->sol->list->vector->the_deno[0]);
128 static void delete_lp(struct pip_lp *lp)
130 if (lp->sol)
131 pip_quast_free(lp->sol);
132 Matrix_Free(lp->eq);
133 free(lp);
136 static int add_lp_row(struct pip_lp *lp, Value *row, int dim)
138 assert(lp->P->Dimension == dim);
139 Vector_Copy(row, lp->eq->p[lp->neq], dim);
140 return lp->neq++;
143 static void get_alpha(struct pip_lp* lp, int row, mpq_t *alpha)
145 PipList *list;
146 int i = 0;
148 assert(lp->sol->next_then);
149 assert(lp->sol->next_then->list);
151 for (i = 0, list = lp->sol->next_then->list; i < row; ++i, list = list->next)
154 value_assign(mpq_numref(*alpha), list->vector->the_vector[0]);
155 value_assign(mpq_denref(*alpha), list->vector->the_deno[0]);
158 static void del_lp_row(struct pip_lp *lp)
160 assert(lp->neq > 0);
161 lp->neq--;