Merge branch 'maint'
[platform/upstream/isl.git] / basis_reduction_tab.c
1 /*
2  * Copyright 2008-2009 Katholieke Universiteit Leuven
3  *
4  * Use of this software is governed by the GNU LGPLv2.1 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/seq.h>
12 #include "isl_tab.h"
13
14 struct tab_lp {
15         struct isl_ctx  *ctx;
16         struct isl_vec  *row;
17         struct isl_tab  *tab;
18         struct isl_tab_undo     **stack;
19         isl_int         *obj;
20         isl_int          opt;
21         isl_int          opt_denom;
22         isl_int          tmp;
23         isl_int          tmp2;
24         int              neq;
25         unsigned         dim;
26         /* number of constraints in initial product tableau */
27         int              con_offset;
28         /* objective function has fixed or no integer value */
29         int              is_fixed;
30 };
31
32 static struct tab_lp *init_lp(struct isl_tab *tab);
33 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
34 static int solve_lp(struct tab_lp *lp);
35 static void get_obj_val(struct tab_lp* lp, mpq_t *F);
36 static void delete_lp(struct tab_lp *lp);
37 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
38 static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha);
39 static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;
40 static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row);
41
42 #define GBR_LP                              struct tab_lp
43 #define GBR_type                            mpq_t
44 #define GBR_init(v)                         mpq_init(v)
45 #define GBR_clear(v)                        mpq_clear(v)
46 #define GBR_set(a,b)                        mpq_set(a,b)
47 #define GBR_set_ui(a,b)                     mpq_set_ui(a,b,1)
48 #define GBR_mul(a,b,c)                      mpq_mul(a,b,c)
49 #define GBR_lt(a,b)                         (mpq_cmp(a,b) < 0)
50 #define GBR_is_zero(a)                      (mpq_sgn(a) == 0)
51 #define GBR_floor(a,b)                      mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b))
52 #define GBR_ceil(a,b)                       mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b))
53 #define GBR_lp_init(P)                      init_lp(P)
54 #define GBR_lp_set_obj(lp, obj, dim)        set_lp_obj(lp, obj, dim)
55 #define GBR_lp_solve(lp)                    solve_lp(lp)
56 #define GBR_lp_get_obj_val(lp, F)           get_obj_val(lp, F)
57 #define GBR_lp_delete(lp)                   delete_lp(lp)
58 #define GBR_lp_next_row(lp)                 lp->neq
59 #define GBR_lp_add_row(lp, row, dim)        add_lp_row(lp, row, dim)
60 #define GBR_lp_get_alpha(lp, row, alpha)    get_alpha(lp, row, alpha)
61 #define GBR_lp_del_row(lp)                  del_lp_row(lp)
62 #define GBR_lp_is_fixed(lp)                 (lp)->is_fixed
63 #define GBR_lp_cut(lp, obj)                 cut_lp_to_hyperplane(lp, obj)
64 #include "basis_reduction_templ.c"
65
66 /* Set up a tableau for the Cartesian product of bset with itself.
67  * This could be optimized by first setting up a tableau for bset
68  * and then performing the Cartesian product on the tableau.
69  */
70 static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row)
71 {
72         unsigned dim;
73         struct isl_tab *prod;
74
75         if (!tab || !row)
76                 return NULL;
77
78         dim = tab->n_var;
79         prod = isl_tab_product(tab, tab);
80         if (isl_tab_extend_cons(prod, 3 * dim + 1) < 0) {
81                 isl_tab_free(prod);
82                 return NULL;
83         }
84         return prod;
85 }
86
87 static struct tab_lp *init_lp(struct isl_tab *tab)
88 {
89         struct tab_lp *lp = NULL;
90
91         if (!tab)
92                 return NULL;
93
94         lp = isl_calloc_type(tab->mat->ctx, struct tab_lp);
95         if (!lp)
96                 return NULL;
97
98         isl_int_init(lp->opt);
99         isl_int_init(lp->opt_denom);
100         isl_int_init(lp->tmp);
101         isl_int_init(lp->tmp2);
102
103         lp->dim = tab->n_var;
104
105         lp->ctx = tab->mat->ctx;
106         isl_ctx_ref(lp->ctx);
107
108         lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);
109
110         lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);
111         if (!lp->row)
112                 goto error;
113         lp->tab = gbr_tab(tab, lp->row);
114         if (!lp->tab)
115                 goto error;
116         lp->con_offset = lp->tab->n_con;
117         lp->obj = NULL;
118         lp->neq = 0;
119
120         return lp;
121 error:
122         delete_lp(lp);
123         return NULL;
124 }
125
126 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
127 {
128         lp->obj = row;
129 }
130
131 static int solve_lp(struct tab_lp *lp)
132 {
133         enum isl_lp_result res;
134         unsigned flags = 0;
135
136         lp->is_fixed = 0;
137
138         isl_int_set_si(lp->row->el[0], 0);
139         isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);
140         isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);
141         if (lp->neq)
142                 flags = ISL_TAB_SAVE_DUAL;
143         res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
144                           &lp->opt, &lp->opt_denom, flags);
145         isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2);
146         if (isl_int_abs_lt(lp->opt, lp->opt_denom)) {
147                 struct isl_vec *sample = isl_tab_get_sample_value(lp->tab);
148                 if (!sample)
149                         return -1;
150                 isl_seq_inner_product(lp->obj, sample->el + 1, lp->dim, &lp->tmp);
151                 isl_seq_inner_product(lp->obj, sample->el + 1 + lp->dim, lp->dim, &lp->tmp2);
152                 isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]);
153                 isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]);
154                 if (isl_int_ge(lp->tmp, lp->tmp2))
155                         lp->is_fixed = 1;
156                 isl_vec_free(sample);
157         }
158         isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2);
159         if (res != isl_lp_ok)
160                 return -1;
161         return 0;
162 }
163
164 /* The current objective function has a fixed (or no) integer value.
165  * Cut the tableau to the hyperplane that fixes this value in
166  * both halves of the tableau.
167  * Return 1 if the resulting tableau is empty.
168  */
169 static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row)
170 {
171         enum isl_lp_result res;
172
173         isl_int_set_si(lp->row->el[0], 0);
174         isl_seq_cpy(lp->row->el + 1, row, lp->dim);
175         isl_seq_clr(lp->row->el + 1 + lp->dim, lp->dim);
176         res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
177                           &lp->tmp, NULL, 0);
178         if (res != isl_lp_ok)
179                 return -1;
180
181         isl_int_neg(lp->row->el[0], lp->tmp);
182         if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)
183                 return -1;
184
185         isl_seq_cpy(lp->row->el + 1 + lp->dim, row, lp->dim);
186         isl_seq_clr(lp->row->el + 1, lp->dim);
187         if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)
188                 return -1;
189
190         lp->con_offset += 2;
191
192         return lp->tab->empty;
193 }
194
195 static void get_obj_val(struct tab_lp* lp, mpq_t *F)
196 {
197         isl_int_neg(mpq_numref(*F), lp->opt);
198         isl_int_set(mpq_denref(*F), lp->opt_denom);
199 }
200
201 static void delete_lp(struct tab_lp *lp)
202 {
203         if (!lp)
204                 return;
205
206         isl_int_clear(lp->opt);
207         isl_int_clear(lp->opt_denom);
208         isl_int_clear(lp->tmp);
209         isl_int_clear(lp->tmp2);
210         isl_vec_free(lp->row);
211         free(lp->stack);
212         isl_tab_free(lp->tab);
213         isl_ctx_deref(lp->ctx);
214         free(lp);
215 }
216
217 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)
218 {
219         lp->stack[lp->neq] = isl_tab_snap(lp->tab);
220
221         isl_int_set_si(lp->row->el[0], 0);
222         isl_seq_cpy(lp->row->el + 1, row, lp->dim);
223         isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);
224
225         if (isl_tab_add_valid_eq(lp->tab, lp->row->el) < 0)
226                 return -1;
227
228         return lp->neq++;
229 }
230
231 static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha)
232 {
233         row += lp->con_offset;
234         isl_int_neg(mpq_numref(*alpha), lp->tab->dual->el[1 + row]);
235         isl_int_set(mpq_denref(*alpha), lp->tab->dual->el[0]);
236 }
237
238 static int del_lp_row(struct tab_lp *lp)
239 {
240         lp->neq--;
241         return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);
242 }