update pet for explicitly linking in gmp
[barvinok.git] / basis_reduction_cdd.cc
blobc3dc8b75d39cb82a45ff58719d250ea49623e60b
1 #include <assert.h>
2 #define GMPRATIONAL
3 #include <setoper.h>
4 #include <cdd.h>
5 #include <barvinok/basis_reduction.h>
6 #include "initcdd.h"
8 struct CDD_LP {
9 Polyhedron *P;
10 Value *obj;
11 dd_LPType *lp;
12 Matrix *eq;
13 int neq;
15 CDD_LP(Polyhedron *P);
16 void set_obj(Value *row, int dim);
17 int solve();
18 void get_obj_val(mytype *F);
19 int add_row(Value *row, int dim);
20 void get_alpha(int row, mytype *alpha);
21 void del_row();
22 ~CDD_LP() {
23 if (lp)
24 dd_FreeLPData(lp);
25 Matrix_Free(eq);
29 #define GBR_LP CDD_LP
30 #define GBR_type mytype
31 #define GBR_init(v) dd_init(v)
32 #define GBR_clear(v) dd_clear(v)
33 #define GBR_set(a,b) dd_set(a,b)
34 #define GBR_set_ui(a,b) mpq_set_ui(a,b,1)
35 #define GBR_mul(a,b,c) dd_mul(a,b,c)
36 #define GBR_lt(a,b) (dd_cmp(a,b) < 0)
37 #define GBR_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b))
38 #define GBR_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b))
39 #define GBR_lp_init(P) new CDD_LP(P)
40 #define GBR_lp_set_obj(lp, obj, dim) lp->set_obj(obj, dim)
41 #define GBR_lp_solve(lp) lp->solve()
42 #define GBR_lp_get_obj_val(lp, F) lp->get_obj_val(F)
43 #define GBR_lp_delete(lp) delete lp
44 #define GBR_lp_next_row(lp) lp->neq
45 #define GBR_lp_add_row(lp, row, dim) lp->add_row(row, dim)
46 #define GBR_lp_get_alpha(lp, row, alpha) lp->get_alpha(row, alpha)
47 #define GBR_lp_del_row(lp) lp->del_row()
48 #define Polyhedron_Reduced_Basis cdd_Polyhedron_Reduced_Basis
49 #include "basis_reduction_templ.c"
51 CDD_LP::CDD_LP(Polyhedron *P)
53 INIT_CDD;
55 this->P = P;
56 lp = NULL;
57 eq = Matrix_Alloc(P->Dimension, P->Dimension);
58 neq = 0;
61 void CDD_LP::set_obj(Value *row, int dim)
63 assert(P->Dimension == dim);
64 obj = row;
67 int CDD_LP::solve()
69 if (lp)
70 dd_FreeLPData(lp);
72 dd_rowrange irev = 2 * P->NbConstraints + neq;
73 dd_rowrange rows = irev + 2 * P->NbEq + neq + 1;
74 dd_colrange cols = 1 + 2*P->Dimension;
75 lp = dd_CreateLPData(dd_LPmax, dd_Rational, rows, cols);
76 lp->Homogeneous = dd_FALSE;
77 lp->objective = dd_LPmax;
78 for (int i = 0; i < 2; ++i) {
79 for (dd_rowrange j = 0; j < P->NbConstraints; ++j) {
80 dd_rowrange row = j+i*P->NbConstraints;
81 for (dd_colrange k = 0; k < P->Dimension; ++k) {
82 dd_colrange col = 1+k+i*P->Dimension;
83 mpq_set_z(lp->A[row][col], P->Constraint[j][1+k]);
85 mpq_set_z(lp->A[row][0], P->Constraint[j][1+P->Dimension]);
86 if (j < P->NbEq) {
87 set_addelem(lp->equalityset, row+1);
88 for (dd_colrange k = 0; k < P->Dimension; ++k) {
89 dd_colrange col = 1+k+i*P->Dimension;
90 dd_neg(lp->A[irev][col], lp->A[row][col]);
92 dd_neg(lp->A[irev][0], lp->A[row][0]);
93 ++irev;
97 for (int i = 0; i < neq; ++i) {
98 dd_rowrange row = 2*P->NbConstraints+i;
99 set_addelem(lp->equalityset, row+1);
100 for (int j = 0; j < P->Dimension; ++j) {
101 mpq_set_z(lp->A[row][1+j], eq->p[i][j]);
102 dd_neg(lp->A[row][1+P->Dimension+j], lp->A[row][1+j]);
103 dd_neg(lp->A[irev][1+j], lp->A[row][1+j]);
104 dd_set(lp->A[irev][1+P->Dimension+j], lp->A[row][1+j]);
106 ++irev;
108 /* objective function */
109 for (dd_colrange k = 0; k < P->Dimension; ++k) {
110 mpq_set_z(lp->A[rows-1][1+k], obj[k]);
111 dd_neg(lp->A[rows-1][1+P->Dimension+k], lp->A[rows-1][1+k]);
113 dd_ErrorType err = dd_NoError;
114 dd_LPSolve(lp, dd_DualSimplex, &err);
115 assert(err == dd_NoError);
116 /* We only call this function on a polytope that is known
117 * to be (rationally) non-empty.
118 * A bug in cddlib 0.94d can sometimes make cddlib claim
119 * that the polytope is empty. We can't continue if we
120 * hit this bug.
122 assert(lp->LPS != dd_Inconsistent);
124 return lp->LPS == dd_DualInconsistent;
127 void CDD_LP::get_obj_val(mytype *F)
129 dd_set(*F, lp->optvalue);
130 assert(dd_sgn(*F) >= 0);
133 int CDD_LP::add_row(Value *row, int dim)
135 assert(P->Dimension == dim);
136 Vector_Copy(row, eq->p[neq], dim);
137 return neq++;
140 void CDD_LP::get_alpha(int row, mytype *alpha)
142 row += 2 * P->NbConstraints;
143 /* Value for basic variable is zero, right ?
145 dd_set(*alpha, dd_zero);
146 for (int j = 1; j < lp->d; j++) {
147 if (lp->nbindex[j+1] == row+1) {
148 dd_set(*alpha, lp->dsol[j]);
149 break;
150 } else if (lp->nbindex[j+1] == row+neq+2*P->NbEq+1) {
151 dd_neg(*alpha, lp->dsol[j]);
152 break;
157 void CDD_LP::del_row()
159 assert(neq > 0);
160 neq--;