iscc: add operations on int objects
[barvinok.git] / basis_reduction_pip_dual.c
blob2c403f8ef012a905e0e26e976b8728245b47fb2d
1 #include <assert.h>
2 #include <piplib/piplibMP.h>
3 #include <barvinok/basis_reduction.h>
5 struct pip_lp {
6 Polyhedron *P;
7 Value *obj;
8 Matrix *eq;
9 int neq;
10 PipQuast *sol;
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);
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"
44 #define ALLOC(type) (type*)malloc(sizeof(type))
46 static struct pip_lp *init_lp(Polyhedron *P)
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;
54 return lp;
57 static void set_lp_obj(struct pip_lp *lp, Value *row, int dim)
59 assert(lp->P->Dimension == dim);
60 lp->obj = row;
63 static int solve_lp(struct pip_lp *lp)
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;
72 if (lp->sol)
73 pip_quast_free(lp->sol);
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]);
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]);
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]);
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);
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);
109 pip_options_free(options);
111 return value_zero_p(lp->sol->list->vector->the_deno[0]);
114 static void get_obj_val(struct pip_lp* lp, mpq_t *F)
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]);
120 static void delete_lp(struct pip_lp *lp)
122 if (lp->sol)
123 pip_quast_free(lp->sol);
124 Matrix_Free(lp->eq);
125 free(lp);
128 static int add_lp_row(struct pip_lp *lp, Value *row, int dim)
130 assert(lp->P->Dimension == dim);
131 Vector_Copy(row, lp->eq->p[lp->neq], dim);
132 return lp->neq++;
135 static void get_alpha(struct pip_lp* lp, int row, mpq_t *alpha)
137 PipList *list;
138 int i = 0;
140 for (i = 0, list = lp->sol->list->next; i < row; ++i, list = list->next)
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]);
152 static void del_lp_row(struct pip_lp *lp)
154 assert(lp->neq > 0);
155 lp->neq--;