basis_reduction_pip.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_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     /* We only call this function on a polytope that is known
 112      * to be (rationally) non-empty.
 113      * However, if the test for rational emptiness was inaccurate
 114      * (e.g., when using GLPK with big coefficients), the input
 115      * polytope may be empty.  It's better to just abort here
 116      * so the user can select a different linear solver.
 117      */
 118     assert(lp->sol->list);
 119     return value_zero_p(lp->sol->list->vector->the_deno[0]);
 120 }
 121
 122 static void get_obj_val(struct pip_lp* lp, mpq_t *F)
 123 {
 124     value_assign(mpq_numref(*F), lp->sol->list->vector->the_vector[0]);
 125     value_assign(mpq_denref(*F), lp->sol->list->vector->the_deno[0]);
 126 }
 127
 128 static void delete_lp(struct pip_lp *lp)
 129 {
 130     if (lp->sol)
 131         pip_quast_free(lp->sol);
 132     Matrix_Free(lp->eq);
 133     free(lp);
 134 }
 135
 136 static int add_lp_row(struct pip_lp *lp, Value *row, int dim)
 137 {
 138     assert(lp->P->Dimension == dim);
 139     Vector_Copy(row, lp->eq->p[lp->neq], dim);
 140     return lp->neq++;
 141 }
 142
 143 static void get_alpha(struct pip_lp* lp, int row, mpq_t *alpha)
 144 {
 145     PipList *list;
 146     int i = 0;
 147
 148     assert(lp->sol->next_then);
 149     assert(lp->sol->next_then->list);
 150
 151     for (i = 0, list = lp->sol->next_then->list; i < row; ++i, list = list->next)
 152         ;
 153
 154     value_assign(mpq_numref(*alpha), list->vector->the_vector[0]);
 155     value_assign(mpq_denref(*alpha), list->vector->the_deno[0]);
 156 }
 157
 158 static void del_lp_row(struct pip_lp *lp)
 159 {
 160     assert(lp->neq > 0);
 161     lp->neq--;
 162 }