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/union_set.h>
14 #include "isl_sample.h"
16 #include "isl_equalities.h"
17 #include <isl_aff_private.h>
18 #include <isl_local_space_private.h>
19 #include <isl_mat_private.h>
20 #include <isl_val_private.h>
21 #include <isl_vec_private.h>
22 #include <isl_lp_private.h>
23 #include <isl_ilp_private.h>
24 #include <isl/deprecated/ilp_int.h>
26 /* Given a basic set "bset", construct a basic set U such that for
27 * each element x in U, the whole unit box positioned at x is inside
28 * the given basic set.
29 * Note that U may not contain all points that satisfy this property.
31 * We simply add the sum of all negative coefficients to the constant
32 * term. This ensures that if x satisfies the resulting constraints,
33 * then x plus any sum of unit vectors satisfies the original constraints.
35 static __isl_give isl_basic_set
*unit_box_base_points(
36 __isl_take isl_basic_set
*bset
)
39 struct isl_basic_set
*unit_box
= NULL
;
45 if (bset
->n_eq
!= 0) {
46 isl_space
*space
= isl_basic_set_get_space(bset
);
47 isl_basic_set_free(bset
);
48 return isl_basic_set_empty(space
);
51 total
= isl_basic_set_total_dim(bset
);
52 unit_box
= isl_basic_set_alloc_space(isl_basic_set_get_space(bset
),
55 for (i
= 0; i
< bset
->n_ineq
; ++i
) {
56 k
= isl_basic_set_alloc_inequality(unit_box
);
59 isl_seq_cpy(unit_box
->ineq
[k
], bset
->ineq
[i
], 1 + total
);
60 for (j
= 0; j
< total
; ++j
) {
61 if (isl_int_is_nonneg(unit_box
->ineq
[k
][1 + j
]))
63 isl_int_add(unit_box
->ineq
[k
][0],
64 unit_box
->ineq
[k
][0], unit_box
->ineq
[k
][1 + j
]);
68 isl_basic_set_free(bset
);
71 isl_basic_set_free(bset
);
72 isl_basic_set_free(unit_box
);
76 /* Find an integer point in "bset", preferably one that is
77 * close to minimizing "f".
79 * We first check if we can easily put unit boxes inside bset.
80 * If so, we take the best base point of any of the unit boxes we can find
81 * and round it up to the nearest integer.
82 * If not, we simply pick any integer point in "bset".
84 static __isl_give isl_vec
*initial_solution(__isl_keep isl_basic_set
*bset
,
87 enum isl_lp_result res
;
88 struct isl_basic_set
*unit_box
;
91 unit_box
= unit_box_base_points(isl_basic_set_copy(bset
));
93 res
= isl_basic_set_solve_lp(unit_box
, 0, f
, bset
->ctx
->one
,
95 if (res
== isl_lp_ok
) {
96 isl_basic_set_free(unit_box
);
97 return isl_vec_ceil(sol
);
100 isl_basic_set_free(unit_box
);
102 return isl_basic_set_sample_vec(isl_basic_set_copy(bset
));
105 /* Restrict "bset" to those points with values for f in the interval [l, u].
107 static __isl_give isl_basic_set
*add_bounds(__isl_take isl_basic_set
*bset
,
108 isl_int
*f
, isl_int l
, isl_int u
)
113 total
= isl_basic_set_total_dim(bset
);
114 bset
= isl_basic_set_extend_constraints(bset
, 0, 2);
116 k
= isl_basic_set_alloc_inequality(bset
);
119 isl_seq_cpy(bset
->ineq
[k
], f
, 1 + total
);
120 isl_int_sub(bset
->ineq
[k
][0], bset
->ineq
[k
][0], l
);
122 k
= isl_basic_set_alloc_inequality(bset
);
125 isl_seq_neg(bset
->ineq
[k
], f
, 1 + total
);
126 isl_int_add(bset
->ineq
[k
][0], bset
->ineq
[k
][0], u
);
130 isl_basic_set_free(bset
);
134 /* Find an integer point in "bset" that minimizes f (in any) such that
135 * the value of f lies inside the interval [l, u].
136 * Return this integer point if it can be found.
137 * Otherwise, return sol.
139 * We perform a number of steps until l > u.
140 * In each step, we look for an integer point with value in either
141 * the whole interval [l, u] or half of the interval [l, l+floor(u-l-1/2)].
142 * The choice depends on whether we have found an integer point in the
143 * previous step. If so, we look for the next point in half of the remaining
145 * If we find a point, the current solution is updated and u is set
146 * to its value minus 1.
147 * If no point can be found, we update l to the upper bound of the interval
148 * we checked (u or l+floor(u-l-1/2)) plus 1.
150 static __isl_give isl_vec
*solve_ilp_search(__isl_keep isl_basic_set
*bset
,
151 isl_int
*f
, isl_int
*opt
, __isl_take isl_vec
*sol
, isl_int l
, isl_int u
)
158 while (isl_int_le(l
, u
)) {
159 struct isl_basic_set
*slice
;
160 struct isl_vec
*sample
;
165 isl_int_sub(tmp
, u
, l
);
166 isl_int_fdiv_q_ui(tmp
, tmp
, 2);
167 isl_int_add(tmp
, tmp
, l
);
169 slice
= add_bounds(isl_basic_set_copy(bset
), f
, l
, tmp
);
170 sample
= isl_basic_set_sample_vec(slice
);
176 if (sample
->size
> 0) {
179 isl_seq_inner_product(f
, sol
->el
, sol
->size
, opt
);
180 isl_int_sub_ui(u
, *opt
, 1);
183 isl_vec_free(sample
);
186 isl_int_add_ui(l
, tmp
, 1);
196 /* Find an integer point in "bset" that minimizes f (if any).
197 * If sol_p is not NULL then the integer point is returned in *sol_p.
198 * The optimal value of f is returned in *opt.
200 * The algorithm maintains a currently best solution and an interval [l, u]
201 * of values of f for which integer solutions could potentially still be found.
202 * The initial value of the best solution so far is any solution.
203 * The initial value of l is minimal value of f over the rationals
204 * (rounded up to the nearest integer).
205 * The initial value of u is the value of f at the initial solution minus 1.
207 * We then call solve_ilp_search to perform a binary search on the interval.
209 static enum isl_lp_result
solve_ilp(__isl_keep isl_basic_set
*bset
,
210 isl_int
*f
, isl_int
*opt
, __isl_give isl_vec
**sol_p
)
212 enum isl_lp_result res
;
216 res
= isl_basic_set_solve_lp(bset
, 0, f
, bset
->ctx
->one
,
218 if (res
== isl_lp_ok
&& isl_int_is_one(sol
->el
[0])) {
226 if (res
== isl_lp_error
|| res
== isl_lp_empty
)
229 sol
= initial_solution(bset
, f
);
232 if (sol
->size
== 0) {
236 if (res
== isl_lp_unbounded
) {
238 return isl_lp_unbounded
;
244 isl_int_set(l
, *opt
);
246 isl_seq_inner_product(f
, sol
->el
, sol
->size
, opt
);
247 isl_int_sub_ui(u
, *opt
, 1);
249 sol
= solve_ilp_search(bset
, f
, opt
, sol
, l
, u
);
264 static enum isl_lp_result
solve_ilp_with_eq(__isl_keep isl_basic_set
*bset
,
265 int max
, isl_int
*f
, isl_int
*opt
, __isl_give isl_vec
**sol_p
)
268 enum isl_lp_result res
;
269 struct isl_mat
*T
= NULL
;
272 bset
= isl_basic_set_copy(bset
);
273 dim
= isl_basic_set_total_dim(bset
);
274 v
= isl_vec_alloc(bset
->ctx
, 1 + dim
);
277 isl_seq_cpy(v
->el
, f
, 1 + dim
);
278 bset
= isl_basic_set_remove_equalities(bset
, &T
, NULL
);
279 v
= isl_vec_mat_product(v
, isl_mat_copy(T
));
282 res
= isl_basic_set_solve_ilp(bset
, max
, v
->el
, opt
, sol_p
);
284 if (res
== isl_lp_ok
&& sol_p
) {
285 *sol_p
= isl_mat_vec_product(T
, *sol_p
);
290 isl_basic_set_free(bset
);
294 isl_basic_set_free(bset
);
298 /* Find an integer point in "bset" that minimizes (or maximizes if max is set)
300 * If sol_p is not NULL then the integer point is returned in *sol_p.
301 * The optimal value of f is returned in *opt.
303 * If there is any equality among the points in "bset", then we first
304 * project it out. Otherwise, we continue with solve_ilp above.
306 enum isl_lp_result
isl_basic_set_solve_ilp(__isl_keep isl_basic_set
*bset
,
307 int max
, isl_int
*f
, isl_int
*opt
, __isl_give isl_vec
**sol_p
)
310 enum isl_lp_result res
;
317 isl_assert(bset
->ctx
, isl_basic_set_n_param(bset
) == 0,
318 return isl_lp_error
);
320 if (isl_basic_set_plain_is_empty(bset
))
324 return solve_ilp_with_eq(bset
, max
, f
, opt
, sol_p
);
326 dim
= isl_basic_set_total_dim(bset
);
329 isl_seq_neg(f
, f
, 1 + dim
);
331 res
= solve_ilp(bset
, f
, opt
, sol_p
);
334 isl_seq_neg(f
, f
, 1 + dim
);
335 isl_int_neg(*opt
, *opt
);
341 static enum isl_lp_result
basic_set_opt(__isl_keep isl_basic_set
*bset
, int max
,
342 __isl_keep isl_aff
*obj
, isl_int
*opt
)
344 enum isl_lp_result res
;
348 bset
= isl_basic_set_copy(bset
);
349 bset
= isl_basic_set_underlying_set(bset
);
350 res
= isl_basic_set_solve_ilp(bset
, max
, obj
->v
->el
+ 1, opt
, NULL
);
351 isl_basic_set_free(bset
);
355 static __isl_give isl_mat
*extract_divs(__isl_keep isl_basic_set
*bset
)
358 isl_ctx
*ctx
= isl_basic_set_get_ctx(bset
);
361 div
= isl_mat_alloc(ctx
, bset
->n_div
,
362 1 + 1 + isl_basic_set_total_dim(bset
));
366 for (i
= 0; i
< bset
->n_div
; ++i
)
367 isl_seq_cpy(div
->row
[i
], bset
->div
[i
], div
->n_col
);
372 enum isl_lp_result
isl_basic_set_opt(__isl_keep isl_basic_set
*bset
, int max
,
373 __isl_keep isl_aff
*obj
, isl_int
*opt
)
378 isl_mat
*bset_div
= NULL
;
380 enum isl_lp_result res
;
381 int bset_n_div
, obj_n_div
;
386 ctx
= isl_aff_get_ctx(obj
);
387 if (!isl_space_is_equal(bset
->dim
, obj
->ls
->dim
))
388 isl_die(ctx
, isl_error_invalid
,
389 "spaces don't match", return isl_lp_error
);
390 if (!isl_int_is_one(obj
->v
->el
[0]))
391 isl_die(ctx
, isl_error_unsupported
,
392 "expecting integer affine expression",
393 return isl_lp_error
);
395 bset_n_div
= isl_basic_set_dim(bset
, isl_dim_div
);
396 obj_n_div
= isl_aff_dim(obj
, isl_dim_div
);
397 if (bset_n_div
== 0 && obj_n_div
== 0)
398 return basic_set_opt(bset
, max
, obj
, opt
);
400 bset
= isl_basic_set_copy(bset
);
401 obj
= isl_aff_copy(obj
);
403 bset_div
= extract_divs(bset
);
404 exp1
= isl_alloc_array(ctx
, int, bset_n_div
);
405 exp2
= isl_alloc_array(ctx
, int, obj_n_div
);
406 if (!bset_div
|| (bset_n_div
&& !exp1
) || (obj_n_div
&& !exp2
))
409 div
= isl_merge_divs(bset_div
, obj
->ls
->div
, exp1
, exp2
);
411 bset
= isl_basic_set_expand_divs(bset
, isl_mat_copy(div
), exp1
);
412 obj
= isl_aff_expand_divs(obj
, isl_mat_copy(div
), exp2
);
414 res
= basic_set_opt(bset
, max
, obj
, opt
);
416 isl_mat_free(bset_div
);
420 isl_basic_set_free(bset
);
426 isl_mat_free(bset_div
);
429 isl_basic_set_free(bset
);
434 /* Compute the minimum (maximum if max is set) of the integer affine
435 * expression obj over the points in set and put the result in *opt.
437 * The parameters are assumed to have been aligned.
439 static enum isl_lp_result
isl_set_opt_aligned(__isl_keep isl_set
*set
, int max
,
440 __isl_keep isl_aff
*obj
, isl_int
*opt
)
443 enum isl_lp_result res
;
452 res
= isl_basic_set_opt(set
->p
[0], max
, obj
, opt
);
453 if (res
== isl_lp_error
|| res
== isl_lp_unbounded
)
457 if (res
== isl_lp_ok
)
461 for (i
= 1; i
< set
->n
; ++i
) {
462 res
= isl_basic_set_opt(set
->p
[i
], max
, obj
, &opt_i
);
463 if (res
== isl_lp_error
|| res
== isl_lp_unbounded
) {
464 isl_int_clear(opt_i
);
467 if (res
== isl_lp_empty
)
470 if (max
? isl_int_gt(opt_i
, *opt
) : isl_int_lt(opt_i
, *opt
))
471 isl_int_set(*opt
, opt_i
);
473 isl_int_clear(opt_i
);
475 return empty
? isl_lp_empty
: isl_lp_ok
;
478 /* Compute the minimum (maximum if max is set) of the integer affine
479 * expression obj over the points in set and put the result in *opt.
481 enum isl_lp_result
isl_set_opt(__isl_keep isl_set
*set
, int max
,
482 __isl_keep isl_aff
*obj
, isl_int
*opt
)
484 enum isl_lp_result res
;
490 aligned
= isl_set_space_has_equal_params(set
, obj
->ls
->dim
);
494 return isl_set_opt_aligned(set
, max
, obj
, opt
);
496 set
= isl_set_copy(set
);
497 obj
= isl_aff_copy(obj
);
498 set
= isl_set_align_params(set
, isl_aff_get_domain_space(obj
));
499 obj
= isl_aff_align_params(obj
, isl_set_get_space(set
));
501 res
= isl_set_opt_aligned(set
, max
, obj
, opt
);
509 enum isl_lp_result
isl_basic_set_max(__isl_keep isl_basic_set
*bset
,
510 __isl_keep isl_aff
*obj
, isl_int
*opt
)
512 return isl_basic_set_opt(bset
, 1, obj
, opt
);
515 enum isl_lp_result
isl_set_max(__isl_keep isl_set
*set
,
516 __isl_keep isl_aff
*obj
, isl_int
*opt
)
518 return isl_set_opt(set
, 1, obj
, opt
);
521 enum isl_lp_result
isl_set_min(__isl_keep isl_set
*set
,
522 __isl_keep isl_aff
*obj
, isl_int
*opt
)
524 return isl_set_opt(set
, 0, obj
, opt
);
527 /* Convert the result of a function that returns an isl_lp_result
528 * to an isl_val. The numerator of "v" is set to the optimal value
529 * if lp_res is isl_lp_ok. "max" is set if a maximum was computed.
531 * Return "v" with denominator set to 1 if lp_res is isl_lp_ok.
532 * Return NULL on error.
533 * Return a NaN if lp_res is isl_lp_empty.
534 * Return infinity or negative infinity if lp_res is isl_lp_unbounded,
535 * depending on "max".
537 static __isl_give isl_val
*convert_lp_result(enum isl_lp_result lp_res
,
538 __isl_take isl_val
*v
, int max
)
542 if (lp_res
== isl_lp_ok
) {
543 isl_int_set_si(v
->d
, 1);
544 return isl_val_normalize(v
);
546 ctx
= isl_val_get_ctx(v
);
548 if (lp_res
== isl_lp_error
)
550 if (lp_res
== isl_lp_empty
)
551 return isl_val_nan(ctx
);
553 return isl_val_infty(ctx
);
555 return isl_val_neginfty(ctx
);
558 /* Return the minimum (maximum if max is set) of the integer affine
559 * expression "obj" over the points in "bset".
561 * Return infinity or negative infinity if the optimal value is unbounded and
562 * NaN if "bset" is empty.
564 * Call isl_basic_set_opt and translate the results.
566 __isl_give isl_val
*isl_basic_set_opt_val(__isl_keep isl_basic_set
*bset
,
567 int max
, __isl_keep isl_aff
*obj
)
571 enum isl_lp_result lp_res
;
576 ctx
= isl_aff_get_ctx(obj
);
577 res
= isl_val_alloc(ctx
);
580 lp_res
= isl_basic_set_opt(bset
, max
, obj
, &res
->n
);
581 return convert_lp_result(lp_res
, res
, max
);
584 /* Return the maximum of the integer affine
585 * expression "obj" over the points in "bset".
587 * Return infinity or negative infinity if the optimal value is unbounded and
588 * NaN if "bset" is empty.
590 __isl_give isl_val
*isl_basic_set_max_val(__isl_keep isl_basic_set
*bset
,
591 __isl_keep isl_aff
*obj
)
593 return isl_basic_set_opt_val(bset
, 1, obj
);
596 /* Return the minimum (maximum if max is set) of the integer affine
597 * expression "obj" over the points in "set".
599 * Return infinity or negative infinity if the optimal value is unbounded and
600 * NaN if "set" is empty.
602 * Call isl_set_opt and translate the results.
604 __isl_give isl_val
*isl_set_opt_val(__isl_keep isl_set
*set
, int max
,
605 __isl_keep isl_aff
*obj
)
609 enum isl_lp_result lp_res
;
614 ctx
= isl_aff_get_ctx(obj
);
615 res
= isl_val_alloc(ctx
);
618 lp_res
= isl_set_opt(set
, max
, obj
, &res
->n
);
619 return convert_lp_result(lp_res
, res
, max
);
622 /* Return the minimum of the integer affine
623 * expression "obj" over the points in "set".
625 * Return infinity or negative infinity if the optimal value is unbounded and
626 * NaN if "set" is empty.
628 __isl_give isl_val
*isl_set_min_val(__isl_keep isl_set
*set
,
629 __isl_keep isl_aff
*obj
)
631 return isl_set_opt_val(set
, 0, obj
);
634 /* Return the maximum of the integer affine
635 * expression "obj" over the points in "set".
637 * Return infinity or negative infinity if the optimal value is unbounded and
638 * NaN if "set" is empty.
640 __isl_give isl_val
*isl_set_max_val(__isl_keep isl_set
*set
,
641 __isl_keep isl_aff
*obj
)
643 return isl_set_opt_val(set
, 1, obj
);
646 /* Return the optimum (min or max depending on "max") of "v1" and "v2",
647 * where either may be NaN, signifying an uninitialized value.
648 * That is, if either is NaN, then return the other one.
650 static __isl_give isl_val
*val_opt(__isl_take isl_val
*v1
,
651 __isl_take isl_val
*v2
, int max
)
655 if (isl_val_is_nan(v1
)) {
659 if (isl_val_is_nan(v2
)) {
664 return isl_val_max(v1
, v2
);
666 return isl_val_min(v1
, v2
);
673 /* Internal data structure for isl_set_opt_pw_aff.
675 * "max" is set if the maximum should be computed.
676 * "set" is the set over which the optimum should be computed.
677 * "res" contains the current optimum and is initialized to NaN.
679 struct isl_set_opt_data
{
686 /* Update the optimum in data->res with respect to the affine function
687 * "aff" defined over "set".
689 static isl_stat
piece_opt(__isl_take isl_set
*set
, __isl_take isl_aff
*aff
,
692 struct isl_set_opt_data
*data
= user
;
695 set
= isl_set_intersect(set
, isl_set_copy(data
->set
));
696 opt
= isl_set_opt_val(set
, data
->max
, aff
);
700 data
->res
= val_opt(data
->res
, opt
, data
->max
);
702 return isl_stat_error
;
707 /* Return the minimum (maximum if "max" is set) of the integer piecewise affine
708 * expression "obj" over the points in "set".
710 * Return infinity or negative infinity if the optimal value is unbounded and
711 * NaN if the intersection of "set" with the domain of "obj" is empty.
713 * Initialize the result to NaN and then update it for each of the pieces
716 static __isl_give isl_val
*isl_set_opt_pw_aff(__isl_keep isl_set
*set
, int max
,
717 __isl_keep isl_pw_aff
*obj
)
719 struct isl_set_opt_data data
= { max
, set
};
721 data
.res
= isl_val_nan(isl_set_get_ctx(set
));
722 if (isl_pw_aff_foreach_piece(obj
, &piece_opt
, &data
) < 0)
723 return isl_val_free(data
.res
);
728 /* Internal data structure for isl_union_set_opt_union_pw_aff.
730 * "max" is set if the maximum should be computed.
731 * "obj" is the objective function that needs to be optimized.
732 * "res" contains the current optimum and is initialized to NaN.
734 struct isl_union_set_opt_data
{
736 isl_union_pw_aff
*obj
;
741 /* Update the optimum in data->res with the optimum over "set".
742 * Do so by first extracting the matching objective function
745 static isl_stat
set_opt(__isl_take isl_set
*set
, void *user
)
747 struct isl_union_set_opt_data
*data
= user
;
752 space
= isl_set_get_space(set
);
753 space
= isl_space_from_domain(space
);
754 space
= isl_space_add_dims(space
, isl_dim_out
, 1);
755 pa
= isl_union_pw_aff_extract_pw_aff(data
->obj
, space
);
756 opt
= isl_set_opt_pw_aff(set
, data
->max
, pa
);
760 data
->res
= val_opt(data
->res
, opt
, data
->max
);
762 return isl_stat_error
;
767 /* Return the minimum (maximum if "max" is set) of the integer piecewise affine
768 * expression "obj" over the points in "uset".
770 * Return infinity or negative infinity if the optimal value is unbounded and
771 * NaN if the intersection of "uset" with the domain of "obj" is empty.
773 * Initialize the result to NaN and then update it for each of the sets
776 static __isl_give isl_val
*isl_union_set_opt_union_pw_aff(
777 __isl_keep isl_union_set
*uset
, int max
,
778 __isl_keep isl_union_pw_aff
*obj
)
780 struct isl_union_set_opt_data data
= { max
, obj
};
782 data
.res
= isl_val_nan(isl_union_set_get_ctx(uset
));
783 if (isl_union_set_foreach_set(uset
, &set_opt
, &data
) < 0)
784 return isl_val_free(data
.res
);
789 /* Return a list of minima (maxima if "max" is set) over the points in "uset"
790 * for each of the expressions in "obj".
792 * An element in the list is infinity or negative infinity if the optimal
793 * value of the corresponding expression is unbounded and
794 * NaN if the intersection of "uset" with the domain of the expression
797 * Iterate over all the expressions in "obj" and collect the results.
799 static __isl_give isl_multi_val
*isl_union_set_opt_multi_union_pw_aff(
800 __isl_keep isl_union_set
*uset
, int max
,
801 __isl_keep isl_multi_union_pw_aff
*obj
)
809 n
= isl_multi_union_pw_aff_dim(obj
, isl_dim_set
);
810 mv
= isl_multi_val_zero(isl_multi_union_pw_aff_get_space(obj
));
812 for (i
= 0; i
< n
; ++i
) {
814 isl_union_pw_aff
*upa
;
816 upa
= isl_multi_union_pw_aff_get_union_pw_aff(obj
, i
);
817 v
= isl_union_set_opt_union_pw_aff(uset
, max
, upa
);
818 isl_union_pw_aff_free(upa
);
819 mv
= isl_multi_val_set_val(mv
, i
, v
);
825 /* Return a list of minima over the points in "uset"
826 * for each of the expressions in "obj".
828 * An element in the list is infinity or negative infinity if the optimal
829 * value of the corresponding expression is unbounded and
830 * NaN if the intersection of "uset" with the domain of the expression
833 __isl_give isl_multi_val
*isl_union_set_min_multi_union_pw_aff(
834 __isl_keep isl_union_set
*uset
, __isl_keep isl_multi_union_pw_aff
*obj
)
836 return isl_union_set_opt_multi_union_pw_aff(uset
, 0, obj
);