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->opt->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
;
68 gbr_only_first
= ctx
->opt
->gbr_only_first
;
74 n_bounded
= dim
- tab
->n_unbounded
;
75 if (n_bounded
<= tab
->n_zero
+ 1)
91 b_tmp
= isl_vec_alloc(ctx
, dim
);
95 F
= isl_alloc_array(ctx
, GBR_type
, n_bounded
);
96 alpha_buffer
[0] = isl_alloc_array(ctx
, GBR_type
, n_bounded
);
97 alpha_buffer
[1] = isl_alloc_array(ctx
, GBR_type
, n_bounded
);
98 alpha_saved
= alpha_buffer
[0];
100 if (!F
|| !alpha_buffer
[0] || !alpha_buffer
[1])
103 for (i
= 0; i
< n_bounded
; ++i
) {
105 GBR_init(alpha_buffer
[0][i
]);
106 GBR_init(alpha_buffer
[1][i
]);
112 lp
= GBR_lp_init(tab
);
118 GBR_lp_set_obj(lp
, B
->row
[1+i
]+1, dim
);
119 ctx
->stats
->gbr_solved_lps
++;
120 unbounded
= GBR_lp_solve(lp
);
121 isl_assert(ctx
, !unbounded
, goto error
);
122 GBR_lp_get_obj_val(lp
, &F
[i
]);
124 if (GBR_lt(F
[i
], one
)) {
125 if (!GBR_is_zero(F
[i
])) {
126 empty
= GBR_lp_cut(lp
, B
->row
[1+i
]+1);
135 if (i
+1 == tab
->n_zero
) {
136 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
137 ctx
->stats
->gbr_solved_lps
++;
138 unbounded
= GBR_lp_solve(lp
);
139 isl_assert(ctx
, !unbounded
, goto error
);
140 GBR_lp_get_obj_val(lp
, &F_new
);
141 fixed
= GBR_lp_is_fixed(lp
);
142 GBR_set_ui(alpha
, 0);
145 row
= GBR_lp_next_row(lp
);
146 GBR_set(F_new
, F_saved
);
148 GBR_set(alpha
, alpha_saved
[i
]);
150 row
= GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
151 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
152 ctx
->stats
->gbr_solved_lps
++;
153 unbounded
= GBR_lp_solve(lp
);
154 isl_assert(ctx
, !unbounded
, goto error
);
155 GBR_lp_get_obj_val(lp
, &F_new
);
156 fixed
= GBR_lp_is_fixed(lp
);
158 GBR_lp_get_alpha(lp
, row
, &alpha
);
161 save_alpha(lp
, row
-i
, i
, alpha_saved
);
163 if (GBR_lp_del_row(lp
) < 0)
166 GBR_set(F
[i
+1], F_new
);
168 GBR_floor(mu
[0], alpha
);
169 GBR_ceil(mu
[1], alpha
);
171 if (isl_int_eq(mu
[0], mu
[1]))
172 isl_int_set(tmp
, mu
[0]);
176 for (j
= 0; j
<= 1; ++j
) {
177 isl_int_set(tmp
, mu
[j
]);
178 isl_seq_combine(b_tmp
->el
,
179 ctx
->one
, B
->row
[1+i
+1]+1,
180 tmp
, B
->row
[1+i
]+1, dim
);
181 GBR_lp_set_obj(lp
, b_tmp
->el
, dim
);
182 ctx
->stats
->gbr_solved_lps
++;
183 unbounded
= GBR_lp_solve(lp
);
184 isl_assert(ctx
, !unbounded
, goto error
);
185 GBR_lp_get_obj_val(lp
, &mu_F
[j
]);
186 mu_fixed
[j
] = GBR_lp_is_fixed(lp
);
188 save_alpha(lp
, row
-i
, i
, alpha_buffer
[j
]);
191 if (GBR_lt(mu_F
[0], mu_F
[1]))
196 isl_int_set(tmp
, mu
[j
]);
197 GBR_set(F_new
, mu_F
[j
]);
199 alpha_saved
= alpha_buffer
[j
];
201 isl_seq_combine(B
->row
[1+i
+1]+1, ctx
->one
, B
->row
[1+i
+1]+1,
202 tmp
, B
->row
[1+i
]+1, dim
);
204 if (i
+1 == tab
->n_zero
&& fixed
) {
205 if (!GBR_is_zero(F
[i
+1])) {
206 empty
= GBR_lp_cut(lp
, B
->row
[1+i
+1]+1);
209 GBR_set_ui(F
[i
+1], 0);
214 GBR_set(F_old
, F
[i
]);
217 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
218 GBR_set_ui(mu_F
[0], 4);
219 GBR_mul(mu_F
[0], mu_F
[0], F_new
);
220 GBR_set_ui(mu_F
[1], 3);
221 GBR_mul(mu_F
[1], mu_F
[1], F_old
);
222 if (GBR_lt(mu_F
[0], mu_F
[1])) {
223 B
= isl_mat_swap_rows(B
, 1 + i
, 1 + i
+ 1);
224 if (i
> tab
->n_zero
) {
226 GBR_set(F_saved
, F_new
);
228 if (GBR_lp_del_row(lp
) < 0)
232 GBR_set(F
[tab
->n_zero
], F_new
);
233 if (gbr_only_first
&& GBR_lt(F
[tab
->n_zero
], two
))
237 if (!GBR_is_zero(F
[tab
->n_zero
])) {
238 empty
= GBR_lp_cut(lp
, B
->row
[1+tab
->n_zero
]+1);
241 GBR_set_ui(F
[tab
->n_zero
], 0);
247 GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
250 } while (i
< n_bounded
- 1);
264 for (i
= 0; i
< n_bounded
; ++i
) {
266 GBR_clear(alpha_buffer
[0][i
]);
267 GBR_clear(alpha_buffer
[1][i
]);
270 free(alpha_buffer
[0]);
271 free(alpha_buffer
[1]);
285 isl_int_clear(mu
[0]);
286 isl_int_clear(mu
[1]);
293 struct isl_mat
*isl_basic_set_reduced_basis(struct isl_basic_set
*bset
)
295 struct isl_mat
*basis
;
298 isl_assert(bset
->ctx
, bset
->n_eq
== 0, return NULL
);
300 tab
= isl_tab_from_basic_set(bset
);
301 tab
->basis
= isl_mat_identity(bset
->ctx
, 1 + tab
->n_var
);
302 tab
= isl_tab_compute_reduced_basis(tab
);
306 basis
= isl_mat_copy(tab
->basis
);