basis_reduction_cdd.cc

   1 #define GMPRATIONAL
   2 #include <setoper.h>
   3 #include <cdd.h>
   4 #include <barvinok/basis_reduction.h>
   5 #include "initcdd.h"
   6
   7 struct CDD_LP {
   8     Polyhedron *P;
   9     Value      *obj;
  10     dd_LPType  *lp;
  11     Matrix     *eq;
  12     int         neq;
  13
  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);
  25     }
  26 };
  27
  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"
  49
  50 CDD_LP::CDD_LP(Polyhedron *P)
  51 {
  52     INIT_CDD;
  53
  54     this->P = P;
  55     lp = NULL;
  56     eq = Matrix_Alloc(P->Dimension, P->Dimension);
  57     neq = 0;
  58 }
  59
  60 void CDD_LP::set_obj(Value *row, int dim)
  61 {
  62     assert(P->Dimension == dim);
  63     obj = row;
  64 }
  65
  66 void CDD_LP::solve()
  67 {
  68     if (lp)
  69         dd_FreeLPData(lp);
  70
  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]);
  83             }
  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]);
  90                 }
  91                 dd_neg(lp->A[irev][0], lp->A[row][0]);
  92                 ++irev;
  93             }
  94         }
  95     }
  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]);
 104         }
 105         ++irev;
 106     }
 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]);
 111     }
 112     dd_ErrorType err = dd_NoError;
 113     dd_LPSolve(lp, dd_DualSimplex, &err);
 114     assert(err == dd_NoError);
 115 }
 116
 117 void CDD_LP::get_obj_val(mytype *F)
 118 {
 119     dd_set(*F, lp->optvalue);
 120     assert(dd_sgn(*F) >= 0);
 121 }
 122
 123 int CDD_LP::add_row(Value *row, int dim)
 124 {
 125     assert(P->Dimension == dim);
 126     Vector_Copy(row, eq->p[neq], dim);
 127     return neq++;
 128 }
 129
 130 void CDD_LP::get_alpha(int row, mytype *alpha)
 131 {
 132     row += 2 * P->NbConstraints;
 133     /* Value for basic variable is zero, right ?
 134      */
 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+2*P->NbEq+1) {
 141             dd_neg(*alpha, lp->dsol[j]);
 142             break;
 143         }
 144     }
 145 }
 146
 147 void CDD_LP::del_row()
 148 {
 149     assert(neq > 0);
 150     neq--;
 151 }