basis_reduction_cdd.cc

   1 #include <assert.h>
   2 #define GMPRATIONAL
   3 #include <cddlib/setoper.h>
   4 #include <cddlib/cdd.h>
   5 #include <barvinok/basis_reduction.h>
   6 #include "initcdd.h"
   7
   8 struct CDD_LP {
   9     Polyhedron *P;
  10     Value      *obj;
  11     dd_LPType  *lp;
  12     Matrix     *eq;
  13     int         neq;
  14
  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);
  26     }
  27 };
  28
  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"
  50
  51 CDD_LP::CDD_LP(Polyhedron *P)
  52 {
  53     INIT_CDD;
  54
  55     this->P = P;
  56     lp = NULL;
  57     eq = Matrix_Alloc(P->Dimension, P->Dimension);
  58     neq = 0;
  59 }
  60
  61 void CDD_LP::set_obj(Value *row, int dim)
  62 {
  63     assert(P->Dimension == dim);
  64     obj = row;
  65 }
  66
  67 int CDD_LP::solve()
  68 {
  69     if (lp)
  70         dd_FreeLPData(lp);
  71
  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]);
  84             }
  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]);
  91                 }
  92                 dd_neg(lp->A[irev][0], lp->A[row][0]);
  93                 ++irev;
  94             }
  95         }
  96     }
  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]);
 105         }
 106         ++irev;
 107     }
 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]);
 112     }
 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.
 121      */
 122     assert(lp->LPS != dd_Inconsistent);
 123
 124     return lp->LPS == dd_DualInconsistent;
 125 }
 126
 127 void CDD_LP::get_obj_val(mytype *F)
 128 {
 129     dd_set(*F, lp->optvalue);
 130     assert(dd_sgn(*F) >= 0);
 131 }
 132
 133 int CDD_LP::add_row(Value *row, int dim)
 134 {
 135     assert(P->Dimension == dim);
 136     Vector_Copy(row, eq->p[neq], dim);
 137     return neq++;
 138 }
 139
 140 void CDD_LP::get_alpha(int row, mytype *alpha)
 141 {
 142     row += 2 * P->NbConstraints;
 143     /* Value for basic variable is zero, right ?
 144      */
 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;
 153         }
 154     }
 155 }
 156
 157 void CDD_LP::del_row()
 158 {
 159     assert(neq > 0);
 160     neq--;
 161 }