basis_reduction_pip_dual.c

   1 #include <assert.h>
   2 #include <piplib/piplibMP.h>
   3 #include <barvinok/basis_reduction.h>
   4
   5 struct pip_lp {
   6     Polyhedron *P;
   7     Value      *obj;
   8     Matrix     *eq;
   9     int         neq;
  10     PipQuast   *sol;
  11 };
  12
  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);
  21
  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_dual_Polyhedron_Reduced_Basis
  42 #include "basis_reduction_templ.c"
  43
  44 #define ALLOC(type) (type*)malloc(sizeof(type))
  45
  46 static struct pip_lp *init_lp(Polyhedron *P)
  47 {
  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;
  53
  54     return lp;
  55 }
  56
  57 static void set_lp_obj(struct pip_lp *lp, Value *row, int dim)
  58 {
  59     assert(lp->P->Dimension == dim);
  60     lp->obj = row;
  61 }
  62
  63 static int solve_lp(struct pip_lp *lp)
  64 {
  65     PipOptions  *options;
  66     PipMatrix   *domain;
  67     unsigned     dim = lp->P->Dimension;
  68     int          rows = 1 + 2*dim + 2*lp->P->NbConstraints;
  69     int          cols = 1 + 1 + lp->neq + 2*lp->P->NbConstraints + 1;
  70     int          i, j;
  71
  72     if (lp->sol)
  73         pip_quast_free(lp->sol);
  74
  75     assert(lp->P->NbEq == 0);
  76     domain = pip_matrix_alloc(rows, cols);
  77     for (i = 0; i < lp->neq; ++i)
  78         for (j = 0; j < dim; ++j) {
  79             value_assign(domain->p[j][1+1+i], lp->eq->p[i][j]);
  80             value_oppose(domain->p[dim+j][1+1+i], lp->eq->p[i][j]);
  81         }
  82     value_set_si(domain->p[2*dim][1], 1);
  83     for (j = 0; j < dim; ++j) {
  84         value_oppose(domain->p[j][cols-1], lp->obj[j]);
  85         value_assign(domain->p[dim+j][cols-1], lp->obj[j]);
  86     }
  87     for (i = 0; i < lp->P->NbConstraints; ++i) {
  88         for (j = 0; j < dim; ++j) {
  89             value_assign(domain->p[j][1+1+lp->neq+i], lp->P->Constraint[i][1+j]);
  90             value_assign(domain->p[dim+j][1+1+lp->neq+lp->P->NbConstraints+i],
  91                          lp->P->Constraint[i][1+j]);
  92         }
  93         value_assign(domain->p[2*dim][1+1+lp->neq+i], lp->P->Constraint[i][1+dim]);
  94         value_assign(domain->p[2*dim][1+1+lp->neq+lp->P->NbConstraints+i],
  95                      lp->P->Constraint[i][1+dim]);
  96         value_set_si(domain->p[2*dim+1+i][0], 1);
  97         value_set_si(domain->p[2*dim+1+i][1+1+lp->neq+i], -1);
  98         value_set_si(domain->p[2*dim+1+lp->P->NbConstraints+i][0], 1);
  99         value_set_si(domain->p[2*dim+1+lp->P->NbConstraints+i]
 100                                 [1+1+lp->neq+lp->P->NbConstraints+i], -1);
 101     }
 102
 103     options = pip_options_init();
 104     options->Urs_unknowns = -1;
 105     options->Nq = 0;
 106     lp->sol = pip_solve(domain, NULL, -1, options);
 107     pip_matrix_free(domain);
 108
 109     pip_options_free(options);
 110
 111     return value_zero_p(lp->sol->list->vector->the_deno[0]);
 112 }
 113
 114 static void get_obj_val(struct pip_lp* lp, mpq_t *F)
 115 {
 116     value_assign(mpq_numref(*F), lp->sol->list->vector->the_vector[0]);
 117     value_assign(mpq_denref(*F), lp->sol->list->vector->the_deno[0]);
 118 }
 119
 120 static void delete_lp(struct pip_lp *lp)
 121 {
 122     if (lp->sol)
 123         pip_quast_free(lp->sol);
 124     Matrix_Free(lp->eq);
 125     free(lp);
 126 }
 127
 128 static int add_lp_row(struct pip_lp *lp, Value *row, int dim)
 129 {
 130     assert(lp->P->Dimension == dim);
 131     Vector_Copy(row, lp->eq->p[lp->neq], dim);
 132     return lp->neq++;
 133 }
 134
 135 static void get_alpha(struct pip_lp* lp, int row, mpq_t *alpha)
 136 {
 137     PipList *list;
 138     int i = 0;
 139
 140     for (i = 0, list = lp->sol->list->next; i < row; ++i, list = list->next)
 141         ;
 142
 143     if (value_zero_p(list->vector->the_deno[0])) {
 144         value_set_si(mpq_numref(*alpha), 0);
 145         value_set_si(mpq_denref(*alpha), 1);
 146     } else {
 147         value_oppose(mpq_numref(*alpha), list->vector->the_vector[0]);
 148         value_assign(mpq_denref(*alpha), list->vector->the_deno[0]);
 149     }
 150 }
 151
 152 static void del_lp_row(struct pip_lp *lp)
 153 {
 154     assert(lp->neq > 0);
 155     lp->neq--;
 156 }