nexmon – Blame information for rev 1

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