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
17struct 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
90static struct tab_lp *init_lp(struct isl_tab *tab);
91static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
92static int solve_lp(struct tab_lp *lp);
93static void get_obj_val(struct tab_lp* lp, GBR_type *F);
94static void delete_lp(struct tab_lp *lp);
95static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
96static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha);
97static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;
98static 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 */
118static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row)
119{
120 unsigned dim;
121 struct isl_tab *prod;
122
123 if (!tab || !row)
124 return NULL;
125
126 dim = tab->n_var;
127 prod = isl_tab_product(tab1: tab, tab2: tab);
128 if (isl_tab_extend_cons(tab: prod, n_new: 3 * dim + 1) < 0) {
129 isl_tab_free(tab: prod);
130 return NULL;
131 }
132 return prod;
133}
134
135static struct tab_lp *init_lp(struct isl_tab *tab)
136{
137 struct tab_lp *lp = NULL;
138
139 if (!tab)
140 return NULL;
141
142 lp = isl_calloc_type(tab->mat->ctx, struct tab_lp);
143 if (!lp)
144 return NULL;
145
146 isl_int_init(lp->opt);
147 isl_int_init(lp->opt_denom);
148 isl_int_init(lp->tmp);
149 isl_int_init(lp->tmp2);
150
151 lp->dim = tab->n_var;
152
153 lp->ctx = tab->mat->ctx;
154 isl_ctx_ref(ctx: lp->ctx);
155
156 lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);
157
158 lp->row = isl_vec_alloc(ctx: lp->ctx, size: 1 + 2 * lp->dim);
159 if (!lp->row)
160 goto error;
161 lp->tab = gbr_tab(tab, row: lp->row);
162 if (!lp->tab)
163 goto error;
164 lp->con_offset = lp->tab->n_con;
165 lp->obj = NULL;
166 lp->neq = 0;
167
168 return lp;
169error:
170 delete_lp(lp);
171 return NULL;
172}
173
174static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
175{
176 lp->obj = row;
177}
178
179static int solve_lp(struct tab_lp *lp)
180{
181 enum isl_lp_result res;
182 unsigned flags = 0;
183
184 lp->is_fixed = 0;
185
186 isl_int_set_si(lp->row->el[0], 0);
187 isl_seq_cpy(dst: lp->row->el + 1, src: lp->obj, len: lp->dim);
188 isl_seq_neg(dst: lp->row->el + 1 + lp->dim, src: lp->obj, len: lp->dim);
189 if (lp->neq)
190 flags = ISL_TAB_SAVE_DUAL;
191 res = isl_tab_min(tab: lp->tab, f: lp->row->el, denom: lp->ctx->one,
192 opt: &lp->opt, opt_denom: &lp->opt_denom, flags);
193 isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2);
194 if (isl_int_abs_lt(lp->opt, lp->opt_denom)) {
195 struct isl_vec *sample = isl_tab_get_sample_value(tab: lp->tab);
196 if (!sample)
197 return -1;
198 isl_seq_inner_product(p1: lp->obj, p2: sample->el + 1, len: lp->dim, prod: &lp->tmp);
199 isl_seq_inner_product(p1: lp->obj, p2: sample->el + 1 + lp->dim, len: lp->dim, prod: &lp->tmp2);
200 isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]);
201 isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]);
202 if (isl_int_ge(lp->tmp, lp->tmp2))
203 lp->is_fixed = 1;
204 isl_vec_free(vec: sample);
205 }
206 isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2);
207 if (res < 0)
208 return -1;
209 if (res != isl_lp_ok)
210 isl_die(lp->ctx, isl_error_internal,
211 "unexpected missing (bounded) solution", return -1);
212 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 */
220static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row)
221{
222 enum isl_lp_result res;
223
224 isl_int_set_si(lp->row->el[0], 0);
225 isl_seq_cpy(dst: lp->row->el + 1, src: row, len: lp->dim);
226 isl_seq_clr(p: lp->row->el + 1 + lp->dim, len: lp->dim);
227 res = isl_tab_min(tab: lp->tab, f: lp->row->el, denom: lp->ctx->one,
228 opt: &lp->tmp, NULL, flags: 0);
229 if (res != isl_lp_ok)
230 return -1;
231
232 isl_int_neg(lp->row->el[0], lp->tmp);
233 if (isl_tab_add_eq(tab: lp->tab, eq: lp->row->el) < 0)
234 return -1;
235
236 isl_seq_cpy(dst: lp->row->el + 1 + lp->dim, src: row, len: lp->dim);
237 isl_seq_clr(p: lp->row->el + 1, len: lp->dim);
238 if (isl_tab_add_eq(tab: lp->tab, eq: lp->row->el) < 0)
239 return -1;
240
241 lp->con_offset += 2;
242
243 return lp->tab->empty;
244}
245
246static void get_obj_val(struct tab_lp* lp, GBR_type *F)
247{
248 GBR_set_num_neg(F, lp->opt);
249 GBR_set_den(F, lp->opt_denom);
250}
251
252static void delete_lp(struct tab_lp *lp)
253{
254 if (!lp)
255 return;
256
257 isl_int_clear(lp->opt);
258 isl_int_clear(lp->opt_denom);
259 isl_int_clear(lp->tmp);
260 isl_int_clear(lp->tmp2);
261 isl_vec_free(vec: lp->row);
262 free(ptr: lp->stack);
263 isl_tab_free(tab: lp->tab);
264 isl_ctx_deref(ctx: lp->ctx);
265 free(ptr: lp);
266}
267
268static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)
269{
270 lp->stack[lp->neq] = isl_tab_snap(tab: lp->tab);
271
272 isl_int_set_si(lp->row->el[0], 0);
273 isl_seq_cpy(dst: lp->row->el + 1, src: row, len: lp->dim);
274 isl_seq_neg(dst: lp->row->el + 1 + lp->dim, src: row, len: lp->dim);
275
276 if (isl_tab_add_valid_eq(tab: lp->tab, eq: lp->row->el) < 0)
277 return -1;
278
279 return lp->neq++;
280}
281
282static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha)
283{
284 row += lp->con_offset;
285 GBR_set_num_neg(alpha, lp->tab->dual->el[1 + row]);
286 GBR_set_den(alpha, lp->tab->dual->el[0]);
287}
288
289static int del_lp_row(struct tab_lp *lp)
290{
291 lp->neq--;
292 return isl_tab_rollback(tab: lp->tab, snap: lp->stack[lp->neq]);
293}
294

source code of polly/lib/External/isl/basis_reduction_tab.c