isl_pw_qpolynomial_sum: handle existentials in wrapped domains
[barvinok.git] / basis_reduction_templ.c
blob3ad26101099d9e9a8960e24fc7d9067fe338ef02
1 #include <stdlib.h>
2 #include <barvinok/basis_reduction.h>
3 #include <barvinok/options.h>
5 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
7 static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
9 int i;
11 for (i = 0; i < n; ++i)
12 GBR_lp_get_alpha(lp, first+i, &alpha[i]);
15 /* This function implements the algorithm described in
16 * "An Implementation of the Generalized Basis Reduction Algorithm
17 * for Integer Programming" of Cook el al. to compute a reduced basis.
18 * We use \epsilon = 1/4.
20 * If options->gbr_only_first is set, the user is only interested
21 * in the first direction. In this case we stop the basis reduction when
22 * - the width in the first direction becomes smaller than 2
23 * or
24 * - we have moved forward all the way to the last direction
25 * and then back again all the way to the first direction.
27 Matrix *Polyhedron_Reduced_Basis(Polyhedron *P,
28 struct barvinok_options *options)
30 int dim = P->Dimension;
31 int i;
32 Matrix *basis = Identity(dim);
33 GBR_LP *lp;
34 GBR_type F_old, alpha, F_new;
35 int row;
36 Value one, tmp;
37 Vector *b_tmp;
38 GBR_type *F;
39 GBR_type *alpha_buffer[2];
40 GBR_type *alpha_saved;
41 GBR_type F_saved;
42 int use_saved = 0;
43 Value mu[2];
44 GBR_type mu_F[2];
45 GBR_type two;
47 if (P->Dimension == 1)
48 return basis;
50 value_init(one);
51 value_init(tmp);
52 value_set_si(one, 1);
53 value_init(mu[0]);
54 value_init(mu[1]);
56 b_tmp = Vector_Alloc(dim);
58 F = ALLOCN(GBR_type, dim);
59 alpha_buffer[0] = ALLOCN(GBR_type, dim);
60 alpha_buffer[1] = ALLOCN(GBR_type, dim);
61 alpha_saved = alpha_buffer[0];
63 for (i = 0; i < dim; ++i) {
64 GBR_init(F[i]);
65 GBR_init(alpha_buffer[0][i]);
66 GBR_init(alpha_buffer[1][i]);
68 GBR_init(alpha);
69 GBR_init(F_old);
70 GBR_init(F_new);
71 GBR_init(F_saved);
72 GBR_init(mu_F[0]);
73 GBR_init(mu_F[1]);
74 GBR_init(two);
76 GBR_set_ui(two, 2);
78 lp = GBR_lp_init(P);
80 i = 0;
82 GBR_lp_set_obj(lp, basis->p[0], dim);
83 options->stats->gbr_solved_lps++;
84 if (GBR_lp_solve(lp))
85 goto unbounded;
86 GBR_lp_get_obj_val(lp, &F[0]);
88 do {
89 if (use_saved) {
90 row = GBR_lp_next_row(lp);
91 GBR_set(F_new, F_saved);
92 GBR_set(alpha, alpha_saved[i]);
93 } else {
94 row = GBR_lp_add_row(lp, basis->p[i], dim);
95 GBR_lp_set_obj(lp, basis->p[i+1], dim);
96 options->stats->gbr_solved_lps++;
97 if (GBR_lp_solve(lp))
98 goto unbounded;
99 GBR_lp_get_obj_val(lp, &F_new);
101 GBR_lp_get_alpha(lp, row, &alpha);
103 if (i > 0)
104 save_alpha(lp, row-i, i, alpha_saved);
106 GBR_lp_del_row(lp);
108 GBR_set(F[i+1], F_new);
110 GBR_floor(mu[0], alpha);
111 GBR_ceil(mu[1], alpha);
113 if (value_eq(mu[0], mu[1]))
114 value_assign(tmp, mu[0]);
115 else {
116 int j;
118 for (j = 0; j <= 1; ++j) {
119 value_assign(tmp, mu[j]);
120 Vector_Combine(basis->p[i+1], basis->p[i], b_tmp->p, one, tmp, dim);
121 GBR_lp_set_obj(lp, b_tmp->p, dim);
122 options->stats->gbr_solved_lps++;
123 if (GBR_lp_solve(lp))
124 goto unbounded;
125 GBR_lp_get_obj_val(lp, &mu_F[j]);
126 if (i > 0)
127 save_alpha(lp, row-i, i, alpha_buffer[j]);
130 if (GBR_lt(mu_F[0], mu_F[1]))
131 j = 0;
132 else
133 j = 1;
135 value_assign(tmp, mu[j]);
136 GBR_set(F_new, mu_F[j]);
137 alpha_saved = alpha_buffer[j];
139 Vector_Combine(basis->p[i+1], basis->p[i], basis->p[i+1], one, tmp, dim);
141 GBR_set(F_old, F[i]);
143 use_saved = 0;
144 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
145 GBR_set_ui(mu_F[0], 4);
146 GBR_mul(mu_F[0], mu_F[0], F_new);
147 GBR_set_ui(mu_F[1], 3);
148 GBR_mul(mu_F[1], mu_F[1], F_old);
149 if (GBR_lt(mu_F[0], mu_F[1])) {
150 Vector_Exchange(basis->p[i], basis->p[i+1], dim);
151 if (i > 0) {
152 use_saved = 1;
153 GBR_set(F_saved, F_new);
154 GBR_lp_del_row(lp);
155 --i;
156 } else {
157 GBR_set(F[0], F_new);
158 if (options->gbr_only_first && GBR_lt(F[0], two))
159 break;
161 } else {
162 GBR_lp_add_row(lp, basis->p[i], dim);
163 ++i;
165 } while (i < dim-1);
167 if (0) {
168 unbounded:
169 Matrix_Free(basis);
170 basis = NULL;
172 Vector_Free(b_tmp);
174 value_clear(one);
175 value_clear(tmp);
176 value_clear(mu[0]);
177 value_clear(mu[1]);
178 for (i = 0; i < dim; ++i) {
179 GBR_clear(F[i]);
180 GBR_clear(alpha_buffer[0][i]);
181 GBR_clear(alpha_buffer[1][i]);
183 free(F);
184 free(alpha_buffer[0]);
185 free(alpha_buffer[1]);
186 GBR_clear(alpha);
187 GBR_clear(F_old);
188 GBR_clear(F_new);
189 GBR_clear(F_saved);
190 GBR_clear(mu_F[0]);
191 GBR_clear(mu_F[1]);
192 GBR_clear(two);
194 GBR_lp_delete(lp);
196 return basis;