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;
44 dim = isl_basic_set_total_dim(bset);
45 basis = isl_mat_identity(bset->ctx, dim);
64 b_tmp = isl_vec_alloc(bset->ctx, dim);
68 F = isl_alloc_array(bset->ctx, GBR_type, dim);
69 alpha_buffer[0] = isl_alloc_array(bset->ctx, GBR_type, dim);
70 alpha_buffer[1] = isl_alloc_array(bset->ctx, GBR_type, dim);
71 alpha_saved = alpha_buffer[0];
73 if (!F || !alpha_buffer[0] || !alpha_buffer[1])
76 for (i = 0; i < dim; ++i) {
78 GBR_init(alpha_buffer[0][i]);
79 GBR_init(alpha_buffer[1][i]);
84 lp = GBR_lp_init(bset);
90 GBR_lp_set_obj(lp, basis->row[0], dim);
91 bset->ctx->stats->gbr_solved_lps++;
92 unbounded = GBR_lp_solve(lp);
93 isl_assert(bset->ctx, !unbounded, goto error);
94 GBR_lp_get_obj_val(lp, &F[0]);
98 row = GBR_lp_next_row(lp);
99 GBR_set(F_new, F_saved);
100 GBR_set(alpha, alpha_saved[i]);
102 row = GBR_lp_add_row(lp, basis->row[i], dim);
103 GBR_lp_set_obj(lp, basis->row[i+1], dim);
104 bset->ctx->stats->gbr_solved_lps++;
105 unbounded = GBR_lp_solve(lp);
106 isl_assert(bset->ctx, !unbounded, goto error);
107 GBR_lp_get_obj_val(lp, &F_new);
109 GBR_lp_get_alpha(lp, row, &alpha);
112 save_alpha(lp, row-i, i, alpha_saved);
116 GBR_set(F[i+1], F_new);
118 GBR_floor(mu[0], alpha);
119 GBR_ceil(mu[1], alpha);
121 if (isl_int_eq(mu[0], mu[1]))
122 isl_int_set(tmp, mu[0]);
126 for (j = 0; j <= 1; ++j) {
127 isl_int_set(tmp, mu[j]);
128 isl_seq_combine(b_tmp->el,
129 bset->ctx->one, basis->row[i+1],
130 tmp, basis->row[i], dim);
131 GBR_lp_set_obj(lp, b_tmp->el, 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, &mu_F[j]);
137 save_alpha(lp, row-i, i, alpha_buffer[j]);
140 if (GBR_lt(mu_F[0], mu_F[1]))
145 isl_int_set(tmp, mu[j]);
146 GBR_set(F_new, mu_F[j]);
147 alpha_saved = alpha_buffer[j];
149 isl_seq_combine(basis->row[i+1],
150 bset->ctx->one, basis->row[i+1],
151 tmp, basis->row[i], dim);
153 GBR_set(F_old, F[i]);
156 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
157 GBR_set_ui(mu_F[0], 4);
158 GBR_mul(mu_F[0], mu_F[0], F_new);
159 GBR_set_ui(mu_F[1], 3);
160 GBR_mul(mu_F[1], mu_F[1], F_old);
161 if (GBR_lt(mu_F[0], mu_F[1])) {
162 basis = isl_mat_swap_rows(basis, i, i + 1);
165 GBR_set(F_saved, F_new);
169 GBR_set(F[0], F_new);
170 if (bset->ctx->gbr_only_first &&
175 GBR_lp_add_row(lp, basis->row[i], dim);
189 for (i = 0; i < dim; ++i) {
191 GBR_clear(alpha_buffer[0][i]);
192 GBR_clear(alpha_buffer[1][i]);
195 free(alpha_buffer[0]);
196 free(alpha_buffer[1]);
209 isl_int_clear(mu[0]);
210 isl_int_clear(mu[1]);