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 <isl_ctx_private.h>
11#include <isl_map_private.h>
12#include <isl/lp.h>
13#include "isl_lp_piplib.h"
14#include <isl/seq.h>
15#include "isl_tab.h"
16#include <isl_options_private.h>
17#include <isl_local_space_private.h>
18#include <isl_aff_private.h>
19#include <isl_mat_private.h>
20#include <isl_val_private.h>
21
22enum isl_lp_result isl_tab_solve_lp(struct isl_basic_map *bmap, int maximize,
23				      isl_int *f, isl_int denom, isl_int *opt,
24				      isl_int *opt_denom,
25				      struct isl_vec **sol)
26{
27	struct isl_tab *tab;
28	enum isl_lp_result res;
29	unsigned dim = isl_basic_map_total_dim(bmap);
30
31	if (maximize)
32		isl_seq_neg(f, f, 1 + dim);
33
34	bmap = isl_basic_map_gauss(bmap, NULL);
35	tab = isl_tab_from_basic_map(bmap, 0);
36	res = isl_tab_min(tab, f, denom, opt, opt_denom, 0);
37	if (res == isl_lp_ok && sol) {
38		*sol = isl_tab_get_sample_value(tab);
39		if (!*sol)
40			res = isl_lp_error;
41	}
42	isl_tab_free(tab);
43
44	if (maximize)
45		isl_seq_neg(f, f, 1 + dim);
46	if (maximize && opt)
47		isl_int_neg(*opt, *opt);
48
49	return res;
50}
51
52/* Given a basic map "bmap" and an affine combination of the variables "f"
53 * with denominator "denom", set *opt / *opt_denom to the minimal
54 * (or maximal if "maximize" is true) value attained by f/d over "bmap",
55 * assuming the basic map is not empty and the expression cannot attain
56 * arbitrarily small (or large) values.
57 * If opt_denom is NULL, then *opt is rounded up (or down)
58 * to the nearest integer.
59 * The return value reflects the nature of the result (empty, unbounded,
60 * minmimal or maximal value returned in *opt).
61 */
62enum isl_lp_result isl_basic_map_solve_lp(struct isl_basic_map *bmap, int max,
63				      isl_int *f, isl_int d, isl_int *opt,
64				      isl_int *opt_denom,
65				      struct isl_vec **sol)
66{
67	if (sol)
68		*sol = NULL;
69
70	if (!bmap)
71		return isl_lp_error;
72
73	switch (bmap->ctx->opt->lp_solver) {
74	case ISL_LP_PIP:
75		return isl_pip_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
76	case ISL_LP_TAB:
77		return isl_tab_solve_lp(bmap, max, f, d, opt, opt_denom, sol);
78	default:
79		return isl_lp_error;
80	}
81}
82
83enum isl_lp_result isl_basic_set_solve_lp(struct isl_basic_set *bset, int max,
84				      isl_int *f, isl_int d, isl_int *opt,
85				      isl_int *opt_denom,
86				      struct isl_vec **sol)
87{
88	return isl_basic_map_solve_lp((struct isl_basic_map *)bset, max,
89					f, d, opt, opt_denom, sol);
90}
91
92enum isl_lp_result isl_map_solve_lp(__isl_keep isl_map *map, int max,
93				      isl_int *f, isl_int d, isl_int *opt,
94				      isl_int *opt_denom,
95				      struct isl_vec **sol)
96{
97	int i;
98	isl_int o;
99	isl_int t;
100	isl_int opt_i;
101	isl_int opt_denom_i;
102	enum isl_lp_result res;
103	int max_div;
104	isl_vec *v = NULL;
105
106	if (!map)
107		return isl_lp_error;
108	if (map->n == 0)
109		return isl_lp_empty;
110
111	max_div = 0;
112	for (i = 0; i < map->n; ++i)
113		if (map->p[i]->n_div > max_div)
114			max_div = map->p[i]->n_div;
115	if (max_div > 0) {
116		unsigned total = isl_space_dim(map->dim, isl_dim_all);
117		v = isl_vec_alloc(map->ctx, 1 + total + max_div);
118		if (!v)
119			return isl_lp_error;
120		isl_seq_cpy(v->el, f, 1 + total);
121		isl_seq_clr(v->el + 1 + total, max_div);
122		f = v->el;
123	}
124
125	if (!opt && map->n > 1 && sol) {
126		isl_int_init(o);
127		opt = &o;
128	}
129	if (map->n > 0)
130		isl_int_init(opt_i);
131	if (map->n > 0 && opt_denom) {
132		isl_int_init(opt_denom_i);
133		isl_int_init(t);
134	}
135
136	res = isl_basic_map_solve_lp(map->p[0], max, f, d,
137					opt, opt_denom, sol);
138	if (res == isl_lp_error || res == isl_lp_unbounded)
139		goto done;
140
141	if (sol)
142		*sol = NULL;
143
144	for (i = 1; i < map->n; ++i) {
145		isl_vec *sol_i = NULL;
146		enum isl_lp_result res_i;
147		int better;
148
149		res_i = isl_basic_map_solve_lp(map->p[i], max, f, d,
150					    &opt_i,
151					    opt_denom ? &opt_denom_i : NULL,
152					    sol ? &sol_i : NULL);
153		if (res_i == isl_lp_error || res_i == isl_lp_unbounded) {
154			res = res_i;
155			goto done;
156		}
157		if (res_i == isl_lp_empty)
158			continue;
159		if (res == isl_lp_empty) {
160			better = 1;
161		} else if (!opt_denom) {
162			if (max)
163				better = isl_int_gt(opt_i, *opt);
164			else
165				better = isl_int_lt(opt_i, *opt);
166		} else {
167			isl_int_mul(t, opt_i, *opt_denom);
168			isl_int_submul(t, *opt, opt_denom_i);
169			if (max)
170				better = isl_int_is_pos(t);
171			else
172				better = isl_int_is_neg(t);
173		}
174		if (better) {
175			res = res_i;
176			if (opt)
177				isl_int_set(*opt, opt_i);
178			if (opt_denom)
179				isl_int_set(*opt_denom, opt_denom_i);
180			if (sol) {
181				isl_vec_free(*sol);
182				*sol = sol_i;
183			}
184		} else
185			isl_vec_free(sol_i);
186	}
187
188done:
189	isl_vec_free(v);
190	if (map->n > 0 && opt_denom) {
191		isl_int_clear(opt_denom_i);
192		isl_int_clear(t);
193	}
194	if (map->n > 0)
195		isl_int_clear(opt_i);
196	if (opt == &o)
197		isl_int_clear(o);
198	return res;
199}
200
201enum isl_lp_result isl_set_solve_lp(__isl_keep isl_set *set, int max,
202				      isl_int *f, isl_int d, isl_int *opt,
203				      isl_int *opt_denom,
204				      struct isl_vec **sol)
205{
206	return isl_map_solve_lp((struct isl_map *)set, max,
207					f, d, opt, opt_denom, sol);
208}
209
210/* Return the optimal (rational) value of "obj" over "bset", assuming
211 * that "obj" and "bset" have aligned parameters and divs.
212 * If "max" is set, then the maximal value is computed.
213 * Otherwise, the minimal value is computed.
214 *
215 * Return infinity or negative infinity if the optimal value is unbounded and
216 * NaN if "bset" is empty.
217 *
218 * Call isl_basic_set_solve_lp and translate the results.
219 */
220static __isl_give isl_val *basic_set_opt_lp(
221	__isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
222{
223	isl_ctx *ctx;
224	isl_val *res;
225	enum isl_lp_result lp_res;
226
227	if (!bset || !obj)
228		return NULL;
229
230	ctx = isl_aff_get_ctx(obj);
231	res = isl_val_alloc(ctx);
232	if (!res)
233		return NULL;
234	lp_res = isl_basic_set_solve_lp(bset, max, obj->v->el + 1,
235					obj->v->el[0], &res->n, &res->d, NULL);
236	if (lp_res == isl_lp_ok)
237		return isl_val_normalize(res);
238	isl_val_free(res);
239	if (lp_res == isl_lp_error)
240		return NULL;
241	if (lp_res == isl_lp_empty)
242		return isl_val_nan(ctx);
243	if (max)
244		return isl_val_infty(ctx);
245	else
246		return isl_val_neginfty(ctx);
247}
248
249/* Return the optimal (rational) value of "obj" over "bset", assuming
250 * that "obj" and "bset" have aligned parameters.
251 * If "max" is set, then the maximal value is computed.
252 * Otherwise, the minimal value is computed.
253 *
254 * Return infinity or negative infinity if the optimal value is unbounded and
255 * NaN if "bset" is empty.
256 *
257 * Align the divs of "bset" and "obj" and call basic_set_opt_lp.
258 */
259static __isl_give isl_val *isl_basic_set_opt_lp_val_aligned(
260	__isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
261{
262	int *exp1 = NULL;
263	int *exp2 = NULL;
264	isl_ctx *ctx;
265	isl_mat *bset_div = NULL;
266	isl_mat *div = NULL;
267	isl_val *res;
268	int bset_n_div, obj_n_div;
269
270	if (!bset || !obj)
271		return NULL;
272
273	ctx = isl_aff_get_ctx(obj);
274	if (!isl_space_is_equal(bset->dim, obj->ls->dim))
275		isl_die(ctx, isl_error_invalid,
276			"spaces don't match", return NULL);
277
278	bset_n_div = isl_basic_set_dim(bset, isl_dim_div);
279	obj_n_div = isl_aff_dim(obj, isl_dim_div);
280	if (bset_n_div == 0 && obj_n_div == 0)
281		return basic_set_opt_lp(bset, max, obj);
282
283	bset = isl_basic_set_copy(bset);
284	obj = isl_aff_copy(obj);
285
286	bset_div = isl_basic_set_get_divs(bset);
287	exp1 = isl_alloc_array(ctx, int, bset_n_div);
288	exp2 = isl_alloc_array(ctx, int, obj_n_div);
289	if (!bset_div || (bset_n_div && !exp1) || (obj_n_div && !exp2))
290		goto error;
291
292	div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2);
293
294	bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1);
295	obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2);
296
297	res = basic_set_opt_lp(bset, max, obj);
298
299	isl_mat_free(bset_div);
300	isl_mat_free(div);
301	free(exp1);
302	free(exp2);
303	isl_basic_set_free(bset);
304	isl_aff_free(obj);
305
306	return res;
307error:
308	isl_mat_free(div);
309	isl_mat_free(bset_div);
310	free(exp1);
311	free(exp2);
312	isl_basic_set_free(bset);
313	isl_aff_free(obj);
314	return NULL;
315}
316
317/* Return the optimal (rational) value of "obj" over "bset".
318 * If "max" is set, then the maximal value is computed.
319 * Otherwise, the minimal value is computed.
320 *
321 * Return infinity or negative infinity if the optimal value is unbounded and
322 * NaN if "bset" is empty.
323 */
324static __isl_give isl_val *isl_basic_set_opt_lp_val(
325	__isl_keep isl_basic_set *bset, int max, __isl_keep isl_aff *obj)
326{
327	isl_val *res;
328
329	if (!bset || !obj)
330		return NULL;
331
332	if (isl_space_match(bset->dim, isl_dim_param,
333			    obj->ls->dim, isl_dim_param))
334		return isl_basic_set_opt_lp_val_aligned(bset, max, obj);
335
336	bset = isl_basic_set_copy(bset);
337	obj = isl_aff_copy(obj);
338	bset = isl_basic_set_align_params(bset, isl_aff_get_domain_space(obj));
339	obj = isl_aff_align_params(obj, isl_basic_set_get_space(bset));
340
341	res = isl_basic_set_opt_lp_val_aligned(bset, max, obj);
342
343	isl_basic_set_free(bset);
344	isl_aff_free(obj);
345
346	return res;
347}
348
349/* Return the minimal (rational) value of "obj" over "bset".
350 *
351 * Return negative infinity if the minimal value is unbounded and
352 * NaN if "bset" is empty.
353 */
354__isl_give isl_val *isl_basic_set_min_lp_val(__isl_keep isl_basic_set *bset,
355	__isl_keep isl_aff *obj)
356{
357	return isl_basic_set_opt_lp_val(bset, 0, obj);
358}
359
360/* Return the maximal (rational) value of "obj" over "bset".
361 *
362 * Return infinity if the maximal value is unbounded and
363 * NaN if "bset" is empty.
364 */
365__isl_give isl_val *isl_basic_set_max_lp_val(__isl_keep isl_basic_set *bset,
366	__isl_keep isl_aff *obj)
367{
368	return isl_basic_set_opt_lp_val(bset, 1, obj);
369}
370