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 /* Compute a reduced basis for the set represented by the tableau "tab".
13 * tab->basis, must be initialized by the calling function to an affine
14 * unimodular basis, is updated to reflect the reduced basis.
15 * The first tab->n_zero rows of the basis (ignoring the constant row)
16 * are assumed to correspond to equalities and are left untouched.
17 * tab->n_zero is updated to reflect any additional equalities that
18 * have been detected in the first rows of the new basis.
19 * The final tab->n_unbounded rows of the basis are assumed to correspond
20 * to unbounded directions and are also left untouched.
21 * In particular this means that the remaining rows are assumed to
22 * correspond to bounded directions.
24 * This function implements the algorithm described in
25 * "An Implementation of the Generalized Basis Reduction Algorithm
26 * for Integer Programming" of Cook el al. to compute a reduced basis.
27 * We use \epsilon = 1/4.
29 * If ctx->gbr_only_first is set, the user is only interested
30 * in the first direction. In this case we stop the basis reduction when
31 * the width in the first direction becomes smaller than 2.
33 struct isl_tab
*isl_tab_compute_reduced_basis(struct isl_tab
*tab
)
41 GBR_type F_old
, alpha
, F_new
;
44 struct isl_vec
*b_tmp
;
46 GBR_type
*alpha_buffer
[2] = { NULL
, NULL
};
47 GBR_type
*alpha_saved
;
72 n_bounded
= dim
- tab
->n_unbounded
;
73 if (n_bounded
<= tab
->n_zero
+ 1)
89 b_tmp
= isl_vec_alloc(ctx
, dim
);
93 F
= isl_alloc_array(ctx
, GBR_type
, n_bounded
);
94 alpha_buffer
[0] = isl_alloc_array(ctx
, GBR_type
, n_bounded
);
95 alpha_buffer
[1] = isl_alloc_array(ctx
, GBR_type
, n_bounded
);
96 alpha_saved
= alpha_buffer
[0];
98 if (!F
|| !alpha_buffer
[0] || !alpha_buffer
[1])
101 for (i
= 0; i
< n_bounded
; ++i
) {
103 GBR_init(alpha_buffer
[0][i
]);
104 GBR_init(alpha_buffer
[1][i
]);
110 lp
= GBR_lp_init(tab
);
116 GBR_lp_set_obj(lp
, B
->row
[1+i
]+1, dim
);
117 ctx
->stats
->gbr_solved_lps
++;
118 unbounded
= GBR_lp_solve(lp
);
119 isl_assert(ctx
, !unbounded
, goto error
);
120 GBR_lp_get_obj_val(lp
, &F
[i
]);
122 if (GBR_lt(F
[i
], one
)) {
123 if (!GBR_is_zero(F
[i
])) {
124 empty
= GBR_lp_cut(lp
, B
->row
[1+i
]+1);
133 if (i
+1 == tab
->n_zero
) {
134 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
135 ctx
->stats
->gbr_solved_lps
++;
136 unbounded
= GBR_lp_solve(lp
);
137 isl_assert(ctx
, !unbounded
, goto error
);
138 GBR_lp_get_obj_val(lp
, &F_new
);
139 fixed
= GBR_lp_is_fixed(lp
);
140 GBR_set_ui(alpha
, 0);
143 row
= GBR_lp_next_row(lp
);
144 GBR_set(F_new
, F_saved
);
146 GBR_set(alpha
, alpha_saved
[i
]);
148 row
= GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
149 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
150 ctx
->stats
->gbr_solved_lps
++;
151 unbounded
= GBR_lp_solve(lp
);
152 isl_assert(ctx
, !unbounded
, goto error
);
153 GBR_lp_get_obj_val(lp
, &F_new
);
154 fixed
= GBR_lp_is_fixed(lp
);
156 GBR_lp_get_alpha(lp
, row
, &alpha
);
159 save_alpha(lp
, row
-i
, i
, alpha_saved
);
161 if (GBR_lp_del_row(lp
) < 0)
164 GBR_set(F
[i
+1], F_new
);
166 GBR_floor(mu
[0], alpha
);
167 GBR_ceil(mu
[1], alpha
);
169 if (isl_int_eq(mu
[0], mu
[1]))
170 isl_int_set(tmp
, mu
[0]);
174 for (j
= 0; j
<= 1; ++j
) {
175 isl_int_set(tmp
, mu
[j
]);
176 isl_seq_combine(b_tmp
->el
,
177 ctx
->one
, B
->row
[1+i
+1]+1,
178 tmp
, B
->row
[1+i
]+1, dim
);
179 GBR_lp_set_obj(lp
, b_tmp
->el
, dim
);
180 ctx
->stats
->gbr_solved_lps
++;
181 unbounded
= GBR_lp_solve(lp
);
182 isl_assert(ctx
, !unbounded
, goto error
);
183 GBR_lp_get_obj_val(lp
, &mu_F
[j
]);
184 mu_fixed
[j
] = GBR_lp_is_fixed(lp
);
186 save_alpha(lp
, row
-i
, i
, alpha_buffer
[j
]);
189 if (GBR_lt(mu_F
[0], mu_F
[1]))
194 isl_int_set(tmp
, mu
[j
]);
195 GBR_set(F_new
, mu_F
[j
]);
197 alpha_saved
= alpha_buffer
[j
];
199 isl_seq_combine(B
->row
[1+i
+1]+1, ctx
->one
, B
->row
[1+i
+1]+1,
200 tmp
, B
->row
[1+i
]+1, dim
);
202 if (i
+1 == tab
->n_zero
&& fixed
) {
203 if (!GBR_is_zero(F
[i
+1])) {
204 empty
= GBR_lp_cut(lp
, B
->row
[1+i
+1]+1);
207 GBR_set_ui(F
[i
+1], 0);
212 GBR_set(F_old
, F
[i
]);
215 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
216 GBR_set_ui(mu_F
[0], 4);
217 GBR_mul(mu_F
[0], mu_F
[0], F_new
);
218 GBR_set_ui(mu_F
[1], 3);
219 GBR_mul(mu_F
[1], mu_F
[1], F_old
);
220 if (GBR_lt(mu_F
[0], mu_F
[1])) {
221 B
= isl_mat_swap_rows(B
, 1 + i
, 1 + i
+ 1);
222 if (i
> tab
->n_zero
) {
224 GBR_set(F_saved
, F_new
);
226 if (GBR_lp_del_row(lp
) < 0)
230 GBR_set(F
[tab
->n_zero
], F_new
);
231 if (ctx
->gbr_only_first
&& GBR_lt(F
[tab
->n_zero
], two
))
235 if (!GBR_is_zero(F
[tab
->n_zero
])) {
236 empty
= GBR_lp_cut(lp
, B
->row
[1+tab
->n_zero
]+1);
239 GBR_set_ui(F
[tab
->n_zero
], 0);
245 GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
248 } while (i
< n_bounded
- 1);
262 for (i
= 0; i
< n_bounded
; ++i
) {
264 GBR_clear(alpha_buffer
[0][i
]);
265 GBR_clear(alpha_buffer
[1][i
]);
268 free(alpha_buffer
[0]);
269 free(alpha_buffer
[1]);
283 isl_int_clear(mu
[0]);
284 isl_int_clear(mu
[1]);
291 struct isl_mat
*isl_basic_set_reduced_basis(struct isl_basic_set
*bset
)
293 struct isl_mat
*basis
;
296 isl_assert(bset
->ctx
, bset
->n_eq
== 0, return NULL
);
298 tab
= isl_tab_from_basic_set(bset
);
299 tab
->basis
= isl_mat_identity(bset
->ctx
, 1 + tab
->n_var
);
300 tab
= isl_tab_compute_reduced_basis(tab
);
304 basis
= isl_mat_copy(tab
->basis
);