2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the MIT license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10 #include <isl_ctx_private.h>
11 #include <isl_map_private.h>
13 #include "isl_sample.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>
19 #include <isl_val_private.h>
21 /* Given a basic set "bset", construct a basic set U such that for
22 * each element x in U, the whole unit box positioned at x is inside
23 * the given basic set.
24 * Note that U may not contain all points that satisfy this property.
26 * We simply add the sum of all negative coefficients to the constant
27 * term. This ensures that if x satisfies the resulting constraints,
28 * then x plus any sum of unit vectors satisfies the original constraints.
30 static struct isl_basic_set
*unit_box_base_points(struct isl_basic_set
*bset
)
33 struct isl_basic_set
*unit_box
= NULL
;
39 if (bset
->n_eq
!= 0) {
40 unit_box
= isl_basic_set_empty_like(bset
);
41 isl_basic_set_free(bset
);
45 total
= isl_basic_set_total_dim(bset
);
46 unit_box
= isl_basic_set_alloc_space(isl_basic_set_get_space(bset
),
49 for (i
= 0; i
< bset
->n_ineq
; ++i
) {
50 k
= isl_basic_set_alloc_inequality(unit_box
);
53 isl_seq_cpy(unit_box
->ineq
[k
], bset
->ineq
[i
], 1 + total
);
54 for (j
= 0; j
< total
; ++j
) {
55 if (isl_int_is_nonneg(unit_box
->ineq
[k
][1 + j
]))
57 isl_int_add(unit_box
->ineq
[k
][0],
58 unit_box
->ineq
[k
][0], unit_box
->ineq
[k
][1 + j
]);
62 isl_basic_set_free(bset
);
65 isl_basic_set_free(bset
);
66 isl_basic_set_free(unit_box
);
70 /* Find an integer point in "bset", preferably one that is
71 * close to minimizing "f".
73 * We first check if we can easily put unit boxes inside bset.
74 * If so, we take the best base point of any of the unit boxes we can find
75 * and round it up to the nearest integer.
76 * If not, we simply pick any integer point in "bset".
78 static struct isl_vec
*initial_solution(struct isl_basic_set
*bset
, isl_int
*f
)
80 enum isl_lp_result res
;
81 struct isl_basic_set
*unit_box
;
84 unit_box
= unit_box_base_points(isl_basic_set_copy(bset
));
86 res
= isl_basic_set_solve_lp(unit_box
, 0, f
, bset
->ctx
->one
,
88 if (res
== isl_lp_ok
) {
89 isl_basic_set_free(unit_box
);
90 return isl_vec_ceil(sol
);
93 isl_basic_set_free(unit_box
);
95 return isl_basic_set_sample_vec(isl_basic_set_copy(bset
));
98 /* Restrict "bset" to those points with values for f in the interval [l, u].
100 static struct isl_basic_set
*add_bounds(struct isl_basic_set
*bset
,
101 isl_int
*f
, isl_int l
, isl_int u
)
106 total
= isl_basic_set_total_dim(bset
);
107 bset
= isl_basic_set_extend_constraints(bset
, 0, 2);
109 k
= isl_basic_set_alloc_inequality(bset
);
112 isl_seq_cpy(bset
->ineq
[k
], f
, 1 + total
);
113 isl_int_sub(bset
->ineq
[k
][0], bset
->ineq
[k
][0], l
);
115 k
= isl_basic_set_alloc_inequality(bset
);
118 isl_seq_neg(bset
->ineq
[k
], f
, 1 + total
);
119 isl_int_add(bset
->ineq
[k
][0], bset
->ineq
[k
][0], u
);
123 isl_basic_set_free(bset
);
127 /* Find an integer point in "bset" that minimizes f (in any) such that
128 * the value of f lies inside the interval [l, u].
129 * Return this integer point if it can be found.
130 * Otherwise, return sol.
132 * We perform a number of steps until l > u.
133 * In each step, we look for an integer point with value in either
134 * the whole interval [l, u] or half of the interval [l, l+floor(u-l-1/2)].
135 * The choice depends on whether we have found an integer point in the
136 * previous step. If so, we look for the next point in half of the remaining
138 * If we find a point, the current solution is updated and u is set
139 * to its value minus 1.
140 * If no point can be found, we update l to the upper bound of the interval
141 * we checked (u or l+floor(u-l-1/2)) plus 1.
143 static struct isl_vec
*solve_ilp_search(struct isl_basic_set
*bset
,
144 isl_int
*f
, isl_int
*opt
, struct isl_vec
*sol
, isl_int l
, isl_int u
)
151 while (isl_int_le(l
, u
)) {
152 struct isl_basic_set
*slice
;
153 struct isl_vec
*sample
;
158 isl_int_sub(tmp
, u
, l
);
159 isl_int_fdiv_q_ui(tmp
, tmp
, 2);
160 isl_int_add(tmp
, tmp
, l
);
162 slice
= add_bounds(isl_basic_set_copy(bset
), f
, l
, tmp
);
163 sample
= isl_basic_set_sample_vec(slice
);
169 if (sample
->size
> 0) {
172 isl_seq_inner_product(f
, sol
->el
, sol
->size
, opt
);
173 isl_int_sub_ui(u
, *opt
, 1);
176 isl_vec_free(sample
);
179 isl_int_add_ui(l
, tmp
, 1);
189 /* Find an integer point in "bset" that minimizes f (if any).
190 * If sol_p is not NULL then the integer point is returned in *sol_p.
191 * The optimal value of f is returned in *opt.
193 * The algorithm maintains a currently best solution and an interval [l, u]
194 * of values of f for which integer solutions could potentially still be found.
195 * The initial value of the best solution so far is any solution.
196 * The initial value of l is minimal value of f over the rationals
197 * (rounded up to the nearest integer).
198 * The initial value of u is the value of f at the initial solution minus 1.
200 * We then call solve_ilp_search to perform a binary search on the interval.
202 static enum isl_lp_result
solve_ilp(struct isl_basic_set
*bset
,
203 isl_int
*f
, isl_int
*opt
,
204 struct isl_vec
**sol_p
)
206 enum isl_lp_result res
;
210 res
= isl_basic_set_solve_lp(bset
, 0, f
, bset
->ctx
->one
,
212 if (res
== isl_lp_ok
&& isl_int_is_one(sol
->el
[0])) {
220 if (res
== isl_lp_error
|| res
== isl_lp_empty
)
223 sol
= initial_solution(bset
, f
);
226 if (sol
->size
== 0) {
230 if (res
== isl_lp_unbounded
) {
232 return isl_lp_unbounded
;
238 isl_int_set(l
, *opt
);
240 isl_seq_inner_product(f
, sol
->el
, sol
->size
, opt
);
241 isl_int_sub_ui(u
, *opt
, 1);
243 sol
= solve_ilp_search(bset
, f
, opt
, sol
, l
, u
);
258 static enum isl_lp_result
solve_ilp_with_eq(struct isl_basic_set
*bset
, int max
,
259 isl_int
*f
, isl_int
*opt
,
260 struct isl_vec
**sol_p
)
263 enum isl_lp_result res
;
264 struct isl_mat
*T
= NULL
;
267 bset
= isl_basic_set_copy(bset
);
268 dim
= isl_basic_set_total_dim(bset
);
269 v
= isl_vec_alloc(bset
->ctx
, 1 + dim
);
272 isl_seq_cpy(v
->el
, f
, 1 + dim
);
273 bset
= isl_basic_set_remove_equalities(bset
, &T
, NULL
);
274 v
= isl_vec_mat_product(v
, isl_mat_copy(T
));
277 res
= isl_basic_set_solve_ilp(bset
, max
, v
->el
, opt
, sol_p
);
279 if (res
== isl_lp_ok
&& sol_p
) {
280 *sol_p
= isl_mat_vec_product(T
, *sol_p
);
285 isl_basic_set_free(bset
);
289 isl_basic_set_free(bset
);
293 /* Find an integer point in "bset" that minimizes (or maximizes if max is set)
295 * If sol_p is not NULL then the integer point is returned in *sol_p.
296 * The optimal value of f is returned in *opt.
298 * If there is any equality among the points in "bset", then we first
299 * project it out. Otherwise, we continue with solve_ilp above.
301 enum isl_lp_result
isl_basic_set_solve_ilp(struct isl_basic_set
*bset
, int max
,
302 isl_int
*f
, isl_int
*opt
,
303 struct isl_vec
**sol_p
)
306 enum isl_lp_result res
;
313 isl_assert(bset
->ctx
, isl_basic_set_n_param(bset
) == 0, goto error
);
315 if (isl_basic_set_plain_is_empty(bset
))
319 return solve_ilp_with_eq(bset
, max
, f
, opt
, sol_p
);
321 dim
= isl_basic_set_total_dim(bset
);
324 isl_seq_neg(f
, f
, 1 + dim
);
326 res
= solve_ilp(bset
, f
, opt
, sol_p
);
329 isl_seq_neg(f
, f
, 1 + dim
);
330 isl_int_neg(*opt
, *opt
);
335 isl_basic_set_free(bset
);
339 static enum isl_lp_result
basic_set_opt(__isl_keep isl_basic_set
*bset
, int max
,
340 __isl_keep isl_aff
*obj
, isl_int
*opt
)
342 enum isl_lp_result res
;
346 bset
= isl_basic_set_copy(bset
);
347 bset
= isl_basic_set_underlying_set(bset
);
348 res
= isl_basic_set_solve_ilp(bset
, max
, obj
->v
->el
+ 1, opt
, NULL
);
349 isl_basic_set_free(bset
);
353 static __isl_give isl_mat
*extract_divs(__isl_keep isl_basic_set
*bset
)
356 isl_ctx
*ctx
= isl_basic_set_get_ctx(bset
);
359 div
= isl_mat_alloc(ctx
, bset
->n_div
,
360 1 + 1 + isl_basic_set_total_dim(bset
));
364 for (i
= 0; i
< bset
->n_div
; ++i
)
365 isl_seq_cpy(div
->row
[i
], bset
->div
[i
], div
->n_col
);
370 enum isl_lp_result
isl_basic_set_opt(__isl_keep isl_basic_set
*bset
, int max
,
371 __isl_keep isl_aff
*obj
, isl_int
*opt
)
376 isl_mat
*bset_div
= NULL
;
378 enum isl_lp_result res
;
379 int bset_n_div
, obj_n_div
;
384 ctx
= isl_aff_get_ctx(obj
);
385 if (!isl_space_is_equal(bset
->dim
, obj
->ls
->dim
))
386 isl_die(ctx
, isl_error_invalid
,
387 "spaces don't match", return isl_lp_error
);
388 if (!isl_int_is_one(obj
->v
->el
[0]))
389 isl_die(ctx
, isl_error_unsupported
,
390 "expecting integer affine expression",
391 return isl_lp_error
);
393 bset_n_div
= isl_basic_set_dim(bset
, isl_dim_div
);
394 obj_n_div
= isl_aff_dim(obj
, isl_dim_div
);
395 if (bset_n_div
== 0 && obj_n_div
== 0)
396 return basic_set_opt(bset
, max
, obj
, opt
);
398 bset
= isl_basic_set_copy(bset
);
399 obj
= isl_aff_copy(obj
);
401 bset_div
= extract_divs(bset
);
402 exp1
= isl_alloc_array(ctx
, int, bset_n_div
);
403 exp2
= isl_alloc_array(ctx
, int, obj_n_div
);
404 if (!bset_div
|| (bset_n_div
&& !exp1
) || (obj_n_div
&& !exp2
))
407 div
= isl_merge_divs(bset_div
, obj
->ls
->div
, exp1
, exp2
);
409 bset
= isl_basic_set_expand_divs(bset
, isl_mat_copy(div
), exp1
);
410 obj
= isl_aff_expand_divs(obj
, isl_mat_copy(div
), exp2
);
412 res
= basic_set_opt(bset
, max
, obj
, opt
);
414 isl_mat_free(bset_div
);
418 isl_basic_set_free(bset
);
424 isl_mat_free(bset_div
);
427 isl_basic_set_free(bset
);
432 /* Compute the minimum (maximum if max is set) of the integer affine
433 * expression obj over the points in set and put the result in *opt.
435 * The parameters are assumed to have been aligned.
437 static enum isl_lp_result
isl_set_opt_aligned(__isl_keep isl_set
*set
, int max
,
438 __isl_keep isl_aff
*obj
, isl_int
*opt
)
441 enum isl_lp_result res
;
450 res
= isl_basic_set_opt(set
->p
[0], max
, obj
, opt
);
451 if (res
== isl_lp_error
|| res
== isl_lp_unbounded
)
455 if (res
== isl_lp_ok
)
459 for (i
= 1; i
< set
->n
; ++i
) {
460 res
= isl_basic_set_opt(set
->p
[i
], max
, obj
, &opt_i
);
461 if (res
== isl_lp_error
|| res
== isl_lp_unbounded
) {
462 isl_int_clear(opt_i
);
465 if (res
== isl_lp_ok
)
467 if (max
? isl_int_gt(opt_i
, *opt
) : isl_int_lt(opt_i
, *opt
))
468 isl_int_set(*opt
, opt_i
);
470 isl_int_clear(opt_i
);
472 return empty
? isl_lp_empty
: isl_lp_ok
;
475 /* Compute the minimum (maximum if max is set) of the integer affine
476 * expression obj over the points in set and put the result in *opt.
478 enum isl_lp_result
isl_set_opt(__isl_keep isl_set
*set
, int max
,
479 __isl_keep isl_aff
*obj
, isl_int
*opt
)
481 enum isl_lp_result res
;
486 if (isl_space_match(set
->dim
, isl_dim_param
,
487 obj
->ls
->dim
, isl_dim_param
))
488 return isl_set_opt_aligned(set
, max
, obj
, opt
);
490 set
= isl_set_copy(set
);
491 obj
= isl_aff_copy(obj
);
492 set
= isl_set_align_params(set
, isl_aff_get_domain_space(obj
));
493 obj
= isl_aff_align_params(obj
, isl_set_get_space(set
));
495 res
= isl_set_opt_aligned(set
, max
, obj
, opt
);
503 enum isl_lp_result
isl_basic_set_max(__isl_keep isl_basic_set
*bset
,
504 __isl_keep isl_aff
*obj
, isl_int
*opt
)
506 return isl_basic_set_opt(bset
, 1, obj
, opt
);
509 enum isl_lp_result
isl_set_max(__isl_keep isl_set
*set
,
510 __isl_keep isl_aff
*obj
, isl_int
*opt
)
512 return isl_set_opt(set
, 1, obj
, opt
);
515 enum isl_lp_result
isl_set_min(__isl_keep isl_set
*set
,
516 __isl_keep isl_aff
*obj
, isl_int
*opt
)
518 return isl_set_opt(set
, 0, obj
, opt
);
521 /* Convert the result of a function that returns an isl_lp_result
522 * to an isl_val. The numerator of "v" is set to the optimal value
523 * if lp_res is isl_lp_ok. "max" is set if a maximum was computed.
525 * Return "v" with denominator set to 1 if lp_res is isl_lp_ok.
526 * Return NULL on error.
527 * Return a NaN if lp_res is isl_lp_empty.
528 * Return infinity or negative infinity if lp_res is isl_lp_unbounded,
529 * depending on "max".
531 static __isl_give isl_val
*convert_lp_result(enum isl_lp_result lp_res
,
532 __isl_take isl_val
*v
, int max
)
536 if (lp_res
== isl_lp_ok
) {
537 isl_int_set_si(v
->d
, 1);
538 return isl_val_normalize(v
);
540 ctx
= isl_val_get_ctx(v
);
542 if (lp_res
== isl_lp_error
)
544 if (lp_res
== isl_lp_empty
)
545 return isl_val_nan(ctx
);
547 return isl_val_infty(ctx
);
549 return isl_val_neginfty(ctx
);
552 /* Return the minimum (maximum if max is set) of the integer affine
553 * expression "obj" over the points in "bset".
555 * Return infinity or negative infinity if the optimal value is unbounded and
556 * NaN if "bset" is empty.
558 * Call isl_basic_set_opt and translate the results.
560 __isl_give isl_val
*isl_basic_set_opt_val(__isl_keep isl_basic_set
*bset
,
561 int max
, __isl_keep isl_aff
*obj
)
565 enum isl_lp_result lp_res
;
570 ctx
= isl_aff_get_ctx(obj
);
571 res
= isl_val_alloc(ctx
);
574 lp_res
= isl_basic_set_opt(bset
, max
, obj
, &res
->n
);
575 return convert_lp_result(lp_res
, res
, max
);
578 /* Return the maximum of the integer affine
579 * expression "obj" over the points in "bset".
581 * Return infinity or negative infinity if the optimal value is unbounded and
582 * NaN if "bset" is empty.
584 __isl_give isl_val
*isl_basic_set_max_val(__isl_keep isl_basic_set
*bset
,
585 __isl_keep isl_aff
*obj
)
587 return isl_basic_set_opt_val(bset
, 1, obj
);
590 /* Return the minimum (maximum if max is set) of the integer affine
591 * expression "obj" over the points in "set".
593 * Return infinity or negative infinity if the optimal value is unbounded and
594 * NaN if "bset" is empty.
596 * Call isl_set_opt and translate the results.
598 __isl_give isl_val
*isl_set_opt_val(__isl_keep isl_set
*set
, int max
,
599 __isl_keep isl_aff
*obj
)
603 enum isl_lp_result lp_res
;
608 ctx
= isl_aff_get_ctx(obj
);
609 res
= isl_val_alloc(ctx
);
612 lp_res
= isl_set_opt(set
, max
, obj
, &res
->n
);
613 return convert_lp_result(lp_res
, res
, max
);
616 /* Return the minimum of the integer affine
617 * expression "obj" over the points in "set".
619 * Return infinity or negative infinity if the optimal value is unbounded and
620 * NaN if "bset" is empty.
622 __isl_give isl_val
*isl_set_min_val(__isl_keep isl_set
*set
,
623 __isl_keep isl_aff
*obj
)
625 return isl_set_opt_val(set
, 0, obj
);
628 /* Return the maximum of the integer affine
629 * expression "obj" over the points in "set".
631 * Return infinity or negative infinity if the optimal value is unbounded and
632 * NaN if "bset" is empty.
634 __isl_give isl_val
*isl_set_max_val(__isl_keep isl_set
*set
,
635 __isl_keep isl_aff
*obj
)
637 return isl_set_opt_val(set
, 1, obj
);