basis_reduction_pip.c

   1 #include <stdio.h>  /* needed for piplib ! */
   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_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*lp->P->NbConstraints + lp->neq;
  69     int          cols = 2*dim + 3;
  70     int          i, j;
  71
  72     if (lp->sol)
  73         pip_quast_free(lp->sol);
  74
  75     domain = pip_matrix_alloc(rows, cols);
  76     for (i = 0; i < lp->neq; ++i)
  77         for (j = 0; j < dim; ++j) {
  78             value_assign(domain->p[i][1+1+j], lp->eq->p[i][j]);
  79             value_oppose(domain->p[i][1+1+dim+j], lp->eq->p[i][j]);
  80         }
  81     value_set_si(domain->p[lp->neq][1], -1);
  82     for (j = 0; j < dim; ++j) {
  83         value_assign(domain->p[lp->neq][1+1+j], lp->obj[j]);
  84         value_oppose(domain->p[lp->neq][1+1+dim+j], lp->obj[j]);
  85     }
  86     for (i = 0; i < lp->P->NbConstraints; ++i) {
  87         value_assign(domain->p[lp->neq+1+i][0], lp->P->Constraint[i][0]);
  88         value_assign(domain->p[lp->neq+1+lp->P->NbConstraints+i][0],
  89                      lp->P->Constraint[i][0]);
  90         for (j = 0; j < dim; ++j) {
  91             value_assign(domain->p[lp->neq+1+i][2+j], lp->P->Constraint[i][1+j]);
  92             value_assign(domain->p[lp->neq+1+lp->P->NbConstraints+i][2+dim+j],
  93                          lp->P->Constraint[i][1+j]);
  94         }
  95         value_assign(domain->p[lp->neq+1+i][cols-1], lp->P->Constraint[i][1+dim]);
  96         value_assign(domain->p[lp->neq+1+lp->P->NbConstraints+i][cols-1],
  97                      lp->P->Constraint[i][1+dim]);
  98     }
  99
 100     options = pip_options_init();
 101     options->Urs_unknowns = -1;
 102     options->Maximize = 1;
 103     options->Nq = 0;
 104     if (lp->neq)
 105         options->Compute_dual = 1;
 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     assert(lp->sol->next_then);
 141     assert(lp->sol->next_then->list);
 142
 143     for (i = 0, list = lp->sol->next_then->list; i < row; ++i, list = list->next)
 144         ;
 145
 146     value_assign(mpq_numref(*alpha), list->vector->the_vector[0]);
 147     value_assign(mpq_denref(*alpha), list->vector->the_deno[0]);
 148 }
 149
 150 static void del_lp_row(struct pip_lp *lp)
 151 {
 152     assert(lp->neq > 0);
 153     lp->neq--;
 154 }