2 #include "isl_basis_reduction.h"
4 static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha)
8 for (i = 0; i < n; ++i)
9 GBR_lp_get_alpha(lp, first + i, &alpha[i]);
12 /* This function implements the algorithm described in
13 * "An Implementation of the Generalized Basis Reduction Algorithm
14 * for Integer Programming" of Cook el al. to compute a reduced basis.
15 * We use \epsilon = 1/4.
17 * If ctx->gbr_only_first is set, the user is only interested
18 * in the first direction. In this case we stop the basis reduction when
19 * the width in the first direction becomes smaller than 2.
21 struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset)
24 struct isl_mat *basis;
28 GBR_type F_old, alpha, F_new;
31 struct isl_vec *b_tmp;
33 GBR_type *alpha_buffer[2] = { NULL, NULL };
34 GBR_type *alpha_saved;
50 dim = isl_basic_set_total_dim(bset);
51 basis = isl_mat_identity(bset->ctx, dim);
71 b_tmp = isl_vec_alloc(bset->ctx, dim);
75 F = isl_alloc_array(bset->ctx, GBR_type, dim);
76 alpha_buffer[0] = isl_alloc_array(bset->ctx, GBR_type, dim);
77 alpha_buffer[1] = isl_alloc_array(bset->ctx, GBR_type, dim);
78 alpha_saved = alpha_buffer[0];
80 if (!F || !alpha_buffer[0] || !alpha_buffer[1])
83 for (i = 0; i < dim; ++i) {
85 GBR_init(alpha_buffer[0][i]);
86 GBR_init(alpha_buffer[1][i]);
92 lp = GBR_lp_init(bset);
98 GBR_lp_set_obj(lp, basis->row[0], dim);
99 bset->ctx->stats->gbr_solved_lps++;
100 unbounded = GBR_lp_solve(lp);
101 isl_assert(bset->ctx, !unbounded, goto error);
102 GBR_lp_get_obj_val(lp, &F[0]);
104 if (GBR_lt(F[0], one)) {
105 if (!GBR_is_zero(F[0])) {
106 empty = GBR_lp_cut(lp, basis->row[0]);
116 GBR_lp_set_obj(lp, basis->row[i+1], dim);
117 bset->ctx->stats->gbr_solved_lps++;
118 unbounded = GBR_lp_solve(lp);
119 isl_assert(bset->ctx, !unbounded, goto error);
120 GBR_lp_get_obj_val(lp, &F_new);
121 fixed = GBR_lp_is_fixed(lp);
122 GBR_set_ui(alpha, 0);
125 row = GBR_lp_next_row(lp);
126 GBR_set(F_new, F_saved);
128 GBR_set(alpha, alpha_saved[i]);
130 row = GBR_lp_add_row(lp, basis->row[i], dim);
131 GBR_lp_set_obj(lp, basis->row[i+1], dim);
132 bset->ctx->stats->gbr_solved_lps++;
133 unbounded = GBR_lp_solve(lp);
134 isl_assert(bset->ctx, !unbounded, goto error);
135 GBR_lp_get_obj_val(lp, &F_new);
136 fixed = GBR_lp_is_fixed(lp);
138 GBR_lp_get_alpha(lp, row, &alpha);
141 save_alpha(lp, row-i, i, alpha_saved);
145 GBR_set(F[i+1], F_new);
147 GBR_floor(mu[0], alpha);
148 GBR_ceil(mu[1], alpha);
150 if (isl_int_eq(mu[0], mu[1]))
151 isl_int_set(tmp, mu[0]);
155 for (j = 0; j <= 1; ++j) {
156 isl_int_set(tmp, mu[j]);
157 isl_seq_combine(b_tmp->el,
158 bset->ctx->one, basis->row[i+1],
159 tmp, basis->row[i], dim);
160 GBR_lp_set_obj(lp, b_tmp->el, dim);
161 bset->ctx->stats->gbr_solved_lps++;
162 unbounded = GBR_lp_solve(lp);
163 isl_assert(bset->ctx, !unbounded, goto error);
164 GBR_lp_get_obj_val(lp, &mu_F[j]);
165 mu_fixed[j] = GBR_lp_is_fixed(lp);
167 save_alpha(lp, row-i, i, alpha_buffer[j]);
170 if (GBR_lt(mu_F[0], mu_F[1]))
175 isl_int_set(tmp, mu[j]);
176 GBR_set(F_new, mu_F[j]);
178 alpha_saved = alpha_buffer[j];
180 isl_seq_combine(basis->row[i+1],
181 bset->ctx->one, basis->row[i+1],
182 tmp, basis->row[i], dim);
184 if (i+1 == n_zero && fixed) {
185 if (!GBR_is_zero(F[i+1])) {
186 empty = GBR_lp_cut(lp, basis->row[i+1]);
189 GBR_set_ui(F[i+1], 0);
194 GBR_set(F_old, F[i]);
197 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
198 GBR_set_ui(mu_F[0], 4);
199 GBR_mul(mu_F[0], mu_F[0], F_new);
200 GBR_set_ui(mu_F[1], 3);
201 GBR_mul(mu_F[1], mu_F[1], F_old);
202 if (GBR_lt(mu_F[0], mu_F[1])) {
203 basis = isl_mat_swap_rows(basis, i, i + 1);
206 GBR_set(F_saved, F_new);
211 GBR_set(F[0], F_new);
212 if (bset->ctx->gbr_only_first &&
217 if (!GBR_is_zero(F[0])) {
218 empty = GBR_lp_cut(lp, basis->row[0]);
227 GBR_lp_add_row(lp, basis->row[i], dim);
244 for (i = 0; i < dim; ++i) {
246 GBR_clear(alpha_buffer[0][i]);
247 GBR_clear(alpha_buffer[1][i]);
250 free(alpha_buffer[0]);
251 free(alpha_buffer[1]);
265 isl_int_clear(mu[0]);
266 isl_int_clear(mu[1]);