Line data Source code
1 : /*
2 : * Copyright 2008-2009 Katholieke Universiteit Leuven
3 : *
4 : * Use of this software is governed by the MIT license
5 : *
6 : * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 : * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 : */
9 :
10 : #include <assert.h>
11 : #include <isl_map_private.h>
12 : #include <isl_seq.h>
13 : #include "isl_tab.h"
14 : #include <isl_int.h>
15 : #include <isl_config.h>
16 :
17 : struct tab_lp {
18 : struct isl_ctx *ctx;
19 : struct isl_vec *row;
20 : struct isl_tab *tab;
21 : struct isl_tab_undo **stack;
22 : isl_int *obj;
23 : isl_int opt;
24 : isl_int opt_denom;
25 : isl_int tmp;
26 : isl_int tmp2;
27 : int neq;
28 : unsigned dim;
29 : /* number of constraints in initial product tableau */
30 : int con_offset;
31 : /* objective function has fixed or no integer value */
32 : int is_fixed;
33 : };
34 :
35 : #ifdef USE_GMP_FOR_MP
36 : #define GBR_type mpq_t
37 : #define GBR_init(v) mpq_init(v)
38 : #define GBR_clear(v) mpq_clear(v)
39 : #define GBR_set(a,b) mpq_set(a,b)
40 : #define GBR_set_ui(a,b) mpq_set_ui(a,b,1)
41 : #define GBR_mul(a,b,c) mpq_mul(a,b,c)
42 : #define GBR_lt(a,b) (mpq_cmp(a,b) < 0)
43 : #define GBR_is_zero(a) (mpq_sgn(a) == 0)
44 : #define GBR_numref(a) mpq_numref(a)
45 : #define GBR_denref(a) mpq_denref(a)
46 : #define GBR_floor(a,b) mpz_fdiv_q(a,GBR_numref(b),GBR_denref(b))
47 : #define GBR_ceil(a,b) mpz_cdiv_q(a,GBR_numref(b),GBR_denref(b))
48 : #define GBR_set_num_neg(a, b) mpz_neg(GBR_numref(*a), b);
49 : #define GBR_set_den(a, b) mpz_set(GBR_denref(*a), b);
50 : #endif /* USE_GMP_FOR_MP */
51 :
52 : #ifdef USE_IMATH_FOR_MP
53 : #include <imrat.h>
54 :
55 : #define GBR_type mp_rat
56 : #define GBR_init(v) v = mp_rat_alloc()
57 : #define GBR_clear(v) mp_rat_free(v)
58 : #define GBR_set(a,b) mp_rat_copy(b,a)
59 : #define GBR_set_ui(a,b) mp_rat_set_uvalue(a,b,1)
60 : #define GBR_mul(a,b,c) mp_rat_mul(b,c,a)
61 : #define GBR_lt(a,b) (mp_rat_compare(a,b) < 0)
62 : #define GBR_is_zero(a) (mp_rat_compare_zero(a) == 0)
63 : #ifdef USE_SMALL_INT_OPT
64 : #define GBR_numref(a) isl_sioimath_encode_big(mp_rat_numer_ref(a))
65 : #define GBR_denref(a) isl_sioimath_encode_big(mp_rat_denom_ref(a))
66 : #define GBR_floor(a, b) isl_sioimath_fdiv_q((a), GBR_numref(b), GBR_denref(b))
67 : #define GBR_ceil(a, b) isl_sioimath_cdiv_q((a), GBR_numref(b), GBR_denref(b))
68 : #define GBR_set_num_neg(a, b) \
69 : do { \
70 : isl_sioimath_scratchspace_t scratch; \
71 : impz_neg(mp_rat_numer_ref(*a), \
72 : isl_sioimath_bigarg_src(*b, &scratch));\
73 : } while (0)
74 : #define GBR_set_den(a, b) \
75 : do { \
76 : isl_sioimath_scratchspace_t scratch; \
77 : impz_set(mp_rat_denom_ref(*a), \
78 : isl_sioimath_bigarg_src(*b, &scratch));\
79 : } while (0)
80 : #else /* USE_SMALL_INT_OPT */
81 : #define GBR_numref(a) mp_rat_numer_ref(a)
82 : #define GBR_denref(a) mp_rat_denom_ref(a)
83 : #define GBR_floor(a,b) impz_fdiv_q(a,GBR_numref(b),GBR_denref(b))
84 : #define GBR_ceil(a,b) impz_cdiv_q(a,GBR_numref(b),GBR_denref(b))
85 : #define GBR_set_num_neg(a, b) impz_neg(GBR_numref(*a), b)
86 : #define GBR_set_den(a, b) impz_set(GBR_denref(*a), b)
87 : #endif /* USE_SMALL_INT_OPT */
88 : #endif /* USE_IMATH_FOR_MP */
89 :
90 : static struct tab_lp *init_lp(struct isl_tab *tab);
91 : static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
92 : static int solve_lp(struct tab_lp *lp);
93 : static void get_obj_val(struct tab_lp* lp, GBR_type *F);
94 : static void delete_lp(struct tab_lp *lp);
95 : static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
96 : static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha);
97 : static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;
98 : static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row);
99 :
100 : #define GBR_LP struct tab_lp
101 : #define GBR_lp_init(P) init_lp(P)
102 : #define GBR_lp_set_obj(lp, obj, dim) set_lp_obj(lp, obj, dim)
103 : #define GBR_lp_solve(lp) solve_lp(lp)
104 : #define GBR_lp_get_obj_val(lp, F) get_obj_val(lp, F)
105 : #define GBR_lp_delete(lp) delete_lp(lp)
106 : #define GBR_lp_next_row(lp) lp->neq
107 : #define GBR_lp_add_row(lp, row, dim) add_lp_row(lp, row, dim)
108 : #define GBR_lp_get_alpha(lp, row, alpha) get_alpha(lp, row, alpha)
109 : #define GBR_lp_del_row(lp) del_lp_row(lp)
110 : #define GBR_lp_is_fixed(lp) (lp)->is_fixed
111 : #define GBR_lp_cut(lp, obj) cut_lp_to_hyperplane(lp, obj)
112 : #include "basis_reduction_templ.c"
113 :
114 : /* Set up a tableau for the Cartesian product of bset with itself.
115 : * This could be optimized by first setting up a tableau for bset
116 : * and then performing the Cartesian product on the tableau.
117 : */
118 3355293 : static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row)
119 : {
120 : unsigned dim;
121 : struct isl_tab *prod;
122 :
123 3355293 : if (!tab || !row)
124 0 : return NULL;
125 :
126 3355293 : dim = tab->n_var;
127 3355293 : prod = isl_tab_product(tab, tab);
128 3355293 : if (isl_tab_extend_cons(prod, 3 * dim + 1) < 0) {
129 0 : isl_tab_free(prod);
130 0 : return NULL;
131 : }
132 3355293 : return prod;
133 : }
134 :
135 3355293 : static struct tab_lp *init_lp(struct isl_tab *tab)
136 : {
137 3355293 : struct tab_lp *lp = NULL;
138 :
139 3355293 : if (!tab)
140 0 : return NULL;
141 :
142 3355293 : lp = isl_calloc_type(tab->mat->ctx, struct tab_lp);
143 3355293 : if (!lp)
144 0 : return NULL;
145 :
146 3355293 : isl_int_init(lp->opt);
147 3355293 : isl_int_init(lp->opt_denom);
148 3355293 : isl_int_init(lp->tmp);
149 3355293 : isl_int_init(lp->tmp2);
150 :
151 3355293 : lp->dim = tab->n_var;
152 :
153 3355293 : lp->ctx = tab->mat->ctx;
154 3355293 : isl_ctx_ref(lp->ctx);
155 :
156 3355293 : lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);
157 :
158 3355293 : lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);
159 3355293 : if (!lp->row)
160 0 : goto error;
161 3355293 : lp->tab = gbr_tab(tab, lp->row);
162 3355293 : if (!lp->tab)
163 0 : goto error;
164 3355293 : lp->con_offset = lp->tab->n_con;
165 3355293 : lp->obj = NULL;
166 3355293 : lp->neq = 0;
167 :
168 3355293 : return lp;
169 : error:
170 0 : delete_lp(lp);
171 0 : return NULL;
172 : }
173 :
174 71963008 : static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
175 : {
176 71963008 : lp->obj = row;
177 71963008 : }
178 :
179 71963008 : static int solve_lp(struct tab_lp *lp)
180 : {
181 : enum isl_lp_result res;
182 71963008 : unsigned flags = 0;
183 :
184 71963008 : lp->is_fixed = 0;
185 :
186 71963008 : isl_int_set_si(lp->row->el[0], 0);
187 71963008 : isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);
188 71963008 : isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);
189 71963008 : if (lp->neq)
190 55872265 : flags = ISL_TAB_SAVE_DUAL;
191 71963008 : res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
192 : &lp->opt, &lp->opt_denom, flags);
193 71963008 : isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2);
194 71963008 : if (isl_int_abs_lt(lp->opt, lp->opt_denom)) {
195 8863211 : struct isl_vec *sample = isl_tab_get_sample_value(lp->tab);
196 8863211 : if (!sample)
197 0 : return -1;
198 8863211 : isl_seq_inner_product(lp->obj, sample->el + 1, lp->dim, &lp->tmp);
199 8863211 : isl_seq_inner_product(lp->obj, sample->el + 1 + lp->dim, lp->dim, &lp->tmp2);
200 8863211 : isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]);
201 8863211 : isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]);
202 8863211 : if (isl_int_ge(lp->tmp, lp->tmp2))
203 6234584 : lp->is_fixed = 1;
204 8863211 : isl_vec_free(sample);
205 : }
206 71963008 : isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2);
207 71963008 : if (res < 0)
208 0 : return -1;
209 71963008 : if (res != isl_lp_ok)
210 0 : isl_die(lp->ctx, isl_error_internal,
211 : "unexpected missing (bounded) solution", return -1);
212 71963008 : return 0;
213 : }
214 :
215 : /* The current objective function has a fixed (or no) integer value.
216 : * Cut the tableau to the hyperplane that fixes this value in
217 : * both halves of the tableau.
218 : * Return 1 if the resulting tableau is empty.
219 : */
220 0 : static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row)
221 : {
222 : enum isl_lp_result res;
223 :
224 0 : isl_int_set_si(lp->row->el[0], 0);
225 0 : isl_seq_cpy(lp->row->el + 1, row, lp->dim);
226 0 : isl_seq_clr(lp->row->el + 1 + lp->dim, lp->dim);
227 0 : res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
228 : &lp->tmp, NULL, 0);
229 0 : if (res != isl_lp_ok)
230 0 : return -1;
231 :
232 0 : isl_int_neg(lp->row->el[0], lp->tmp);
233 0 : if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)
234 0 : return -1;
235 :
236 0 : isl_seq_cpy(lp->row->el + 1 + lp->dim, row, lp->dim);
237 0 : isl_seq_clr(lp->row->el + 1, lp->dim);
238 0 : if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)
239 0 : return -1;
240 :
241 0 : lp->con_offset += 2;
242 :
243 0 : return lp->tab->empty;
244 : }
245 :
246 71963008 : static void get_obj_val(struct tab_lp* lp, GBR_type *F)
247 : {
248 71963008 : GBR_set_num_neg(F, lp->opt);
249 71963008 : GBR_set_den(F, lp->opt_denom);
250 71963008 : }
251 :
252 3355293 : static void delete_lp(struct tab_lp *lp)
253 : {
254 3355293 : if (!lp)
255 0 : return;
256 :
257 3355293 : isl_int_clear(lp->opt);
258 3355293 : isl_int_clear(lp->opt_denom);
259 3355293 : isl_int_clear(lp->tmp);
260 3355293 : isl_int_clear(lp->tmp2);
261 3355293 : isl_vec_free(lp->row);
262 3355293 : free(lp->stack);
263 3355293 : isl_tab_free(lp->tab);
264 3355293 : isl_ctx_deref(lp->ctx);
265 3355293 : free(lp);
266 : }
267 :
268 38971061 : static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)
269 : {
270 38971061 : lp->stack[lp->neq] = isl_tab_snap(lp->tab);
271 :
272 38971061 : isl_int_set_si(lp->row->el[0], 0);
273 38971061 : isl_seq_cpy(lp->row->el + 1, row, lp->dim);
274 38971061 : isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);
275 :
276 38971061 : if (isl_tab_add_valid_eq(lp->tab, lp->row->el) < 0)
277 0 : return -1;
278 :
279 38971061 : return lp->neq++;
280 : }
281 :
282 133503935 : static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha)
283 : {
284 133503935 : row += lp->con_offset;
285 133503935 : GBR_set_num_neg(alpha, lp->tab->dual->el[1 + row]);
286 133503935 : GBR_set_den(alpha, lp->tab->dual->el[0]);
287 133503935 : }
288 :
289 32810237 : static int del_lp_row(struct tab_lp *lp)
290 : {
291 32810237 : lp->neq--;
292 32810237 : return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);
293 : }
|
