| blob | eb644e5ac58c5e74b6d98c6e84a3b8fa5fafde7c |
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 */
10 #include <isl_ctx_private.h>
11 #include <isl_map_private.h>
13 #include <isl/vec.h>
14 #include <isl/mat.h>
15 #include <isl_seq.h>
19 #include <isl_factorization.h>
20 #include <isl_point_private.h>
21 #include <isl_options_private.h>
22 #include <isl_vec_private.h>
24 #include <bset_from_bmap.c>
27 {
33 }
35 /* Construct a zero sample of the same dimension as bset.
36 * As a special case, if bset is zero-dimensional, this
37 * function creates a zero-dimensional sample point.
38 */
40 {
49 }
52 }
55 {
57 isl_int t;
84 }
87 }
92 else
99 }
104 }
108 error:
112 }
114 /* Find a sample integer point, if any, in bset, which is known
115 * to have equalities. If bset contains no integer points, then
116 * return a zero-length vector.
117 * We simply remove the known equalities, compute a sample
118 * in the resulting bset, using the specified recurse function,
119 * and then transform the sample back to the original space.
120 */
123 {
134 else
137 }
139 /* Return a matrix containing the equalities of the tableau
140 * in constraint form. The tableau is assumed to have
141 * an associated bset that has been kept up-to-date.
142 */
144 {
172 else
176 }
179 error:
182 }
184 /* Compute and return an initial basis for the bounded tableau "tab".
185 *
186 * If the tableau is either full-dimensional or zero-dimensional,
187 * the we simply return an identity matrix.
188 * Otherwise, we construct a basis whose first directions correspond
189 * to equalities.
190 */
192 {
210 }
212 /* Compute the minimum of the current ("level") basis row over "tab"
213 * and store the result in position "level" of "min".
214 *
215 * This function assumes that at least one more row and at least
216 * one more element in the constraint array are available in the tableau.
217 */
220 {
223 }
225 /* Compute the maximum of the current ("level") basis row over "tab"
226 * and store the result in position "level" of "max".
227 *
228 * This function assumes that at least one more row and at least
229 * one more element in the constraint array are available in the tableau.
230 */
233 {
246 }
248 /* Perform a greedy search for an integer point in the set represented
249 * by "tab", given that the minimal rational value (rounded up to the
250 * nearest integer) at "level" is smaller than the maximal rational
251 * value (rounded down to the nearest integer).
252 *
253 * Return 1 if we have found an integer point (if tab->n_unbounded > 0
254 * then we may have only found integer values for the bounded dimensions
255 * and it is the responsibility of the caller to extend this solution
256 * to the unbounded dimensions).
257 * Return 0 if greedy search did not result in a solution.
258 * Return -1 if some error occurred.
259 *
260 * We assign a value half-way between the minimum and the maximum
261 * to the current dimension and check if the minimal value of the
262 * next dimension is still smaller than (or equal) to the maximal value.
263 * We continue this process until either
264 * - the minimal value (rounded up) is greater than the maximal value
265 * (rounded down). In this case, greedy search has failed.
266 * - we have exhausted all bounded dimensions, meaning that we have
267 * found a solution.
268 * - the sample value of the tableau is integral.
269 * - some error has occurred.
270 */
273 {
315 }
317 /* Given a tableau representing a set, find and return
318 * an integer point in the set, if there is any.
319 *
320 * We perform a depth first search
321 * for an integer point, by scanning all possible values in the range
322 * attained by a basis vector, where an initial basis may have been set
323 * by the calling function. Otherwise an initial basis that exploits
324 * the equalities in the tableau is created.
325 * tab->n_zero is currently ignored and is clobbered by this function.
326 *
327 * The tableau is allowed to have unbounded direction, but then
328 * the calling function needs to set an initial basis, with the
329 * unbounded directions last and with tab->n_unbounded set
330 * to the number of unbounded directions.
331 * Furthermore, the calling functions needs to add shifted copies
332 * of all constraints involving unbounded directions to ensure
333 * that any feasible rational value in these directions can be rounded
334 * up to yield a feasible integer value.
335 * In particular, let B define the given basis x' = B x
336 * and let T be the inverse of B, i.e., X = T x'.
337 * Let a x + c >= 0 be a constraint of the set represented by the tableau,
338 * or a T x' + c >= 0 in terms of the given basis. Assume that
339 * the bounded directions have an integer value, then we can safely
340 * round up the values for the unbounded directions if we make sure
341 * that x' not only satisfies the original constraint, but also
342 * the constraint "a T x' + c + s >= 0" with s the sum of all
343 * negative values in the last n_unbounded entries of "a T".
344 * The calling function therefore needs to add the constraint
345 * a x + c + s >= 0. The current function then scans the first
346 * directions for an integer value and once those have been found,
347 * it can compute "T ceil(B x)" to yield an integer point in the set.
348 * Note that during the search, the first rows of B may be changed
349 * by a basis reduction, but the last n_unbounded rows of B remain
350 * unaltered and are also not mixed into the first rows.
351 *
352 * The search is implemented iteratively. "level" identifies the current
353 * basis vector. "init" is true if we want the first value at the current
354 * level and false if we want the next value.
355 *
356 * At the start of each level, we first check if we can find a solution
357 * using greedy search. If not, we continue with the exhaustive search.
358 *
359 * The initial basis is the identity matrix. If the range in some direction
360 * contains more than one integer value, we perform basis reduction based
361 * on the value of ctx->opt->gbr
362 * - ISL_GBR_NEVER: never perform basis reduction
363 * - ISL_GBR_ONCE: only perform basis reduction the first
364 * time such a range is encountered
365 * - ISL_GBR_ALWAYS: always perform basis reduction when
366 * such a range is encountered
367 *
368 * When ctx->opt->gbr is set to ISL_GBR_ALWAYS, then we allow the basis
369 * reduction computation to return early. That is, as soon as it
370 * finds a reasonable first direction.
371 */
373 {
409 sample);
411 }
457 }
473 }
480 level--;
486 }
495 }
497 }
505 sample);
509 }
518 error:
524 }
528 /* Compute a sample point of the given basic set, based on the given,
529 * non-trivial factorization.
530 */
533 {
575 }
581 }
589 error:
594 }
596 /* Given a basic set that is known to be bounded, find and return
597 * an integer point in the basic set, if there is any.
598 *
599 * After handling some trivial cases, we construct a tableau
600 * and then use isl_tab_sample to find a sample, passing it
601 * the identity matrix as initial basis.
602 */
604 {
638 }
651 }
656 error:
660 }
662 /* Given a basic set "bset" and a value "sample" for the first coordinates
663 * of bset, plug in these values and drop the corresponding coordinates.
664 *
665 * We do this by computing the preimage of the transformation
666 *
667 * [ 1 0 ]
668 * x = [ s 0 ] x'
669 * [ 0 I ]
670 *
671 * where [1 s] is the sample value and I is the identity matrix of the
672 * appropriate dimension.
673 */
676 {
692 }
696 }
701 error:
705 }
707 /* Given a basic set "bset", return any (possibly non-integer) point
708 * in the basic set.
709 */
711 {
725 }
727 /* Given a linear cone "cone" and a rational point "vec",
728 * construct a polyhedron with shifted copies of the constraints in "cone",
729 * i.e., a polyhedron with "cone" as its recession cone, such that each
730 * point x in this polyhedron is such that the unit box positioned at x
731 * lies entirely inside the affine cone 'vec + cone'.
732 * Any rational point in this polyhedron may therefore be rounded up
733 * to yield an integer point that lies inside said affine cone.
734 *
735 * Denote the constraints of cone by "<a_i, x> >= 0" and the rational
736 * point "vec" by v/d.
737 * Let b_i = <a_i, v>. Then the affine cone 'vec + cone' is given
738 * by <a_i, x> - b/d >= 0.
739 * The polyhedron <a_i, x> - ceil{b/d} >= 0 is a subset of this affine cone.
740 * We prefer this polyhedron over the actual affine cone because it doesn't
741 * require a scaling of the constraints.
742 * If each of the vertices of the unit cube positioned at x lies inside
743 * this polyhedron, then the whole unit cube at x lies inside the affine cone.
744 * We therefore impose that x' = x + \sum e_i, for any selection of unit
745 * vectors lies inside the polyhedron, i.e.,
746 *
747 * <a_i, x'> - ceil{b/d} = <a_i, x> + sum a_i - ceil{b/d} >= 0
748 *
749 * The most stringent of these constraints is the one that selects
750 * all negative a_i, so the polyhedron we are looking for has constraints
751 *
752 * <a_i, x> + sum_{a_i < 0} a_i - ceil{b/d} >= 0
753 *
754 * Note that if cone were known to have only non-negative rays
755 * (which can be accomplished by a unimodular transformation),
756 * then we would only have to check the points x' = x + e_i
757 * and we only have to add the smallest negative a_i (if any)
758 * instead of the sum of all negative a_i.
759 */
762 {
793 }
794 }
800 error:
805 }
807 /* Given a rational point vec in a (transformed) basic set,
808 * such that cone is the recession cone of the original basic set,
809 * "round up" the rational point to an integer point.
810 *
811 * We first check if the rational point just happens to be integer.
812 * If not, we transform the cone in the same way as the basic set,
813 * pick a point x in this cone shifted to the rational point such that
814 * the whole unit cube at x is also inside this affine cone.
815 * Then we simply round up the coordinates of x and return the
816 * resulting integer point.
817 */
820 {
831 }
843 error:
848 }
850 /* Concatenate two integer vectors, i.e., two vectors with denominator
851 * (stored in element 0) equal to 1.
852 */
854 {
875 error:
879 }
881 /* Give a basic set "bset" with recession cone "cone", compute and
882 * return an integer point in bset, if any.
883 *
884 * If the recession cone is full-dimensional, then we know that
885 * bset contains an infinite number of integer points and it is
886 * fairly easy to pick one of them.
887 * If the recession cone is not full-dimensional, then we first
888 * transform bset such that the bounded directions appear as
889 * the first dimensions of the transformed basic set.
890 * We do this by using a unimodular transformation that transforms
891 * the equalities in the recession cone to equalities on the first
892 * dimensions.
893 *
894 * The transformed set is then projected onto its bounded dimensions.
895 * Note that to compute this projection, we can simply drop all constraints
896 * involving any of the unbounded dimensions since these constraints
897 * cannot be combined to produce a constraint on the bounded dimensions.
898 * To see this, assume that there is such a combination of constraints
899 * that produces a constraint on the bounded dimensions. This means
900 * that some combination of the unbounded dimensions has both an upper
901 * bound and a lower bound in terms of the bounded dimensions, but then
902 * this combination would be a bounded direction too and would have been
903 * transformed into a bounded dimensions.
904 *
905 * We then compute a sample value in the bounded dimensions.
906 * If no such value can be found, then the original set did not contain
907 * any integer points and we are done.
908 * Otherwise, we plug in the value we found in the bounded dimensions,
909 * project out these bounded dimensions and end up with a set with
910 * a full-dimensional recession cone.
911 * A sample point in this set is computed by "rounding up" any
912 * rational point in the set.
913 *
914 * The sample points in the bounded and unbounded dimensions are
915 * then combined into a single sample point and transformed back
916 * to the original space.
917 */
920 {
955 }
962 error:
966 }
969 {
977 }
979 /* Given a tableau "tab", a tableau "tab_cone" that corresponds
980 * to the recession cone and the inverse of a new basis U = inv(B),
981 * with the unbounded directions in B last,
982 * add constraints to "tab" that ensure any rational value
983 * in the unbounded directions can be rounded up to an integer value.
984 *
985 * The new basis is given by x' = B x, i.e., x = U x'.
986 * For any rational value of the last tab->n_unbounded coordinates
987 * in the update tableau, the value that is obtained by rounding
988 * up this value should be contained in the original tableau.
989 * For any constraint "a x + c >= 0", we therefore need to add
990 * a constraint "a x + c + s >= 0", with s the sum of all negative
991 * entries in the last elements of "a U".
992 *
993 * Since we are not interested in the first entries of any of the "a U",
994 * we first drop the columns of U that correpond to bounded directions.
995 */
998 {
1000 isl_int v;
1006 }
1035 }
1040 error:
1044 }
1046 /* Compute and return an initial basis for the possibly
1047 * unbounded tableau "tab". "tab_cone" is a tableau
1048 * for the corresponding recession cone.
1049 * Additionally, add constraints to "tab" that ensure
1050 * that any rational value for the unbounded directions
1051 * can be rounded up to an integer value.
1052 *
1053 * If the tableau is bounded, i.e., if the recession cone
1054 * is zero-dimensional, then we just use inital_basis.
1055 * Otherwise, we construct a basis whose first directions
1056 * correspond to equalities, followed by bounded directions,
1057 * i.e., equalities in the recession cone.
1058 * The remaining directions are then unbounded.
1059 */
1062 {
1073 }
1094 }
1096 /* Compute and return a sample point in bset using generalized basis
1097 * reduction. We first check if the input set has a non-trivial
1098 * recession cone. If so, we perform some extra preprocessing in
1099 * sample_with_cone. Otherwise, we directly perform generalized basis
1100 * reduction.
1101 */
1103 {
1118 error:
1121 }
1124 {
1146 }
1147 }
1160 error:
1163 }
1166 {
1168 }
1170 /* Compute an integer sample in "bset", where the caller guarantees
1171 * that "bset" is bounded.
1172 */
1174 {
1176 }
1179 {
1202 }
1206 error:
1210 }
1213 {
1224 }
1229 error:
1232 }
1235 {
1237 }
1240 {
1254 }
1259 error:
1262 }
1265 {
1267 }
1270 {
1279 }
1282 {
1296 }
1302 error:
1305 }