2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 * Copyright 2010 INRIA Saclay
4 * Copyright 2012 Ecole Normale Superieure
6 * Use of this software is governed by the MIT license
8 * Written by Sven Verdoolaege, K.U.Leuven, Departement
9 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
10 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
11 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
12 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
15 #include <isl_ctx_private.h>
16 #include <isl_map_private.h>
21 #include "isl_equalities.h"
22 #include "isl_sample.h"
24 #include <isl_mat_private.h>
25 #include <isl_vec_private.h>
27 #include <bset_to_bmap.c>
28 #include <bset_from_bmap.c>
29 #include <set_to_map.c>
30 #include <set_from_map.c>
32 __isl_give isl_basic_map
*isl_basic_map_implicit_equalities(
33 __isl_take isl_basic_map
*bmap
)
40 bmap
= isl_basic_map_gauss(bmap
, NULL
);
41 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
43 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
))
45 if (bmap
->n_ineq
<= 1)
48 tab
= isl_tab_from_basic_map(bmap
, 0);
49 if (isl_tab_detect_implicit_equalities(tab
) < 0)
51 bmap
= isl_basic_map_update_from_tab(bmap
, tab
);
53 bmap
= isl_basic_map_gauss(bmap
, NULL
);
54 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
);
58 isl_basic_map_free(bmap
);
62 struct isl_basic_set
*isl_basic_set_implicit_equalities(
63 struct isl_basic_set
*bset
)
65 return bset_from_bmap(
66 isl_basic_map_implicit_equalities(bset_to_bmap(bset
)));
69 /* Make eq[row][col] of both bmaps equal so we can add the row
70 * add the column to the common matrix.
71 * Note that because of the echelon form, the columns of row row
72 * after column col are zero.
74 static void set_common_multiple(
75 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
76 unsigned row
, unsigned col
)
80 if (isl_int_eq(bset1
->eq
[row
][col
], bset2
->eq
[row
][col
]))
85 isl_int_lcm(m
, bset1
->eq
[row
][col
], bset2
->eq
[row
][col
]);
86 isl_int_divexact(c
, m
, bset1
->eq
[row
][col
]);
87 isl_seq_scale(bset1
->eq
[row
], bset1
->eq
[row
], c
, col
+1);
88 isl_int_divexact(c
, m
, bset2
->eq
[row
][col
]);
89 isl_seq_scale(bset2
->eq
[row
], bset2
->eq
[row
], c
, col
+1);
94 /* Delete a given equality, moving all the following equalities one up.
96 static void delete_row(struct isl_basic_set
*bset
, unsigned row
)
103 for (r
= row
; r
< bset
->n_eq
; ++r
)
104 bset
->eq
[r
] = bset
->eq
[r
+1];
105 bset
->eq
[bset
->n_eq
] = t
;
108 /* Make first row entries in column col of bset1 identical to
109 * those of bset2, using the fact that entry bset1->eq[row][col]=a
110 * is non-zero. Initially, these elements of bset1 are all zero.
111 * For each row i < row, we set
112 * A[i] = a * A[i] + B[i][col] * A[row]
115 * A[i][col] = B[i][col] = a * old(B[i][col])
117 static isl_stat
construct_column(
118 __isl_keep isl_basic_set
*bset1
, __isl_keep isl_basic_set
*bset2
,
119 unsigned row
, unsigned col
)
126 total
= isl_basic_set_dim(bset1
, isl_dim_set
);
128 return isl_stat_error
;
132 for (r
= 0; r
< row
; ++r
) {
133 if (isl_int_is_zero(bset2
->eq
[r
][col
]))
135 isl_int_gcd(b
, bset2
->eq
[r
][col
], bset1
->eq
[row
][col
]);
136 isl_int_divexact(a
, bset1
->eq
[row
][col
], b
);
137 isl_int_divexact(b
, bset2
->eq
[r
][col
], b
);
138 isl_seq_combine(bset1
->eq
[r
], a
, bset1
->eq
[r
],
139 b
, bset1
->eq
[row
], 1 + total
);
140 isl_seq_scale(bset2
->eq
[r
], bset2
->eq
[r
], a
, 1 + total
);
144 delete_row(bset1
, row
);
149 /* Make first row entries in column col of bset1 identical to
150 * those of bset2, using only these entries of the two matrices.
151 * Let t be the last row with different entries.
152 * For each row i < t, we set
153 * A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
154 * B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
156 * A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
158 static isl_bool
transform_column(
159 __isl_keep isl_basic_set
*bset1
, __isl_keep isl_basic_set
*bset2
,
160 unsigned row
, unsigned col
)
166 for (t
= row
-1; t
>= 0; --t
)
167 if (isl_int_ne(bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]))
170 return isl_bool_false
;
172 total
= isl_basic_set_dim(bset1
, isl_dim_set
);
174 return isl_bool_error
;
178 isl_int_sub(b
, bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]);
179 for (i
= 0; i
< t
; ++i
) {
180 isl_int_sub(a
, bset2
->eq
[i
][col
], bset1
->eq
[i
][col
]);
181 isl_int_gcd(g
, a
, b
);
182 isl_int_divexact(a
, a
, g
);
183 isl_int_divexact(g
, b
, g
);
184 isl_seq_combine(bset1
->eq
[i
], g
, bset1
->eq
[i
], a
, bset1
->eq
[t
],
186 isl_seq_combine(bset2
->eq
[i
], g
, bset2
->eq
[i
], a
, bset2
->eq
[t
],
192 delete_row(bset1
, t
);
193 delete_row(bset2
, t
);
194 return isl_bool_true
;
197 /* The implementation is based on Section 5.2 of Michael Karr,
198 * "Affine Relationships Among Variables of a Program",
199 * except that the echelon form we use starts from the last column
200 * and that we are dealing with integer coefficients.
202 static __isl_give isl_basic_set
*affine_hull(
203 __isl_take isl_basic_set
*bset1
, __isl_take isl_basic_set
*bset2
)
210 dim
= isl_basic_set_dim(bset1
, isl_dim_set
);
211 if (dim
< 0 || !bset2
)
217 for (col
= total
-1; col
>= 0; --col
) {
218 int is_zero1
= row
>= bset1
->n_eq
||
219 isl_int_is_zero(bset1
->eq
[row
][col
]);
220 int is_zero2
= row
>= bset2
->n_eq
||
221 isl_int_is_zero(bset2
->eq
[row
][col
]);
222 if (!is_zero1
&& !is_zero2
) {
223 set_common_multiple(bset1
, bset2
, row
, col
);
225 } else if (!is_zero1
&& is_zero2
) {
226 if (construct_column(bset1
, bset2
, row
, col
) < 0)
228 } else if (is_zero1
&& !is_zero2
) {
229 if (construct_column(bset2
, bset1
, row
, col
) < 0)
234 transform
= transform_column(bset1
, bset2
, row
, col
);
241 isl_assert(bset1
->ctx
, row
== bset1
->n_eq
, goto error
);
242 isl_basic_set_free(bset2
);
243 bset1
= isl_basic_set_normalize_constraints(bset1
);
246 isl_basic_set_free(bset1
);
247 isl_basic_set_free(bset2
);
251 /* Find an integer point in the set represented by "tab"
252 * that lies outside of the equality "eq" e(x) = 0.
253 * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
254 * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
255 * The point, if found, is returned.
256 * If no point can be found, a zero-length vector is returned.
258 * Before solving an ILP problem, we first check if simply
259 * adding the normal of the constraint to one of the known
260 * integer points in the basic set represented by "tab"
261 * yields another point inside the basic set.
263 * The caller of this function ensures that the tableau is bounded or
264 * that tab->basis and tab->n_unbounded have been set appropriately.
266 static struct isl_vec
*outside_point(struct isl_tab
*tab
, isl_int
*eq
, int up
)
269 struct isl_vec
*sample
= NULL
;
270 struct isl_tab_undo
*snap
;
278 sample
= isl_vec_alloc(ctx
, 1 + dim
);
281 isl_int_set_si(sample
->el
[0], 1);
282 isl_seq_combine(sample
->el
+ 1,
283 ctx
->one
, tab
->bmap
->sample
->el
+ 1,
284 up
? ctx
->one
: ctx
->negone
, eq
+ 1, dim
);
285 if (isl_basic_map_contains(tab
->bmap
, sample
))
287 isl_vec_free(sample
);
290 snap
= isl_tab_snap(tab
);
293 isl_seq_neg(eq
, eq
, 1 + dim
);
294 isl_int_sub_ui(eq
[0], eq
[0], 1);
296 if (isl_tab_extend_cons(tab
, 1) < 0)
298 if (isl_tab_add_ineq(tab
, eq
) < 0)
301 sample
= isl_tab_sample(tab
);
303 isl_int_add_ui(eq
[0], eq
[0], 1);
305 isl_seq_neg(eq
, eq
, 1 + dim
);
307 if (sample
&& isl_tab_rollback(tab
, snap
) < 0)
312 isl_vec_free(sample
);
316 __isl_give isl_basic_set
*isl_basic_set_recession_cone(
317 __isl_take isl_basic_set
*bset
)
322 empty
= isl_basic_set_plain_is_empty(bset
);
324 return isl_basic_set_free(bset
);
328 bset
= isl_basic_set_cow(bset
);
329 if (isl_basic_set_check_no_locals(bset
) < 0)
330 return isl_basic_set_free(bset
);
332 for (i
= 0; i
< bset
->n_eq
; ++i
)
333 isl_int_set_si(bset
->eq
[i
][0], 0);
335 for (i
= 0; i
< bset
->n_ineq
; ++i
)
336 isl_int_set_si(bset
->ineq
[i
][0], 0);
338 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
339 return isl_basic_set_implicit_equalities(bset
);
342 /* Move "sample" to a point that is one up (or down) from the original
343 * point in dimension "pos".
345 static void adjacent_point(__isl_keep isl_vec
*sample
, int pos
, int up
)
348 isl_int_add_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
350 isl_int_sub_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
353 /* Check if any points that are adjacent to "sample" also belong to "bset".
354 * If so, add them to "hull" and return the updated hull.
356 * Before checking whether and adjacent point belongs to "bset", we first
357 * check whether it already belongs to "hull" as this test is typically
360 static __isl_give isl_basic_set
*add_adjacent_points(
361 __isl_take isl_basic_set
*hull
, __isl_take isl_vec
*sample
,
362 __isl_keep isl_basic_set
*bset
)
367 dim
= isl_basic_set_dim(hull
, isl_dim_set
);
368 if (!sample
|| dim
< 0)
371 for (i
= 0; i
< dim
; ++i
) {
372 for (up
= 0; up
<= 1; ++up
) {
374 isl_basic_set
*point
;
376 adjacent_point(sample
, i
, up
);
377 contains
= isl_basic_set_contains(hull
, sample
);
381 adjacent_point(sample
, i
, !up
);
384 contains
= isl_basic_set_contains(bset
, sample
);
388 point
= isl_basic_set_from_vec(
389 isl_vec_copy(sample
));
390 hull
= affine_hull(hull
, point
);
392 adjacent_point(sample
, i
, !up
);
398 isl_vec_free(sample
);
402 isl_vec_free(sample
);
403 isl_basic_set_free(hull
);
407 /* Extend an initial (under-)approximation of the affine hull of basic
408 * set represented by the tableau "tab"
409 * by looking for points that do not satisfy one of the equalities
410 * in the current approximation and adding them to that approximation
411 * until no such points can be found any more.
413 * The caller of this function ensures that "tab" is bounded or
414 * that tab->basis and tab->n_unbounded have been set appropriately.
416 * "bset" may be either NULL or the basic set represented by "tab".
417 * If "bset" is not NULL, we check for any point we find if any
418 * of its adjacent points also belong to "bset".
420 static __isl_give isl_basic_set
*extend_affine_hull(struct isl_tab
*tab
,
421 __isl_take isl_basic_set
*hull
, __isl_keep isl_basic_set
*bset
)
431 if (isl_tab_extend_cons(tab
, 2 * dim
+ 1) < 0)
434 for (i
= 0; i
< dim
; ++i
) {
435 struct isl_vec
*sample
;
436 struct isl_basic_set
*point
;
437 for (j
= 0; j
< hull
->n_eq
; ++j
) {
438 sample
= outside_point(tab
, hull
->eq
[j
], 1);
441 if (sample
->size
> 0)
443 isl_vec_free(sample
);
444 sample
= outside_point(tab
, hull
->eq
[j
], 0);
447 if (sample
->size
> 0)
449 isl_vec_free(sample
);
451 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
457 isl_tab_add_sample(tab
, isl_vec_copy(sample
)) < 0)
458 hull
= isl_basic_set_free(hull
);
460 hull
= add_adjacent_points(hull
, isl_vec_copy(sample
),
462 point
= isl_basic_set_from_vec(sample
);
463 hull
= affine_hull(hull
, point
);
470 isl_basic_set_free(hull
);
474 /* Construct an initial underapproximation of the hull of "bset"
475 * from "sample" and any of its adjacent points that also belong to "bset".
477 static __isl_give isl_basic_set
*initialize_hull(__isl_keep isl_basic_set
*bset
,
478 __isl_take isl_vec
*sample
)
482 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
483 hull
= add_adjacent_points(hull
, sample
, bset
);
488 /* Look for all equalities satisfied by the integer points in bset,
489 * which is assumed to be bounded.
491 * The equalities are obtained by successively looking for
492 * a point that is affinely independent of the points found so far.
493 * In particular, for each equality satisfied by the points so far,
494 * we check if there is any point on a hyperplane parallel to the
495 * corresponding hyperplane shifted by at least one (in either direction).
497 static __isl_give isl_basic_set
*uset_affine_hull_bounded(
498 __isl_take isl_basic_set
*bset
)
500 struct isl_vec
*sample
= NULL
;
501 struct isl_basic_set
*hull
;
502 struct isl_tab
*tab
= NULL
;
505 if (isl_basic_set_plain_is_empty(bset
))
508 dim
= isl_basic_set_dim(bset
, isl_dim_set
);
510 return isl_basic_set_free(bset
);
512 if (bset
->sample
&& bset
->sample
->size
== 1 + dim
) {
513 int contains
= isl_basic_set_contains(bset
, bset
->sample
);
519 sample
= isl_vec_copy(bset
->sample
);
521 isl_vec_free(bset
->sample
);
526 tab
= isl_tab_from_basic_set(bset
, 1);
531 isl_vec_free(sample
);
532 return isl_basic_set_set_to_empty(bset
);
536 struct isl_tab_undo
*snap
;
537 snap
= isl_tab_snap(tab
);
538 sample
= isl_tab_sample(tab
);
539 if (isl_tab_rollback(tab
, snap
) < 0)
541 isl_vec_free(tab
->bmap
->sample
);
542 tab
->bmap
->sample
= isl_vec_copy(sample
);
547 if (sample
->size
== 0) {
549 isl_vec_free(sample
);
550 return isl_basic_set_set_to_empty(bset
);
553 hull
= initialize_hull(bset
, sample
);
555 hull
= extend_affine_hull(tab
, hull
, bset
);
556 isl_basic_set_free(bset
);
561 isl_vec_free(sample
);
563 isl_basic_set_free(bset
);
567 /* Given an unbounded tableau and an integer point satisfying the tableau,
568 * construct an initial affine hull containing the recession cone
569 * shifted to the given point.
571 * The unbounded directions are taken from the last rows of the basis,
572 * which is assumed to have been initialized appropriately.
574 static __isl_give isl_basic_set
*initial_hull(struct isl_tab
*tab
,
575 __isl_take isl_vec
*vec
)
579 struct isl_basic_set
*bset
= NULL
;
586 isl_assert(ctx
, vec
->size
!= 0, goto error
);
588 bset
= isl_basic_set_alloc(ctx
, 0, vec
->size
- 1, 0, vec
->size
- 1, 0);
589 dim
= isl_basic_set_dim(bset
, isl_dim_set
);
592 dim
-= tab
->n_unbounded
;
593 for (i
= 0; i
< dim
; ++i
) {
594 k
= isl_basic_set_alloc_equality(bset
);
597 isl_seq_cpy(bset
->eq
[k
] + 1, tab
->basis
->row
[1 + i
] + 1,
599 isl_seq_inner_product(bset
->eq
[k
] + 1, vec
->el
+1,
600 vec
->size
- 1, &bset
->eq
[k
][0]);
601 isl_int_neg(bset
->eq
[k
][0], bset
->eq
[k
][0]);
604 bset
= isl_basic_set_gauss(bset
, NULL
);
608 isl_basic_set_free(bset
);
613 /* Given a tableau of a set and a tableau of the corresponding
614 * recession cone, detect and add all equalities to the tableau.
615 * If the tableau is bounded, then we can simply keep the
616 * tableau in its state after the return from extend_affine_hull.
617 * However, if the tableau is unbounded, then
618 * isl_tab_set_initial_basis_with_cone will add some additional
619 * constraints to the tableau that have to be removed again.
620 * In this case, we therefore rollback to the state before
621 * any constraints were added and then add the equalities back in.
623 struct isl_tab
*isl_tab_detect_equalities(struct isl_tab
*tab
,
624 struct isl_tab
*tab_cone
)
627 struct isl_vec
*sample
;
628 struct isl_basic_set
*hull
= NULL
;
629 struct isl_tab_undo
*snap
;
631 if (!tab
|| !tab_cone
)
634 snap
= isl_tab_snap(tab
);
636 isl_mat_free(tab
->basis
);
639 isl_assert(tab
->mat
->ctx
, tab
->bmap
, goto error
);
640 isl_assert(tab
->mat
->ctx
, tab
->samples
, goto error
);
641 isl_assert(tab
->mat
->ctx
, tab
->samples
->n_col
== 1 + tab
->n_var
, goto error
);
642 isl_assert(tab
->mat
->ctx
, tab
->n_sample
> tab
->n_outside
, goto error
);
644 if (isl_tab_set_initial_basis_with_cone(tab
, tab_cone
) < 0)
647 sample
= isl_vec_alloc(tab
->mat
->ctx
, 1 + tab
->n_var
);
651 isl_seq_cpy(sample
->el
, tab
->samples
->row
[tab
->n_outside
], sample
->size
);
653 isl_vec_free(tab
->bmap
->sample
);
654 tab
->bmap
->sample
= isl_vec_copy(sample
);
656 if (tab
->n_unbounded
== 0)
657 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
659 hull
= initial_hull(tab
, isl_vec_copy(sample
));
661 for (j
= tab
->n_outside
+ 1; j
< tab
->n_sample
; ++j
) {
662 isl_seq_cpy(sample
->el
, tab
->samples
->row
[j
], sample
->size
);
663 hull
= affine_hull(hull
,
664 isl_basic_set_from_vec(isl_vec_copy(sample
)));
667 isl_vec_free(sample
);
669 hull
= extend_affine_hull(tab
, hull
, NULL
);
673 if (tab
->n_unbounded
== 0) {
674 isl_basic_set_free(hull
);
678 if (isl_tab_rollback(tab
, snap
) < 0)
681 if (hull
->n_eq
> tab
->n_zero
) {
682 for (j
= 0; j
< hull
->n_eq
; ++j
) {
683 isl_seq_normalize(tab
->mat
->ctx
, hull
->eq
[j
], 1 + tab
->n_var
);
684 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
689 isl_basic_set_free(hull
);
693 isl_basic_set_free(hull
);
698 /* Compute the affine hull of "bset", where "cone" is the recession cone
701 * We first compute a unimodular transformation that puts the unbounded
702 * directions in the last dimensions. In particular, we take a transformation
703 * that maps all equalities to equalities (in HNF) on the first dimensions.
704 * Let x be the original dimensions and y the transformed, with y_1 bounded
707 * [ y_1 ] [ y_1 ] [ Q_1 ]
708 * x = U [ y_2 ] [ y_2 ] = [ Q_2 ] x
710 * Let's call the input basic set S. We compute S' = preimage(S, U)
711 * and drop the final dimensions including any constraints involving them.
712 * This results in set S''.
713 * Then we compute the affine hull A'' of S''.
714 * Let F y_1 >= g be the constraint system of A''. In the transformed
715 * space the y_2 are unbounded, so we can add them back without any constraints,
719 * [ F 0 ] [ y_2 ] >= g
722 * [ F 0 ] [ Q_2 ] x >= g
726 * The affine hull in the original space is then obtained as
727 * A = preimage(A'', Q_1).
729 static struct isl_basic_set
*affine_hull_with_cone(struct isl_basic_set
*bset
,
730 struct isl_basic_set
*cone
)
734 struct isl_basic_set
*hull
;
735 struct isl_mat
*M
, *U
, *Q
;
737 total
= isl_basic_set_dim(cone
, isl_dim_all
);
738 if (!bset
|| total
< 0)
741 cone_dim
= total
- cone
->n_eq
;
743 M
= isl_mat_sub_alloc6(bset
->ctx
, cone
->eq
, 0, cone
->n_eq
, 1, total
);
744 M
= isl_mat_left_hermite(M
, 0, &U
, &Q
);
749 U
= isl_mat_lin_to_aff(U
);
750 bset
= isl_basic_set_preimage(bset
, isl_mat_copy(U
));
752 bset
= isl_basic_set_drop_constraints_involving(bset
, total
- cone_dim
,
754 bset
= isl_basic_set_drop_dims(bset
, total
- cone_dim
, cone_dim
);
756 Q
= isl_mat_lin_to_aff(Q
);
757 Q
= isl_mat_drop_rows(Q
, 1 + total
- cone_dim
, cone_dim
);
759 if (bset
&& bset
->sample
&& bset
->sample
->size
== 1 + total
)
760 bset
->sample
= isl_mat_vec_product(isl_mat_copy(Q
), bset
->sample
);
762 hull
= uset_affine_hull_bounded(bset
);
768 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
769 U
= isl_mat_drop_cols(U
, 1 + total
- cone_dim
, cone_dim
);
770 if (sample
&& sample
->size
> 0)
771 sample
= isl_mat_vec_product(U
, sample
);
774 hull
= isl_basic_set_preimage(hull
, Q
);
776 isl_vec_free(hull
->sample
);
777 hull
->sample
= sample
;
779 isl_vec_free(sample
);
782 isl_basic_set_free(cone
);
786 isl_basic_set_free(bset
);
787 isl_basic_set_free(cone
);
791 /* Look for all equalities satisfied by the integer points in bset,
792 * which is assumed not to have any explicit equalities.
794 * The equalities are obtained by successively looking for
795 * a point that is affinely independent of the points found so far.
796 * In particular, for each equality satisfied by the points so far,
797 * we check if there is any point on a hyperplane parallel to the
798 * corresponding hyperplane shifted by at least one (in either direction).
800 * Before looking for any outside points, we first compute the recession
801 * cone. The directions of this recession cone will always be part
802 * of the affine hull, so there is no need for looking for any points
803 * in these directions.
804 * In particular, if the recession cone is full-dimensional, then
805 * the affine hull is simply the whole universe.
807 static struct isl_basic_set
*uset_affine_hull(struct isl_basic_set
*bset
)
809 struct isl_basic_set
*cone
;
812 if (isl_basic_set_plain_is_empty(bset
))
815 cone
= isl_basic_set_recession_cone(isl_basic_set_copy(bset
));
818 if (cone
->n_eq
== 0) {
820 space
= isl_basic_set_get_space(bset
);
821 isl_basic_set_free(cone
);
822 isl_basic_set_free(bset
);
823 return isl_basic_set_universe(space
);
826 total
= isl_basic_set_dim(cone
, isl_dim_all
);
828 bset
= isl_basic_set_free(bset
);
829 if (cone
->n_eq
< total
)
830 return affine_hull_with_cone(bset
, cone
);
832 isl_basic_set_free(cone
);
833 return uset_affine_hull_bounded(bset
);
835 isl_basic_set_free(bset
);
839 /* Look for all equalities satisfied by the integer points in bmap
840 * that are independent of the equalities already explicitly available
843 * We first remove all equalities already explicitly available,
844 * then look for additional equalities in the reduced space
845 * and then transform the result to the original space.
846 * The original equalities are _not_ added to this set. This is
847 * the responsibility of the calling function.
848 * The resulting basic set has all meaning about the dimensions removed.
849 * In particular, dimensions that correspond to existential variables
850 * in bmap and that are found to be fixed are not removed.
852 static struct isl_basic_set
*equalities_in_underlying_set(
853 struct isl_basic_map
*bmap
)
855 struct isl_mat
*T1
= NULL
;
856 struct isl_mat
*T2
= NULL
;
857 struct isl_basic_set
*bset
= NULL
;
858 struct isl_basic_set
*hull
= NULL
;
860 bset
= isl_basic_map_underlying_set(bmap
);
864 bset
= isl_basic_set_remove_equalities(bset
, &T1
, &T2
);
868 hull
= uset_affine_hull(bset
);
876 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
877 if (sample
&& sample
->size
> 0)
878 sample
= isl_mat_vec_product(T1
, sample
);
881 hull
= isl_basic_set_preimage(hull
, T2
);
883 isl_vec_free(hull
->sample
);
884 hull
->sample
= sample
;
886 isl_vec_free(sample
);
893 isl_basic_set_free(bset
);
894 isl_basic_set_free(hull
);
898 /* Detect and make explicit all equalities satisfied by the (integer)
901 __isl_give isl_basic_map
*isl_basic_map_detect_equalities(
902 __isl_take isl_basic_map
*bmap
)
906 struct isl_basic_set
*hull
= NULL
;
910 if (bmap
->n_ineq
== 0)
912 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
914 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_ALL_EQUALITIES
))
916 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_RATIONAL
))
917 return isl_basic_map_implicit_equalities(bmap
);
919 hull
= equalities_in_underlying_set(isl_basic_map_copy(bmap
));
922 if (ISL_F_ISSET(hull
, ISL_BASIC_SET_EMPTY
)) {
923 isl_basic_set_free(hull
);
924 return isl_basic_map_set_to_empty(bmap
);
926 bmap
= isl_basic_map_extend(bmap
, 0, hull
->n_eq
, 0);
927 total
= isl_basic_set_dim(hull
, isl_dim_all
);
930 for (i
= 0; i
< hull
->n_eq
; ++i
) {
931 j
= isl_basic_map_alloc_equality(bmap
);
934 isl_seq_cpy(bmap
->eq
[j
], hull
->eq
[i
], 1 + total
);
936 isl_vec_free(bmap
->sample
);
937 bmap
->sample
= isl_vec_copy(hull
->sample
);
938 isl_basic_set_free(hull
);
939 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
| ISL_BASIC_MAP_ALL_EQUALITIES
);
940 bmap
= isl_basic_map_simplify(bmap
);
941 return isl_basic_map_finalize(bmap
);
943 isl_basic_set_free(hull
);
944 isl_basic_map_free(bmap
);
948 __isl_give isl_basic_set
*isl_basic_set_detect_equalities(
949 __isl_take isl_basic_set
*bset
)
951 return bset_from_bmap(
952 isl_basic_map_detect_equalities(bset_to_bmap(bset
)));
955 __isl_give isl_map
*isl_map_detect_equalities(__isl_take isl_map
*map
)
957 return isl_map_inline_foreach_basic_map(map
,
958 &isl_basic_map_detect_equalities
);
961 __isl_give isl_set
*isl_set_detect_equalities(__isl_take isl_set
*set
)
963 return set_from_map(isl_map_detect_equalities(set_to_map(set
)));
966 /* Return the superset of "bmap" described by the equalities
967 * satisfied by "bmap" that are already known.
969 __isl_give isl_basic_map
*isl_basic_map_plain_affine_hull(
970 __isl_take isl_basic_map
*bmap
)
972 bmap
= isl_basic_map_cow(bmap
);
974 isl_basic_map_free_inequality(bmap
, bmap
->n_ineq
);
975 bmap
= isl_basic_map_finalize(bmap
);
979 /* Return the superset of "bset" described by the equalities
980 * satisfied by "bset" that are already known.
982 __isl_give isl_basic_set
*isl_basic_set_plain_affine_hull(
983 __isl_take isl_basic_set
*bset
)
985 return isl_basic_map_plain_affine_hull(bset
);
988 /* After computing the rational affine hull (by detecting the implicit
989 * equalities), we compute the additional equalities satisfied by
990 * the integer points (if any) and add the original equalities back in.
992 __isl_give isl_basic_map
*isl_basic_map_affine_hull(
993 __isl_take isl_basic_map
*bmap
)
995 bmap
= isl_basic_map_detect_equalities(bmap
);
996 bmap
= isl_basic_map_plain_affine_hull(bmap
);
1000 struct isl_basic_set
*isl_basic_set_affine_hull(struct isl_basic_set
*bset
)
1002 return bset_from_bmap(isl_basic_map_affine_hull(bset_to_bmap(bset
)));
1005 /* Given a rational affine matrix "M", add stride constraints to "bmap"
1010 * is an integer vector. The variables x include all the variables
1011 * of "bmap" except the unknown divs.
1013 * If d is the common denominator of M, then we need to impose that
1019 * exists alpha : d M(x) = d alpha
1021 * This function is similar to add_strides in isl_morph.c
1023 static __isl_give isl_basic_map
*add_strides(__isl_take isl_basic_map
*bmap
,
1024 __isl_keep isl_mat
*M
, int n_known
)
1029 if (isl_int_is_one(M
->row
[0][0]))
1032 bmap
= isl_basic_map_extend(bmap
, M
->n_row
- 1, M
->n_row
- 1, 0);
1035 for (i
= 1; i
< M
->n_row
; ++i
) {
1036 isl_seq_gcd(M
->row
[i
], M
->n_col
, &gcd
);
1037 if (isl_int_is_divisible_by(gcd
, M
->row
[0][0]))
1039 div
= isl_basic_map_alloc_div(bmap
);
1042 isl_int_set_si(bmap
->div
[div
][0], 0);
1043 k
= isl_basic_map_alloc_equality(bmap
);
1046 isl_seq_cpy(bmap
->eq
[k
], M
->row
[i
], M
->n_col
);
1047 isl_seq_clr(bmap
->eq
[k
] + M
->n_col
, bmap
->n_div
- n_known
);
1048 isl_int_set(bmap
->eq
[k
][M
->n_col
- n_known
+ div
],
1056 isl_basic_map_free(bmap
);
1060 /* If there are any equalities that involve (multiple) unknown divs,
1061 * then extract the stride information encoded by those equalities
1062 * and make it explicitly available in "bmap".
1064 * We first sort the divs so that the unknown divs appear last and
1065 * then we count how many equalities involve these divs.
1067 * Let these equalities be of the form
1071 * where y represents the unknown divs and x the remaining variables.
1072 * Let [H 0] be the Hermite Normal Form of B, i.e.,
1076 * Then x is a solution of the equalities iff
1078 * H^-1 A(x) (= - [I 0] Q y)
1080 * is an integer vector. Let d be the common denominator of H^-1.
1083 * d H^-1 A(x) = d alpha
1085 * in add_strides, with alpha fresh existentially quantified variables.
1087 static __isl_give isl_basic_map
*isl_basic_map_make_strides_explicit(
1088 __isl_take isl_basic_map
*bmap
)
1097 known
= isl_basic_map_divs_known(bmap
);
1099 return isl_basic_map_free(bmap
);
1102 bmap
= isl_basic_map_sort_divs(bmap
);
1103 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1107 for (n_known
= 0; n_known
< bmap
->n_div
; ++n_known
)
1108 if (isl_int_is_zero(bmap
->div
[n_known
][0]))
1110 ctx
= isl_basic_map_get_ctx(bmap
);
1111 v_div
= isl_basic_map_var_offset(bmap
, isl_dim_div
);
1113 return isl_basic_map_free(bmap
);
1114 for (n
= 0; n
< bmap
->n_eq
; ++n
)
1115 if (isl_seq_first_non_zero(bmap
->eq
[n
] + 1 + v_div
+ n_known
,
1116 bmap
->n_div
- n_known
) == -1)
1120 B
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 0, 1 + v_div
+ n_known
);
1121 n_col
= bmap
->n_div
- n_known
;
1122 A
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 1 + v_div
+ n_known
, n_col
);
1123 A
= isl_mat_left_hermite(A
, 0, NULL
, NULL
);
1124 A
= isl_mat_drop_cols(A
, n
, n_col
- n
);
1125 A
= isl_mat_lin_to_aff(A
);
1126 A
= isl_mat_right_inverse(A
);
1127 B
= isl_mat_insert_zero_rows(B
, 0, 1);
1128 B
= isl_mat_set_element_si(B
, 0, 0, 1);
1129 M
= isl_mat_product(A
, B
);
1131 return isl_basic_map_free(bmap
);
1132 bmap
= add_strides(bmap
, M
, n_known
);
1133 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1139 /* Compute the affine hull of each basic map in "map" separately
1140 * and make all stride information explicit so that we can remove
1141 * all unknown divs without losing this information.
1142 * The result is also guaranteed to be gaussed.
1144 * In simple cases where a div is determined by an equality,
1145 * calling isl_basic_map_gauss is enough to make the stride information
1146 * explicit, as it will derive an explicit representation for the div
1147 * from the equality. If, however, the stride information
1148 * is encoded through multiple unknown divs then we need to make
1149 * some extra effort in isl_basic_map_make_strides_explicit.
1151 static __isl_give isl_map
*isl_map_local_affine_hull(__isl_take isl_map
*map
)
1155 map
= isl_map_cow(map
);
1159 for (i
= 0; i
< map
->n
; ++i
) {
1160 map
->p
[i
] = isl_basic_map_affine_hull(map
->p
[i
]);
1161 map
->p
[i
] = isl_basic_map_gauss(map
->p
[i
], NULL
);
1162 map
->p
[i
] = isl_basic_map_make_strides_explicit(map
->p
[i
]);
1164 return isl_map_free(map
);
1170 static __isl_give isl_set
*isl_set_local_affine_hull(__isl_take isl_set
*set
)
1172 return isl_map_local_affine_hull(set
);
1175 /* Return an empty basic map living in the same space as "map".
1177 static __isl_give isl_basic_map
*replace_map_by_empty_basic_map(
1178 __isl_take isl_map
*map
)
1182 space
= isl_map_get_space(map
);
1184 return isl_basic_map_empty(space
);
1187 /* Compute the affine hull of "map".
1189 * We first compute the affine hull of each basic map separately.
1190 * Then we align the divs and recompute the affine hulls of the basic
1191 * maps since some of them may now have extra divs.
1192 * In order to avoid performing parametric integer programming to
1193 * compute explicit expressions for the divs, possible leading to
1194 * an explosion in the number of basic maps, we first drop all unknown
1195 * divs before aligning the divs. Note that isl_map_local_affine_hull tries
1196 * to make sure that all stride information is explicitly available
1197 * in terms of known divs. This involves calling isl_basic_set_gauss,
1198 * which is also needed because affine_hull assumes its input has been gaussed,
1199 * while isl_map_affine_hull may be called on input that has not been gaussed,
1200 * in particular from initial_facet_constraint.
1201 * Similarly, align_divs may reorder some divs so that we need to
1202 * gauss the result again.
1203 * Finally, we combine the individual affine hulls into a single
1206 __isl_give isl_basic_map
*isl_map_affine_hull(__isl_take isl_map
*map
)
1208 struct isl_basic_map
*model
= NULL
;
1209 struct isl_basic_map
*hull
= NULL
;
1210 struct isl_set
*set
;
1211 isl_basic_set
*bset
;
1213 map
= isl_map_detect_equalities(map
);
1214 map
= isl_map_local_affine_hull(map
);
1215 map
= isl_map_remove_empty_parts(map
);
1216 map
= isl_map_remove_unknown_divs(map
);
1217 map
= isl_map_align_divs_internal(map
);
1223 return replace_map_by_empty_basic_map(map
);
1225 model
= isl_basic_map_copy(map
->p
[0]);
1226 set
= isl_map_underlying_set(map
);
1227 set
= isl_set_cow(set
);
1228 set
= isl_set_local_affine_hull(set
);
1233 set
->p
[0] = affine_hull(set
->p
[0], set
->p
[--set
->n
]);
1235 bset
= isl_basic_set_copy(set
->p
[0]);
1236 hull
= isl_basic_map_overlying_set(bset
, model
);
1238 hull
= isl_basic_map_simplify(hull
);
1239 return isl_basic_map_finalize(hull
);
1241 isl_basic_map_free(model
);
1246 struct isl_basic_set
*isl_set_affine_hull(struct isl_set
*set
)
1248 return bset_from_bmap(isl_map_affine_hull(set_to_map(set
)));