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.
20 * This function implements the algorithm described in
21 * "An Implementation of the Generalized Basis Reduction Algorithm
22 * for Integer Programming" of Cook el al. to compute a reduced basis.
23 * We use \epsilon = 1/4.
25 * If ctx->gbr_only_first is set, the user is only interested
26 * in the first direction. In this case we stop the basis reduction when
27 * the width in the first direction becomes smaller than 2.
29 struct isl_tab
*isl_tab_compute_reduced_basis(struct isl_tab
*tab
)
37 GBR_type F_old
, alpha
, F_new
;
40 struct isl_vec
*b_tmp
;
42 GBR_type
*alpha_buffer
[2] = { NULL
, NULL
};
43 GBR_type
*alpha_saved
;
64 if (dim
<= tab
->n_zero
+ 1)
80 b_tmp
= isl_vec_alloc(ctx
, dim
);
84 F
= isl_alloc_array(ctx
, GBR_type
, dim
);
85 alpha_buffer
[0] = isl_alloc_array(ctx
, GBR_type
, dim
);
86 alpha_buffer
[1] = isl_alloc_array(ctx
, GBR_type
, dim
);
87 alpha_saved
= alpha_buffer
[0];
89 if (!F
|| !alpha_buffer
[0] || !alpha_buffer
[1])
92 for (i
= 0; i
< dim
; ++i
) {
94 GBR_init(alpha_buffer
[0][i
]);
95 GBR_init(alpha_buffer
[1][i
]);
101 lp
= GBR_lp_init(tab
);
107 GBR_lp_set_obj(lp
, B
->row
[1+i
]+1, dim
);
108 ctx
->stats
->gbr_solved_lps
++;
109 unbounded
= GBR_lp_solve(lp
);
110 isl_assert(ctx
, !unbounded
, goto error
);
111 GBR_lp_get_obj_val(lp
, &F
[i
]);
113 if (GBR_lt(F
[i
], one
)) {
114 if (!GBR_is_zero(F
[i
])) {
115 empty
= GBR_lp_cut(lp
, B
->row
[1+i
]+1);
124 if (i
+1 == tab
->n_zero
) {
125 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
126 ctx
->stats
->gbr_solved_lps
++;
127 unbounded
= GBR_lp_solve(lp
);
128 isl_assert(ctx
, !unbounded
, goto error
);
129 GBR_lp_get_obj_val(lp
, &F_new
);
130 fixed
= GBR_lp_is_fixed(lp
);
131 GBR_set_ui(alpha
, 0);
134 row
= GBR_lp_next_row(lp
);
135 GBR_set(F_new
, F_saved
);
137 GBR_set(alpha
, alpha_saved
[i
]);
139 row
= GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
140 GBR_lp_set_obj(lp
, B
->row
[1+i
+1]+1, dim
);
141 ctx
->stats
->gbr_solved_lps
++;
142 unbounded
= GBR_lp_solve(lp
);
143 isl_assert(ctx
, !unbounded
, goto error
);
144 GBR_lp_get_obj_val(lp
, &F_new
);
145 fixed
= GBR_lp_is_fixed(lp
);
147 GBR_lp_get_alpha(lp
, row
, &alpha
);
150 save_alpha(lp
, row
-i
, i
, alpha_saved
);
154 GBR_set(F
[i
+1], F_new
);
156 GBR_floor(mu
[0], alpha
);
157 GBR_ceil(mu
[1], alpha
);
159 if (isl_int_eq(mu
[0], mu
[1]))
160 isl_int_set(tmp
, mu
[0]);
164 for (j
= 0; j
<= 1; ++j
) {
165 isl_int_set(tmp
, mu
[j
]);
166 isl_seq_combine(b_tmp
->el
,
167 ctx
->one
, B
->row
[1+i
+1]+1,
168 tmp
, B
->row
[1+i
]+1, dim
);
169 GBR_lp_set_obj(lp
, b_tmp
->el
, dim
);
170 ctx
->stats
->gbr_solved_lps
++;
171 unbounded
= GBR_lp_solve(lp
);
172 isl_assert(ctx
, !unbounded
, goto error
);
173 GBR_lp_get_obj_val(lp
, &mu_F
[j
]);
174 mu_fixed
[j
] = GBR_lp_is_fixed(lp
);
176 save_alpha(lp
, row
-i
, i
, alpha_buffer
[j
]);
179 if (GBR_lt(mu_F
[0], mu_F
[1]))
184 isl_int_set(tmp
, mu
[j
]);
185 GBR_set(F_new
, mu_F
[j
]);
187 alpha_saved
= alpha_buffer
[j
];
189 isl_seq_combine(B
->row
[1+i
+1]+1, ctx
->one
, B
->row
[1+i
+1]+1,
190 tmp
, B
->row
[1+i
]+1, dim
);
192 if (i
+1 == tab
->n_zero
&& fixed
) {
193 if (!GBR_is_zero(F
[i
+1])) {
194 empty
= GBR_lp_cut(lp
, B
->row
[1+i
+1]+1);
197 GBR_set_ui(F
[i
+1], 0);
202 GBR_set(F_old
, F
[i
]);
205 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
206 GBR_set_ui(mu_F
[0], 4);
207 GBR_mul(mu_F
[0], mu_F
[0], F_new
);
208 GBR_set_ui(mu_F
[1], 3);
209 GBR_mul(mu_F
[1], mu_F
[1], F_old
);
210 if (GBR_lt(mu_F
[0], mu_F
[1])) {
211 B
= isl_mat_swap_rows(B
, 1 + i
, 1 + i
+ 1);
212 if (i
> tab
->n_zero
) {
214 GBR_set(F_saved
, F_new
);
219 GBR_set(F
[tab
->n_zero
], F_new
);
220 if (ctx
->gbr_only_first
&& GBR_lt(F
[tab
->n_zero
], two
))
224 if (!GBR_is_zero(F
[tab
->n_zero
])) {
225 empty
= GBR_lp_cut(lp
, B
->row
[1+tab
->n_zero
]+1);
228 GBR_set_ui(F
[tab
->n_zero
], 0);
234 GBR_lp_add_row(lp
, B
->row
[1+i
]+1, dim
);
251 for (i
= 0; i
< dim
; ++i
) {
253 GBR_clear(alpha_buffer
[0][i
]);
254 GBR_clear(alpha_buffer
[1][i
]);
257 free(alpha_buffer
[0]);
258 free(alpha_buffer
[1]);
272 isl_int_clear(mu
[0]);
273 isl_int_clear(mu
[1]);
280 struct isl_mat
*isl_basic_set_reduced_basis(struct isl_basic_set
*bset
)
282 struct isl_mat
*basis
;
285 isl_assert(bset
->ctx
, bset
->n_eq
== 0, return NULL
);
287 tab
= isl_tab_from_basic_set(bset
);
288 tab
->basis
= isl_mat_identity(bset
->ctx
, 1 + tab
->n_var
);
289 tab
= isl_tab_compute_reduced_basis(tab
);
293 basis
= isl_mat_copy(tab
->basis
);