volume.c: allow computation of lower and upper bound volumes
[barvinok.git] / basis_reduction_cdd.cc
blob64ca5bdb710b2063c7db06ec13d8d285b905773b
1 #define GMPRATIONAL
2 #include <setoper.h>
3 #include <cdd.h>
4 #include <barvinok/basis_reduction.h>
5 #include "initcdd.h"
7 struct CDD_LP {
8 Polyhedron *P;
9 Value *obj;
10 dd_LPType *lp;
11 Matrix *eq;
12 int neq;
14 CDD_LP(Polyhedron *P);
15 void set_obj(Value *row, int dim);
16 void solve();
17 void get_obj_val(mytype *F);
18 int add_row(Value *row, int dim);
19 void get_alpha(int row, mytype *alpha);
20 void del_row();
21 ~CDD_LP() {
22 if (lp)
23 dd_FreeLPData(lp);
24 Matrix_Free(eq);
28 #define GBR_LP CDD_LP
29 #define GBR_type mytype
30 #define GBR_init(v) dd_init(v)
31 #define GBR_clear(v) dd_clear(v)
32 #define GBR_set(a,b) dd_set(a,b)
33 #define GBR_set_ui(a,b) mpq_set_ui(a,b,1)
34 #define GBR_mul(a,b,c) dd_mul(a,b,c)
35 #define GBR_lt(a,b) (dd_cmp(a,b) < 0)
36 #define GBR_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b))
37 #define GBR_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b))
38 #define GBR_lp_init(P) new CDD_LP(P)
39 #define GBR_lp_set_obj(lp, obj, dim) lp->set_obj(obj, dim)
40 #define GBR_lp_solve(lp) lp->solve()
41 #define GBR_lp_get_obj_val(lp, F) lp->get_obj_val(F)
42 #define GBR_lp_delete(lp) delete lp
43 #define GBR_lp_next_row(lp) lp->neq
44 #define GBR_lp_add_row(lp, row, dim) lp->add_row(row, dim)
45 #define GBR_lp_get_alpha(lp, row, alpha) lp->get_alpha(row, alpha)
46 #define GBR_lp_del_row(lp) lp->del_row()
47 #define Polyhedron_Reduced_Basis cdd_Polyhedron_Reduced_Basis
48 #include "basis_reduction_templ.c"
50 CDD_LP::CDD_LP(Polyhedron *P)
52 INIT_CDD;
54 this->P = P;
55 lp = NULL;
56 eq = Matrix_Alloc(P->Dimension, P->Dimension);
57 neq = 0;
60 void CDD_LP::set_obj(Value *row, int dim)
62 assert(P->Dimension == dim);
63 obj = row;
66 void CDD_LP::solve()
68 if (lp)
69 dd_FreeLPData(lp);
71 dd_rowrange irev = 2 * P->NbConstraints + neq;
72 dd_rowrange rows = irev + 2 * P->NbEq + neq + 1;
73 dd_colrange cols = 1 + 2*P->Dimension;
74 lp = dd_CreateLPData(dd_LPmax, dd_Rational, rows, cols);
75 lp->Homogeneous = dd_FALSE;
76 lp->objective = dd_LPmax;
77 for (int i = 0; i < 2; ++i) {
78 for (dd_rowrange j = 0; j < P->NbConstraints; ++j) {
79 dd_rowrange row = j+i*P->NbConstraints;
80 for (dd_colrange k = 0; k < P->Dimension; ++k) {
81 dd_colrange col = 1+k+i*P->Dimension;
82 mpq_set_z(lp->A[row][col], P->Constraint[j][1+k]);
84 mpq_set_z(lp->A[row][0], P->Constraint[j][1+P->Dimension]);
85 if (j < P->NbEq) {
86 set_addelem(lp->equalityset, row+1);
87 for (dd_colrange k = 0; k < P->Dimension; ++k) {
88 dd_colrange col = 1+k+i*P->Dimension;
89 dd_neg(lp->A[irev][col], lp->A[row][col]);
91 dd_neg(lp->A[irev][0], lp->A[row][0]);
92 ++irev;
96 for (int i = 0; i < neq; ++i) {
97 dd_rowrange row = 2*P->NbConstraints+i;
98 set_addelem(lp->equalityset, row+1);
99 for (int j = 0; j < P->Dimension; ++j) {
100 mpq_set_z(lp->A[row][1+j], eq->p[i][j]);
101 dd_neg(lp->A[row][1+P->Dimension+j], lp->A[row][1+j]);
102 dd_neg(lp->A[irev][1+j], lp->A[row][1+j]);
103 dd_set(lp->A[irev][1+P->Dimension+j], lp->A[row][1+j]);
105 ++irev;
107 /* objective function */
108 for (dd_colrange k = 0; k < P->Dimension; ++k) {
109 mpq_set_z(lp->A[rows-1][1+k], obj[k]);
110 dd_neg(lp->A[rows-1][1+P->Dimension+k], lp->A[rows-1][1+k]);
112 dd_ErrorType err = dd_NoError;
113 dd_LPSolve(lp, dd_DualSimplex, &err);
114 assert(err == dd_NoError);
117 void CDD_LP::get_obj_val(mytype *F)
119 dd_set(*F, lp->optvalue);
120 assert(dd_sgn(*F) >= 0);
123 int CDD_LP::add_row(Value *row, int dim)
125 assert(P->Dimension == dim);
126 Vector_Copy(row, eq->p[neq], dim);
127 return neq++;
130 void CDD_LP::get_alpha(int row, mytype *alpha)
132 row += 2 * P->NbConstraints;
133 /* Value for basic variable is zero, right ?
135 dd_set(*alpha, dd_zero);
136 for (int j = 1; j < lp->d; j++) {
137 if (lp->nbindex[j+1] == row+1) {
138 dd_set(*alpha, lp->dsol[j]);
139 break;
140 } else if (lp->nbindex[j+1] == row+neq+1) {
141 dd_neg(*alpha, lp->dsol[j]);
142 break;
147 void CDD_LP::del_row()
149 assert(neq > 0);
150 neq--;