add isl_union_map_intersect_params
[isl.git] / isl_ilp.c
blobabfac265ace6bd6a67b34bc3e2ebf726206a6b81
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
10 #include <isl_ctx_private.h>
11 #include <isl_map_private.h>
12 #include <isl/ilp.h>
13 #include "isl_sample.h"
14 #include <isl/seq.h>
15 #include "isl_equalities.h"
16 #include <isl_aff_private.h>
17 #include <isl_local_space_private.h>
18 #include <isl_mat_private.h>
20 /* Given a basic set "bset", construct a basic set U such that for
21 * each element x in U, the whole unit box positioned at x is inside
22 * the given basic set.
23 * Note that U may not contain all points that satisfy this property.
25 * We simply add the sum of all negative coefficients to the constant
26 * term. This ensures that if x satisfies the resulting constraints,
27 * then x plus any sum of unit vectors satisfies the original constraints.
29 static struct isl_basic_set *unit_box_base_points(struct isl_basic_set *bset)
31 int i, j, k;
32 struct isl_basic_set *unit_box = NULL;
33 unsigned total;
35 if (!bset)
36 goto error;
38 if (bset->n_eq != 0) {
39 unit_box = isl_basic_set_empty_like(bset);
40 isl_basic_set_free(bset);
41 return unit_box;
44 total = isl_basic_set_total_dim(bset);
45 unit_box = isl_basic_set_alloc_space(isl_basic_set_get_space(bset),
46 0, 0, bset->n_ineq);
48 for (i = 0; i < bset->n_ineq; ++i) {
49 k = isl_basic_set_alloc_inequality(unit_box);
50 if (k < 0)
51 goto error;
52 isl_seq_cpy(unit_box->ineq[k], bset->ineq[i], 1 + total);
53 for (j = 0; j < total; ++j) {
54 if (isl_int_is_nonneg(unit_box->ineq[k][1 + j]))
55 continue;
56 isl_int_add(unit_box->ineq[k][0],
57 unit_box->ineq[k][0], unit_box->ineq[k][1 + j]);
61 isl_basic_set_free(bset);
62 return unit_box;
63 error:
64 isl_basic_set_free(bset);
65 isl_basic_set_free(unit_box);
66 return NULL;
69 /* Find an integer point in "bset", preferably one that is
70 * close to minimizing "f".
72 * We first check if we can easily put unit boxes inside bset.
73 * If so, we take the best base point of any of the unit boxes we can find
74 * and round it up to the nearest integer.
75 * If not, we simply pick any integer point in "bset".
77 static struct isl_vec *initial_solution(struct isl_basic_set *bset, isl_int *f)
79 enum isl_lp_result res;
80 struct isl_basic_set *unit_box;
81 struct isl_vec *sol;
83 unit_box = unit_box_base_points(isl_basic_set_copy(bset));
85 res = isl_basic_set_solve_lp(unit_box, 0, f, bset->ctx->one,
86 NULL, NULL, &sol);
87 if (res == isl_lp_ok) {
88 isl_basic_set_free(unit_box);
89 return isl_vec_ceil(sol);
92 isl_basic_set_free(unit_box);
94 return isl_basic_set_sample_vec(isl_basic_set_copy(bset));
97 /* Restrict "bset" to those points with values for f in the interval [l, u].
99 static struct isl_basic_set *add_bounds(struct isl_basic_set *bset,
100 isl_int *f, isl_int l, isl_int u)
102 int k;
103 unsigned total;
105 total = isl_basic_set_total_dim(bset);
106 bset = isl_basic_set_extend_constraints(bset, 0, 2);
108 k = isl_basic_set_alloc_inequality(bset);
109 if (k < 0)
110 goto error;
111 isl_seq_cpy(bset->ineq[k], f, 1 + total);
112 isl_int_sub(bset->ineq[k][0], bset->ineq[k][0], l);
114 k = isl_basic_set_alloc_inequality(bset);
115 if (k < 0)
116 goto error;
117 isl_seq_neg(bset->ineq[k], f, 1 + total);
118 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], u);
120 return bset;
121 error:
122 isl_basic_set_free(bset);
123 return NULL;
126 /* Find an integer point in "bset" that minimizes f (in any) such that
127 * the value of f lies inside the interval [l, u].
128 * Return this integer point if it can be found.
129 * Otherwise, return sol.
131 * We perform a number of steps until l > u.
132 * In each step, we look for an integer point with value in either
133 * the whole interval [l, u] or half of the interval [l, l+floor(u-l-1/2)].
134 * The choice depends on whether we have found an integer point in the
135 * previous step. If so, we look for the next point in half of the remaining
136 * interval.
137 * If we find a point, the current solution is updated and u is set
138 * to its value minus 1.
139 * If no point can be found, we update l to the upper bound of the interval
140 * we checked (u or l+floor(u-l-1/2)) plus 1.
142 static struct isl_vec *solve_ilp_search(struct isl_basic_set *bset,
143 isl_int *f, isl_int *opt, struct isl_vec *sol, isl_int l, isl_int u)
145 isl_int tmp;
146 int divide = 1;
148 isl_int_init(tmp);
150 while (isl_int_le(l, u)) {
151 struct isl_basic_set *slice;
152 struct isl_vec *sample;
154 if (!divide)
155 isl_int_set(tmp, u);
156 else {
157 isl_int_sub(tmp, u, l);
158 isl_int_fdiv_q_ui(tmp, tmp, 2);
159 isl_int_add(tmp, tmp, l);
161 slice = add_bounds(isl_basic_set_copy(bset), f, l, tmp);
162 sample = isl_basic_set_sample_vec(slice);
163 if (!sample) {
164 isl_vec_free(sol);
165 sol = NULL;
166 break;
168 if (sample->size > 0) {
169 isl_vec_free(sol);
170 sol = sample;
171 isl_seq_inner_product(f, sol->el, sol->size, opt);
172 isl_int_sub_ui(u, *opt, 1);
173 divide = 1;
174 } else {
175 isl_vec_free(sample);
176 if (!divide)
177 break;
178 isl_int_add_ui(l, tmp, 1);
179 divide = 0;
183 isl_int_clear(tmp);
185 return sol;
188 /* Find an integer point in "bset" that minimizes f (if any).
189 * If sol_p is not NULL then the integer point is returned in *sol_p.
190 * The optimal value of f is returned in *opt.
192 * The algorithm maintains a currently best solution and an interval [l, u]
193 * of values of f for which integer solutions could potentially still be found.
194 * The initial value of the best solution so far is any solution.
195 * The initial value of l is minimal value of f over the rationals
196 * (rounded up to the nearest integer).
197 * The initial value of u is the value of f at the initial solution minus 1.
199 * We then call solve_ilp_search to perform a binary search on the interval.
201 static enum isl_lp_result solve_ilp(struct isl_basic_set *bset,
202 isl_int *f, isl_int *opt,
203 struct isl_vec **sol_p)
205 enum isl_lp_result res;
206 isl_int l, u;
207 struct isl_vec *sol;
209 res = isl_basic_set_solve_lp(bset, 0, f, bset->ctx->one,
210 opt, NULL, &sol);
211 if (res == isl_lp_ok && isl_int_is_one(sol->el[0])) {
212 if (sol_p)
213 *sol_p = sol;
214 else
215 isl_vec_free(sol);
216 return isl_lp_ok;
218 isl_vec_free(sol);
219 if (res == isl_lp_error || res == isl_lp_empty)
220 return res;
222 sol = initial_solution(bset, f);
223 if (!sol)
224 return isl_lp_error;
225 if (sol->size == 0) {
226 isl_vec_free(sol);
227 return isl_lp_empty;
229 if (res == isl_lp_unbounded) {
230 isl_vec_free(sol);
231 return isl_lp_unbounded;
234 isl_int_init(l);
235 isl_int_init(u);
237 isl_int_set(l, *opt);
239 isl_seq_inner_product(f, sol->el, sol->size, opt);
240 isl_int_sub_ui(u, *opt, 1);
242 sol = solve_ilp_search(bset, f, opt, sol, l, u);
243 if (!sol)
244 res = isl_lp_error;
246 isl_int_clear(l);
247 isl_int_clear(u);
249 if (sol_p)
250 *sol_p = sol;
251 else
252 isl_vec_free(sol);
254 return res;
257 static enum isl_lp_result solve_ilp_with_eq(struct isl_basic_set *bset, int max,
258 isl_int *f, isl_int *opt,
259 struct isl_vec **sol_p)
261 unsigned dim;
262 enum isl_lp_result res;
263 struct isl_mat *T = NULL;
264 struct isl_vec *v;
266 bset = isl_basic_set_copy(bset);
267 dim = isl_basic_set_total_dim(bset);
268 v = isl_vec_alloc(bset->ctx, 1 + dim);
269 if (!v)
270 goto error;
271 isl_seq_cpy(v->el, f, 1 + dim);
272 bset = isl_basic_set_remove_equalities(bset, &T, NULL);
273 v = isl_vec_mat_product(v, isl_mat_copy(T));
274 if (!v)
275 goto error;
276 res = isl_basic_set_solve_ilp(bset, max, v->el, opt, sol_p);
277 isl_vec_free(v);
278 if (res == isl_lp_ok && sol_p) {
279 *sol_p = isl_mat_vec_product(T, *sol_p);
280 if (!*sol_p)
281 res = isl_lp_error;
282 } else
283 isl_mat_free(T);
284 isl_basic_set_free(bset);
285 return res;
286 error:
287 isl_mat_free(T);
288 isl_basic_set_free(bset);
289 return isl_lp_error;
292 /* Find an integer point in "bset" that minimizes (or maximizes if max is set)
293 * f (if any).
294 * If sol_p is not NULL then the integer point is returned in *sol_p.
295 * The optimal value of f is returned in *opt.
297 * If there is any equality among the points in "bset", then we first
298 * project it out. Otherwise, we continue with solve_ilp above.
300 enum isl_lp_result isl_basic_set_solve_ilp(struct isl_basic_set *bset, int max,
301 isl_int *f, isl_int *opt,
302 struct isl_vec **sol_p)
304 unsigned dim;
305 enum isl_lp_result res;
307 if (!bset)
308 return isl_lp_error;
309 if (sol_p)
310 *sol_p = NULL;
312 isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
314 if (isl_basic_set_plain_is_empty(bset))
315 return isl_lp_empty;
317 if (bset->n_eq)
318 return solve_ilp_with_eq(bset, max, f, opt, sol_p);
320 dim = isl_basic_set_total_dim(bset);
322 if (max)
323 isl_seq_neg(f, f, 1 + dim);
325 res = solve_ilp(bset, f, opt, sol_p);
327 if (max) {
328 isl_seq_neg(f, f, 1 + dim);
329 isl_int_neg(*opt, *opt);
332 return res;
333 error:
334 isl_basic_set_free(bset);
335 return isl_lp_error;
338 static enum isl_lp_result basic_set_opt(__isl_keep isl_basic_set *bset, int max,
339 __isl_keep isl_aff *obj, isl_int *opt)
341 enum isl_lp_result res;
343 if (!obj)
344 return isl_lp_error;
345 bset = isl_basic_set_copy(bset);
346 bset = isl_basic_set_underlying_set(bset);
347 res = isl_basic_set_solve_ilp(bset, max, obj->v->el + 1, opt, NULL);
348 isl_basic_set_free(bset);
349 return res;
352 static __isl_give isl_mat *extract_divs(__isl_keep isl_basic_set *bset)
354 int i;
355 isl_ctx *ctx = isl_basic_set_get_ctx(bset);
356 isl_mat *div;
358 div = isl_mat_alloc(ctx, bset->n_div,
359 1 + 1 + isl_basic_set_total_dim(bset));
360 if (!div)
361 return NULL;
363 for (i = 0; i < bset->n_div; ++i)
364 isl_seq_cpy(div->row[i], bset->div[i], div->n_col);
366 return div;
369 enum isl_lp_result isl_basic_set_opt(__isl_keep isl_basic_set *bset, int max,
370 __isl_keep isl_aff *obj, isl_int *opt)
372 int *exp1 = NULL;
373 int *exp2 = NULL;
374 isl_ctx *ctx;
375 isl_mat *bset_div = NULL;
376 isl_mat *div = NULL;
377 enum isl_lp_result res;
379 if (!bset || !obj)
380 return isl_lp_error;
382 ctx = isl_aff_get_ctx(obj);
383 if (!isl_space_is_equal(bset->dim, obj->ls->dim))
384 isl_die(ctx, isl_error_invalid,
385 "spaces don't match", return isl_lp_error);
386 if (!isl_int_is_one(obj->v->el[0]))
387 isl_die(ctx, isl_error_unsupported,
388 "expecting integer affine expression",
389 return isl_lp_error);
391 if (bset->n_div == 0 && obj->ls->div->n_row == 0)
392 return basic_set_opt(bset, max, obj, opt);
394 bset = isl_basic_set_copy(bset);
395 obj = isl_aff_copy(obj);
397 bset_div = extract_divs(bset);
398 exp1 = isl_alloc_array(ctx, int, bset_div->n_row);
399 exp2 = isl_alloc_array(ctx, int, obj->ls->div->n_row);
400 if (!bset_div || !exp1 || !exp2)
401 goto error;
403 div = isl_merge_divs(bset_div, obj->ls->div, exp1, exp2);
405 bset = isl_basic_set_expand_divs(bset, isl_mat_copy(div), exp1);
406 obj = isl_aff_expand_divs(obj, isl_mat_copy(div), exp2);
408 res = basic_set_opt(bset, max, obj, opt);
410 isl_mat_free(bset_div);
411 isl_mat_free(div);
412 free(exp1);
413 free(exp2);
414 isl_basic_set_free(bset);
415 isl_aff_free(obj);
417 return res;
418 error:
419 isl_mat_free(div);
420 isl_mat_free(bset_div);
421 free(exp1);
422 free(exp2);
423 isl_basic_set_free(bset);
424 isl_aff_free(obj);
425 return isl_lp_error;
428 /* Compute the minimum (maximum if max is set) of the integer affine
429 * expression obj over the points in set and put the result in *opt.
431 enum isl_lp_result isl_set_opt(__isl_keep isl_set *set, int max,
432 __isl_keep isl_aff *obj, isl_int *opt)
434 int i;
435 enum isl_lp_result res;
436 int empty = 1;
437 isl_int opt_i;
439 if (!set || !obj)
440 return isl_lp_error;
441 if (set->n == 0)
442 return isl_lp_empty;
444 res = isl_basic_set_opt(set->p[0], max, obj, opt);
445 if (res == isl_lp_error || res == isl_lp_unbounded)
446 return res;
447 if (set->n == 1)
448 return res;
449 if (res == isl_lp_ok)
450 empty = 0;
452 isl_int_init(opt_i);
453 for (i = 1; i < set->n; ++i) {
454 res = isl_basic_set_opt(set->p[i], max, obj, &opt_i);
455 if (res == isl_lp_error || res == isl_lp_unbounded) {
456 isl_int_clear(opt_i);
457 return res;
459 if (res == isl_lp_ok)
460 empty = 0;
461 if (isl_int_gt(opt_i, *opt))
462 isl_int_set(*opt, opt_i);
464 isl_int_clear(opt_i);
466 return empty ? isl_lp_empty : isl_lp_ok;
469 enum isl_lp_result isl_basic_set_max(__isl_keep isl_basic_set *bset,
470 __isl_keep isl_aff *obj, isl_int *opt)
472 return isl_basic_set_opt(bset, 1, obj, opt);
475 enum isl_lp_result isl_set_max(__isl_keep isl_set *set,
476 __isl_keep isl_aff *obj, isl_int *opt)
478 return isl_set_opt(set, 1, obj, opt);
481 enum isl_lp_result isl_set_min(__isl_keep isl_set *set,
482 __isl_keep isl_aff *obj, isl_int *opt)
484 return isl_set_opt(set, 0, obj, opt);