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
;
69 n_bounded
= dim
- tab
->n_unbounded
;
70 if (n_bounded
<= tab
->n_zero
+ 1)
86 b_tmp
= isl_vec_alloc(ctx
, dim
);
90 F
= isl_alloc_array(ctx
, GBR_type
, n_bounded
);
91 alpha_buffer
[0] = isl_alloc_array(ctx
, GBR_type
, n_bounded
);
92 alpha_buffer
[1] = isl_alloc_array(ctx
, GBR_type
, n_bounded
);
93 alpha_saved
= alpha_buffer
[0];
95 if (!F
|| !alpha_buffer
[0] || !alpha_buffer
[1])
98 for (i
= 0; i
< n_bounded
; ++i
) {
100 GBR_init(alpha_buffer
[0][i
]);
101 GBR_init(alpha_buffer
[1][i
]);
107 lp
= GBR_lp_init(tab
);
113 GBR_lp_set_obj(lp
, B
->row
[1+i
]+1, dim
);
114 ctx
->stats
->gbr_solved_lps
++;
115 unbounded
= GBR_lp_solve(lp
);
116 isl_assert(ctx
, !unbounded
, goto error
);
117 GBR_lp_get_obj_val(lp
, &F
[i
]);
119 if (GBR_lt(F
[i
], one
)) {
120 if (!GBR_is_zero(F
[i
])) {
121 empty
= GBR_lp_cut(lp
, B
->row
[1+i
]+1);
130 if (i
+1 == tab
->n_zero
) {
131 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
132 ctx
->stats
->gbr_solved_lps
++;
133 unbounded
= GBR_lp_solve(lp
);
134 isl_assert(ctx
, !unbounded
, goto error
);
135 GBR_lp_get_obj_val(lp
, &F_new
);
136 fixed
= GBR_lp_is_fixed(lp
);
137 GBR_set_ui(alpha
, 0);
140 row
= GBR_lp_next_row(lp
);
141 GBR_set(F_new
, F_saved
);
143 GBR_set(alpha
, alpha_saved
[i
]);
145 row
= GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
146 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
147 ctx
->stats
->gbr_solved_lps
++;
148 unbounded
= GBR_lp_solve(lp
);
149 isl_assert(ctx
, !unbounded
, goto error
);
150 GBR_lp_get_obj_val(lp
, &F_new
);
151 fixed
= GBR_lp_is_fixed(lp
);
153 GBR_lp_get_alpha(lp
, row
, &alpha
);
156 save_alpha(lp
, row
-i
, i
, alpha_saved
);
160 GBR_set(F
[i
+1], F_new
);
162 GBR_floor(mu
[0], alpha
);
163 GBR_ceil(mu
[1], alpha
);
165 if (isl_int_eq(mu
[0], mu
[1]))
166 isl_int_set(tmp
, mu
[0]);
170 for (j
= 0; j
<= 1; ++j
) {
171 isl_int_set(tmp
, mu
[j
]);
172 isl_seq_combine(b_tmp
->el
,
173 ctx
->one
, B
->row
[1+i
+1]+1,
174 tmp
, B
->row
[1+i
]+1, dim
);
175 GBR_lp_set_obj(lp
, b_tmp
->el
, dim
);
176 ctx
->stats
->gbr_solved_lps
++;
177 unbounded
= GBR_lp_solve(lp
);
178 isl_assert(ctx
, !unbounded
, goto error
);
179 GBR_lp_get_obj_val(lp
, &mu_F
[j
]);
180 mu_fixed
[j
] = GBR_lp_is_fixed(lp
);
182 save_alpha(lp
, row
-i
, i
, alpha_buffer
[j
]);
185 if (GBR_lt(mu_F
[0], mu_F
[1]))
190 isl_int_set(tmp
, mu
[j
]);
191 GBR_set(F_new
, mu_F
[j
]);
193 alpha_saved
= alpha_buffer
[j
];
195 isl_seq_combine(B
->row
[1+i
+1]+1, ctx
->one
, B
->row
[1+i
+1]+1,
196 tmp
, B
->row
[1+i
]+1, dim
);
198 if (i
+1 == tab
->n_zero
&& fixed
) {
199 if (!GBR_is_zero(F
[i
+1])) {
200 empty
= GBR_lp_cut(lp
, B
->row
[1+i
+1]+1);
203 GBR_set_ui(F
[i
+1], 0);
208 GBR_set(F_old
, F
[i
]);
211 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
212 GBR_set_ui(mu_F
[0], 4);
213 GBR_mul(mu_F
[0], mu_F
[0], F_new
);
214 GBR_set_ui(mu_F
[1], 3);
215 GBR_mul(mu_F
[1], mu_F
[1], F_old
);
216 if (GBR_lt(mu_F
[0], mu_F
[1])) {
217 B
= isl_mat_swap_rows(B
, 1 + i
, 1 + i
+ 1);
218 if (i
> tab
->n_zero
) {
220 GBR_set(F_saved
, F_new
);
225 GBR_set(F
[tab
->n_zero
], F_new
);
226 if (ctx
->gbr_only_first
&& GBR_lt(F
[tab
->n_zero
], two
))
230 if (!GBR_is_zero(F
[tab
->n_zero
])) {
231 empty
= GBR_lp_cut(lp
, B
->row
[1+tab
->n_zero
]+1);
234 GBR_set_ui(F
[tab
->n_zero
], 0);
240 GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
243 } while (i
< n_bounded
- 1);
257 for (i
= 0; i
< n_bounded
; ++i
) {
259 GBR_clear(alpha_buffer
[0][i
]);
260 GBR_clear(alpha_buffer
[1][i
]);
263 free(alpha_buffer
[0]);
264 free(alpha_buffer
[1]);
278 isl_int_clear(mu
[0]);
279 isl_int_clear(mu
[1]);
286 struct isl_mat
*isl_basic_set_reduced_basis(struct isl_basic_set
*bset
)
288 struct isl_mat
*basis
;
291 isl_assert(bset
->ctx
, bset
->n_eq
== 0, return NULL
);
293 tab
= isl_tab_from_basic_set(bset
);
294 tab
->basis
= isl_mat_identity(bset
->ctx
, 1 + tab
->n_var
);
295 tab
= isl_tab_compute_reduced_basis(tab
);
299 basis
= isl_mat_copy(tab
->basis
);