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 options->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 * - we have moved forward all the way to the last direction
22 * and then back again all the way to the first direction.
24 struct isl_mat
*isl_basic_set_reduced_basis(struct isl_basic_set
*bset
)
27 struct isl_mat
*basis
;
31 GBR_type F_old
, alpha
, F_new
;
34 struct isl_vec
*b_tmp
;
36 GBR_type
*alpha_buffer
[2] = { NULL
, NULL
};
37 GBR_type
*alpha_saved
;
47 dim
= isl_basic_set_total_dim(bset
);
48 basis
= isl_mat_identity(bset
->ctx
, dim
);
67 b_tmp
= isl_vec_alloc(bset
->ctx
, dim
);
71 F
= isl_alloc_array(bset
->ctx
, GBR_type
, dim
);
72 alpha_buffer
[0] = isl_alloc_array(bset
->ctx
, GBR_type
, dim
);
73 alpha_buffer
[1] = isl_alloc_array(bset
->ctx
, GBR_type
, dim
);
74 alpha_saved
= alpha_buffer
[0];
76 if (!F
|| !alpha_buffer
[0] || !alpha_buffer
[1])
79 for (i
= 0; i
< dim
; ++i
) {
81 GBR_init(alpha_buffer
[0][i
]);
82 GBR_init(alpha_buffer
[1][i
]);
87 lp
= GBR_lp_init(bset
);
93 GBR_lp_set_obj(lp
, basis
->row
[0], dim
);
94 bset
->ctx
->stats
->gbr_solved_lps
++;
95 unbounded
= GBR_lp_solve(lp
);
96 isl_assert(bset
->ctx
, !unbounded
, goto error
);
97 GBR_lp_get_obj_val(lp
, &F
[0]);
101 row
= GBR_lp_next_row(lp
);
102 GBR_set(F_new
, F_saved
);
103 GBR_set(alpha
, alpha_saved
[i
]);
105 row
= GBR_lp_add_row(lp
, basis
->row
[i
], dim
);
106 GBR_lp_set_obj(lp
, basis
->row
[i
+1], dim
);
107 bset
->ctx
->stats
->gbr_solved_lps
++;
108 unbounded
= GBR_lp_solve(lp
);
109 isl_assert(bset
->ctx
, !unbounded
, goto error
);
110 GBR_lp_get_obj_val(lp
, &F_new
);
112 GBR_lp_get_alpha(lp
, row
, &alpha
);
115 save_alpha(lp
, row
-i
, i
, alpha_saved
);
119 GBR_set(F
[i
+1], F_new
);
121 GBR_floor(mu
[0], alpha
);
122 GBR_ceil(mu
[1], alpha
);
124 if (isl_int_eq(mu
[0], mu
[1]))
125 isl_int_set(tmp
, mu
[0]);
129 for (j
= 0; j
<= 1; ++j
) {
130 isl_int_set(tmp
, mu
[j
]);
131 isl_seq_combine(b_tmp
->el
,
132 bset
->ctx
->one
, basis
->row
[i
+1],
133 tmp
, basis
->row
[i
], dim
);
134 GBR_lp_set_obj(lp
, b_tmp
->el
, dim
);
135 bset
->ctx
->stats
->gbr_solved_lps
++;
136 unbounded
= GBR_lp_solve(lp
);
137 isl_assert(bset
->ctx
, !unbounded
, goto error
);
138 GBR_lp_get_obj_val(lp
, &mu_F
[j
]);
140 save_alpha(lp
, row
-i
, i
, alpha_buffer
[j
]);
143 if (GBR_lt(mu_F
[0], mu_F
[1]))
148 isl_int_set(tmp
, mu
[j
]);
149 GBR_set(F_new
, mu_F
[j
]);
150 alpha_saved
= alpha_buffer
[j
];
152 isl_seq_combine(basis
->row
[i
+1],
153 bset
->ctx
->one
, basis
->row
[i
+1],
154 tmp
, basis
->row
[i
], dim
);
156 GBR_set(F_old
, F
[i
]);
159 /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */
160 GBR_set_ui(mu_F
[0], 4);
161 GBR_mul(mu_F
[0], mu_F
[0], F_new
);
162 GBR_set_ui(mu_F
[1], 3);
163 GBR_mul(mu_F
[1], mu_F
[1], F_old
);
164 if (GBR_lt(mu_F
[0], mu_F
[1])) {
165 basis
= isl_mat_swap_rows(basis
, i
, i
+ 1);
168 GBR_set(F_saved
, F_new
);
172 GBR_set(F
[0], F_new
);
173 if (bset
->ctx
->gbr_only_first
&&
178 GBR_lp_add_row(lp
, basis
->row
[i
], dim
);
192 for (i
= 0; i
< dim
; ++i
) {
194 GBR_clear(alpha_buffer
[0][i
]);
195 GBR_clear(alpha_buffer
[1][i
]);
198 free(alpha_buffer
[0]);
199 free(alpha_buffer
[1]);
212 isl_int_clear(mu
[0]);
213 isl_int_clear(mu
[1]);