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
15struct 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
33static struct tab_lp *init_lp(struct isl_tab *tab);
34static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
35static int solve_lp(struct tab_lp *lp);
36static void get_obj_val(struct tab_lp* lp, mpq_t *F);
37static void delete_lp(struct tab_lp *lp);
38static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
39static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha);
40static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;
41static 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 */
71static 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
88static 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;
122error:
123	delete_lp(lp);
124	return NULL;
125}
126
127static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
128{
129	lp->obj = row;
130}
131
132static 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 */
170static 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
196static 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
202static 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
218static 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
232static 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
239static int del_lp_row(struct tab_lp *lp)
240{
241	lp->neq--;
242	return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);
243}
244