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
)
128 total
= 1 + isl_basic_set_n_dim(bset1
);
129 for (r
= 0; r
< row
; ++r
) {
130 if (isl_int_is_zero(bset2
->eq
[r
][col
]))
132 isl_int_gcd(b
, bset2
->eq
[r
][col
], bset1
->eq
[row
][col
]);
133 isl_int_divexact(a
, bset1
->eq
[row
][col
], b
);
134 isl_int_divexact(b
, bset2
->eq
[r
][col
], b
);
135 isl_seq_combine(bset1
->eq
[r
], a
, bset1
->eq
[r
],
136 b
, bset1
->eq
[row
], total
);
137 isl_seq_scale(bset2
->eq
[r
], bset2
->eq
[r
], a
, total
);
141 delete_row(bset1
, row
);
146 /* Make first row entries in column col of bset1 identical to
147 * those of bset2, using only these entries of the two matrices.
148 * Let t be the last row with different entries.
149 * For each row i < t, we set
150 * A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
151 * B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
153 * A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
155 static isl_bool
transform_column(
156 __isl_keep isl_basic_set
*bset1
, __isl_keep isl_basic_set
*bset2
,
157 unsigned row
, unsigned col
)
163 for (t
= row
-1; t
>= 0; --t
)
164 if (isl_int_ne(bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]))
167 return isl_bool_false
;
169 total
= 1 + isl_basic_set_n_dim(bset1
);
173 isl_int_sub(b
, bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]);
174 for (i
= 0; i
< t
; ++i
) {
175 isl_int_sub(a
, bset2
->eq
[i
][col
], bset1
->eq
[i
][col
]);
176 isl_int_gcd(g
, a
, b
);
177 isl_int_divexact(a
, a
, g
);
178 isl_int_divexact(g
, b
, g
);
179 isl_seq_combine(bset1
->eq
[i
], g
, bset1
->eq
[i
], a
, bset1
->eq
[t
],
181 isl_seq_combine(bset2
->eq
[i
], g
, bset2
->eq
[i
], a
, bset2
->eq
[t
],
187 delete_row(bset1
, t
);
188 delete_row(bset2
, t
);
189 return isl_bool_true
;
192 /* The implementation is based on Section 5.2 of Michael Karr,
193 * "Affine Relationships Among Variables of a Program",
194 * except that the echelon form we use starts from the last column
195 * and that we are dealing with integer coefficients.
197 static __isl_give isl_basic_set
*affine_hull(
198 __isl_take isl_basic_set
*bset1
, __isl_take isl_basic_set
*bset2
)
204 if (!bset1
|| !bset2
)
207 total
= 1 + isl_basic_set_n_dim(bset1
);
210 for (col
= total
-1; col
>= 0; --col
) {
211 int is_zero1
= row
>= bset1
->n_eq
||
212 isl_int_is_zero(bset1
->eq
[row
][col
]);
213 int is_zero2
= row
>= bset2
->n_eq
||
214 isl_int_is_zero(bset2
->eq
[row
][col
]);
215 if (!is_zero1
&& !is_zero2
) {
216 set_common_multiple(bset1
, bset2
, row
, col
);
218 } else if (!is_zero1
&& is_zero2
) {
219 if (construct_column(bset1
, bset2
, row
, col
) < 0)
221 } else if (is_zero1
&& !is_zero2
) {
222 if (construct_column(bset2
, bset1
, row
, col
) < 0)
227 transform
= transform_column(bset1
, bset2
, row
, col
);
234 isl_assert(bset1
->ctx
, row
== bset1
->n_eq
, goto error
);
235 isl_basic_set_free(bset2
);
236 bset1
= isl_basic_set_normalize_constraints(bset1
);
239 isl_basic_set_free(bset1
);
240 isl_basic_set_free(bset2
);
244 /* Find an integer point in the set represented by "tab"
245 * that lies outside of the equality "eq" e(x) = 0.
246 * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
247 * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
248 * The point, if found, is returned.
249 * If no point can be found, a zero-length vector is returned.
251 * Before solving an ILP problem, we first check if simply
252 * adding the normal of the constraint to one of the known
253 * integer points in the basic set represented by "tab"
254 * yields another point inside the basic set.
256 * The caller of this function ensures that the tableau is bounded or
257 * that tab->basis and tab->n_unbounded have been set appropriately.
259 static struct isl_vec
*outside_point(struct isl_tab
*tab
, isl_int
*eq
, int up
)
262 struct isl_vec
*sample
= NULL
;
263 struct isl_tab_undo
*snap
;
271 sample
= isl_vec_alloc(ctx
, 1 + dim
);
274 isl_int_set_si(sample
->el
[0], 1);
275 isl_seq_combine(sample
->el
+ 1,
276 ctx
->one
, tab
->bmap
->sample
->el
+ 1,
277 up
? ctx
->one
: ctx
->negone
, eq
+ 1, dim
);
278 if (isl_basic_map_contains(tab
->bmap
, sample
))
280 isl_vec_free(sample
);
283 snap
= isl_tab_snap(tab
);
286 isl_seq_neg(eq
, eq
, 1 + dim
);
287 isl_int_sub_ui(eq
[0], eq
[0], 1);
289 if (isl_tab_extend_cons(tab
, 1) < 0)
291 if (isl_tab_add_ineq(tab
, eq
) < 0)
294 sample
= isl_tab_sample(tab
);
296 isl_int_add_ui(eq
[0], eq
[0], 1);
298 isl_seq_neg(eq
, eq
, 1 + dim
);
300 if (sample
&& isl_tab_rollback(tab
, snap
) < 0)
305 isl_vec_free(sample
);
309 __isl_give isl_basic_set
*isl_basic_set_recession_cone(
310 __isl_take isl_basic_set
*bset
)
314 bset
= isl_basic_set_cow(bset
);
315 if (isl_basic_set_check_no_locals(bset
) < 0)
316 return isl_basic_set_free(bset
);
318 for (i
= 0; i
< bset
->n_eq
; ++i
)
319 isl_int_set_si(bset
->eq
[i
][0], 0);
321 for (i
= 0; i
< bset
->n_ineq
; ++i
)
322 isl_int_set_si(bset
->ineq
[i
][0], 0);
324 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
325 return isl_basic_set_implicit_equalities(bset
);
328 /* Move "sample" to a point that is one up (or down) from the original
329 * point in dimension "pos".
331 static void adjacent_point(__isl_keep isl_vec
*sample
, int pos
, int up
)
334 isl_int_add_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
336 isl_int_sub_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
339 /* Check if any points that are adjacent to "sample" also belong to "bset".
340 * If so, add them to "hull" and return the updated hull.
342 * Before checking whether and adjacent point belongs to "bset", we first
343 * check whether it already belongs to "hull" as this test is typically
346 static __isl_give isl_basic_set
*add_adjacent_points(
347 __isl_take isl_basic_set
*hull
, __isl_take isl_vec
*sample
,
348 __isl_keep isl_basic_set
*bset
)
356 dim
= isl_basic_set_dim(hull
, isl_dim_set
);
358 for (i
= 0; i
< dim
; ++i
) {
359 for (up
= 0; up
<= 1; ++up
) {
361 isl_basic_set
*point
;
363 adjacent_point(sample
, i
, up
);
364 contains
= isl_basic_set_contains(hull
, sample
);
368 adjacent_point(sample
, i
, !up
);
371 contains
= isl_basic_set_contains(bset
, sample
);
375 point
= isl_basic_set_from_vec(
376 isl_vec_copy(sample
));
377 hull
= affine_hull(hull
, point
);
379 adjacent_point(sample
, i
, !up
);
385 isl_vec_free(sample
);
389 isl_vec_free(sample
);
390 isl_basic_set_free(hull
);
394 /* Extend an initial (under-)approximation of the affine hull of basic
395 * set represented by the tableau "tab"
396 * by looking for points that do not satisfy one of the equalities
397 * in the current approximation and adding them to that approximation
398 * until no such points can be found any more.
400 * The caller of this function ensures that "tab" is bounded or
401 * that tab->basis and tab->n_unbounded have been set appropriately.
403 * "bset" may be either NULL or the basic set represented by "tab".
404 * If "bset" is not NULL, we check for any point we find if any
405 * of its adjacent points also belong to "bset".
407 static __isl_give isl_basic_set
*extend_affine_hull(struct isl_tab
*tab
,
408 __isl_take isl_basic_set
*hull
, __isl_keep isl_basic_set
*bset
)
418 if (isl_tab_extend_cons(tab
, 2 * dim
+ 1) < 0)
421 for (i
= 0; i
< dim
; ++i
) {
422 struct isl_vec
*sample
;
423 struct isl_basic_set
*point
;
424 for (j
= 0; j
< hull
->n_eq
; ++j
) {
425 sample
= outside_point(tab
, hull
->eq
[j
], 1);
428 if (sample
->size
> 0)
430 isl_vec_free(sample
);
431 sample
= outside_point(tab
, hull
->eq
[j
], 0);
434 if (sample
->size
> 0)
436 isl_vec_free(sample
);
438 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
444 isl_tab_add_sample(tab
, isl_vec_copy(sample
)) < 0)
445 hull
= isl_basic_set_free(hull
);
447 hull
= add_adjacent_points(hull
, isl_vec_copy(sample
),
449 point
= isl_basic_set_from_vec(sample
);
450 hull
= affine_hull(hull
, point
);
457 isl_basic_set_free(hull
);
461 /* Construct an initial underapproximation of the hull of "bset"
462 * from "sample" and any of its adjacent points that also belong to "bset".
464 static __isl_give isl_basic_set
*initialize_hull(__isl_keep isl_basic_set
*bset
,
465 __isl_take isl_vec
*sample
)
469 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
470 hull
= add_adjacent_points(hull
, sample
, bset
);
475 /* Look for all equalities satisfied by the integer points in bset,
476 * which is assumed to be bounded.
478 * The equalities are obtained by successively looking for
479 * a point that is affinely independent of the points found so far.
480 * In particular, for each equality satisfied by the points so far,
481 * we check if there is any point on a hyperplane parallel to the
482 * corresponding hyperplane shifted by at least one (in either direction).
484 static __isl_give isl_basic_set
*uset_affine_hull_bounded(
485 __isl_take isl_basic_set
*bset
)
487 struct isl_vec
*sample
= NULL
;
488 struct isl_basic_set
*hull
;
489 struct isl_tab
*tab
= NULL
;
492 if (isl_basic_set_plain_is_empty(bset
))
495 dim
= isl_basic_set_n_dim(bset
);
497 if (bset
->sample
&& bset
->sample
->size
== 1 + dim
) {
498 int contains
= isl_basic_set_contains(bset
, bset
->sample
);
504 sample
= isl_vec_copy(bset
->sample
);
506 isl_vec_free(bset
->sample
);
511 tab
= isl_tab_from_basic_set(bset
, 1);
516 isl_vec_free(sample
);
517 return isl_basic_set_set_to_empty(bset
);
521 struct isl_tab_undo
*snap
;
522 snap
= isl_tab_snap(tab
);
523 sample
= isl_tab_sample(tab
);
524 if (isl_tab_rollback(tab
, snap
) < 0)
526 isl_vec_free(tab
->bmap
->sample
);
527 tab
->bmap
->sample
= isl_vec_copy(sample
);
532 if (sample
->size
== 0) {
534 isl_vec_free(sample
);
535 return isl_basic_set_set_to_empty(bset
);
538 hull
= initialize_hull(bset
, sample
);
540 hull
= extend_affine_hull(tab
, hull
, bset
);
541 isl_basic_set_free(bset
);
546 isl_vec_free(sample
);
548 isl_basic_set_free(bset
);
552 /* Given an unbounded tableau and an integer point satisfying the tableau,
553 * construct an initial affine hull containing the recession cone
554 * shifted to the given point.
556 * The unbounded directions are taken from the last rows of the basis,
557 * which is assumed to have been initialized appropriately.
559 static __isl_give isl_basic_set
*initial_hull(struct isl_tab
*tab
,
560 __isl_take isl_vec
*vec
)
564 struct isl_basic_set
*bset
= NULL
;
571 isl_assert(ctx
, vec
->size
!= 0, goto error
);
573 bset
= isl_basic_set_alloc(ctx
, 0, vec
->size
- 1, 0, vec
->size
- 1, 0);
576 dim
= isl_basic_set_n_dim(bset
) - tab
->n_unbounded
;
577 for (i
= 0; i
< dim
; ++i
) {
578 k
= isl_basic_set_alloc_equality(bset
);
581 isl_seq_cpy(bset
->eq
[k
] + 1, tab
->basis
->row
[1 + i
] + 1,
583 isl_seq_inner_product(bset
->eq
[k
] + 1, vec
->el
+1,
584 vec
->size
- 1, &bset
->eq
[k
][0]);
585 isl_int_neg(bset
->eq
[k
][0], bset
->eq
[k
][0]);
588 bset
= isl_basic_set_gauss(bset
, NULL
);
592 isl_basic_set_free(bset
);
597 /* Given a tableau of a set and a tableau of the corresponding
598 * recession cone, detect and add all equalities to the tableau.
599 * If the tableau is bounded, then we can simply keep the
600 * tableau in its state after the return from extend_affine_hull.
601 * However, if the tableau is unbounded, then
602 * isl_tab_set_initial_basis_with_cone will add some additional
603 * constraints to the tableau that have to be removed again.
604 * In this case, we therefore rollback to the state before
605 * any constraints were added and then add the equalities back in.
607 struct isl_tab
*isl_tab_detect_equalities(struct isl_tab
*tab
,
608 struct isl_tab
*tab_cone
)
611 struct isl_vec
*sample
;
612 struct isl_basic_set
*hull
= NULL
;
613 struct isl_tab_undo
*snap
;
615 if (!tab
|| !tab_cone
)
618 snap
= isl_tab_snap(tab
);
620 isl_mat_free(tab
->basis
);
623 isl_assert(tab
->mat
->ctx
, tab
->bmap
, goto error
);
624 isl_assert(tab
->mat
->ctx
, tab
->samples
, goto error
);
625 isl_assert(tab
->mat
->ctx
, tab
->samples
->n_col
== 1 + tab
->n_var
, goto error
);
626 isl_assert(tab
->mat
->ctx
, tab
->n_sample
> tab
->n_outside
, goto error
);
628 if (isl_tab_set_initial_basis_with_cone(tab
, tab_cone
) < 0)
631 sample
= isl_vec_alloc(tab
->mat
->ctx
, 1 + tab
->n_var
);
635 isl_seq_cpy(sample
->el
, tab
->samples
->row
[tab
->n_outside
], sample
->size
);
637 isl_vec_free(tab
->bmap
->sample
);
638 tab
->bmap
->sample
= isl_vec_copy(sample
);
640 if (tab
->n_unbounded
== 0)
641 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
643 hull
= initial_hull(tab
, isl_vec_copy(sample
));
645 for (j
= tab
->n_outside
+ 1; j
< tab
->n_sample
; ++j
) {
646 isl_seq_cpy(sample
->el
, tab
->samples
->row
[j
], sample
->size
);
647 hull
= affine_hull(hull
,
648 isl_basic_set_from_vec(isl_vec_copy(sample
)));
651 isl_vec_free(sample
);
653 hull
= extend_affine_hull(tab
, hull
, NULL
);
657 if (tab
->n_unbounded
== 0) {
658 isl_basic_set_free(hull
);
662 if (isl_tab_rollback(tab
, snap
) < 0)
665 if (hull
->n_eq
> tab
->n_zero
) {
666 for (j
= 0; j
< hull
->n_eq
; ++j
) {
667 isl_seq_normalize(tab
->mat
->ctx
, hull
->eq
[j
], 1 + tab
->n_var
);
668 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
673 isl_basic_set_free(hull
);
677 isl_basic_set_free(hull
);
682 /* Compute the affine hull of "bset", where "cone" is the recession cone
685 * We first compute a unimodular transformation that puts the unbounded
686 * directions in the last dimensions. In particular, we take a transformation
687 * that maps all equalities to equalities (in HNF) on the first dimensions.
688 * Let x be the original dimensions and y the transformed, with y_1 bounded
691 * [ y_1 ] [ y_1 ] [ Q_1 ]
692 * x = U [ y_2 ] [ y_2 ] = [ Q_2 ] x
694 * Let's call the input basic set S. We compute S' = preimage(S, U)
695 * and drop the final dimensions including any constraints involving them.
696 * This results in set S''.
697 * Then we compute the affine hull A'' of S''.
698 * Let F y_1 >= g be the constraint system of A''. In the transformed
699 * space the y_2 are unbounded, so we can add them back without any constraints,
703 * [ F 0 ] [ y_2 ] >= g
706 * [ F 0 ] [ Q_2 ] x >= g
710 * The affine hull in the original space is then obtained as
711 * A = preimage(A'', Q_1).
713 static struct isl_basic_set
*affine_hull_with_cone(struct isl_basic_set
*bset
,
714 struct isl_basic_set
*cone
)
718 struct isl_basic_set
*hull
;
719 struct isl_mat
*M
, *U
, *Q
;
724 total
= isl_basic_set_total_dim(cone
);
725 cone_dim
= total
- cone
->n_eq
;
727 M
= isl_mat_sub_alloc6(bset
->ctx
, cone
->eq
, 0, cone
->n_eq
, 1, total
);
728 M
= isl_mat_left_hermite(M
, 0, &U
, &Q
);
733 U
= isl_mat_lin_to_aff(U
);
734 bset
= isl_basic_set_preimage(bset
, isl_mat_copy(U
));
736 bset
= isl_basic_set_drop_constraints_involving(bset
, total
- cone_dim
,
738 bset
= isl_basic_set_drop_dims(bset
, total
- cone_dim
, cone_dim
);
740 Q
= isl_mat_lin_to_aff(Q
);
741 Q
= isl_mat_drop_rows(Q
, 1 + total
- cone_dim
, cone_dim
);
743 if (bset
&& bset
->sample
&& bset
->sample
->size
== 1 + total
)
744 bset
->sample
= isl_mat_vec_product(isl_mat_copy(Q
), bset
->sample
);
746 hull
= uset_affine_hull_bounded(bset
);
752 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
753 U
= isl_mat_drop_cols(U
, 1 + total
- cone_dim
, cone_dim
);
754 if (sample
&& sample
->size
> 0)
755 sample
= isl_mat_vec_product(U
, sample
);
758 hull
= isl_basic_set_preimage(hull
, Q
);
760 isl_vec_free(hull
->sample
);
761 hull
->sample
= sample
;
763 isl_vec_free(sample
);
766 isl_basic_set_free(cone
);
770 isl_basic_set_free(bset
);
771 isl_basic_set_free(cone
);
775 /* Look for all equalities satisfied by the integer points in bset,
776 * which is assumed not to have any explicit equalities.
778 * The equalities are obtained by successively looking for
779 * a point that is affinely independent of the points found so far.
780 * In particular, for each equality satisfied by the points so far,
781 * we check if there is any point on a hyperplane parallel to the
782 * corresponding hyperplane shifted by at least one (in either direction).
784 * Before looking for any outside points, we first compute the recession
785 * cone. The directions of this recession cone will always be part
786 * of the affine hull, so there is no need for looking for any points
787 * in these directions.
788 * In particular, if the recession cone is full-dimensional, then
789 * the affine hull is simply the whole universe.
791 static struct isl_basic_set
*uset_affine_hull(struct isl_basic_set
*bset
)
793 struct isl_basic_set
*cone
;
795 if (isl_basic_set_plain_is_empty(bset
))
798 cone
= isl_basic_set_recession_cone(isl_basic_set_copy(bset
));
801 if (cone
->n_eq
== 0) {
803 space
= isl_basic_set_get_space(bset
);
804 isl_basic_set_free(cone
);
805 isl_basic_set_free(bset
);
806 return isl_basic_set_universe(space
);
809 if (cone
->n_eq
< isl_basic_set_total_dim(cone
))
810 return affine_hull_with_cone(bset
, cone
);
812 isl_basic_set_free(cone
);
813 return uset_affine_hull_bounded(bset
);
815 isl_basic_set_free(bset
);
819 /* Look for all equalities satisfied by the integer points in bmap
820 * that are independent of the equalities already explicitly available
823 * We first remove all equalities already explicitly available,
824 * then look for additional equalities in the reduced space
825 * and then transform the result to the original space.
826 * The original equalities are _not_ added to this set. This is
827 * the responsibility of the calling function.
828 * The resulting basic set has all meaning about the dimensions removed.
829 * In particular, dimensions that correspond to existential variables
830 * in bmap and that are found to be fixed are not removed.
832 static struct isl_basic_set
*equalities_in_underlying_set(
833 struct isl_basic_map
*bmap
)
835 struct isl_mat
*T1
= NULL
;
836 struct isl_mat
*T2
= NULL
;
837 struct isl_basic_set
*bset
= NULL
;
838 struct isl_basic_set
*hull
= NULL
;
840 bset
= isl_basic_map_underlying_set(bmap
);
844 bset
= isl_basic_set_remove_equalities(bset
, &T1
, &T2
);
848 hull
= uset_affine_hull(bset
);
856 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
857 if (sample
&& sample
->size
> 0)
858 sample
= isl_mat_vec_product(T1
, sample
);
861 hull
= isl_basic_set_preimage(hull
, T2
);
863 isl_vec_free(hull
->sample
);
864 hull
->sample
= sample
;
866 isl_vec_free(sample
);
873 isl_basic_set_free(bset
);
874 isl_basic_set_free(hull
);
878 /* Detect and make explicit all equalities satisfied by the (integer)
881 __isl_give isl_basic_map
*isl_basic_map_detect_equalities(
882 __isl_take isl_basic_map
*bmap
)
885 struct isl_basic_set
*hull
= NULL
;
889 if (bmap
->n_ineq
== 0)
891 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
893 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_ALL_EQUALITIES
))
895 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_RATIONAL
))
896 return isl_basic_map_implicit_equalities(bmap
);
898 hull
= equalities_in_underlying_set(isl_basic_map_copy(bmap
));
901 if (ISL_F_ISSET(hull
, ISL_BASIC_SET_EMPTY
)) {
902 isl_basic_set_free(hull
);
903 return isl_basic_map_set_to_empty(bmap
);
905 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
), 0,
907 for (i
= 0; i
< hull
->n_eq
; ++i
) {
908 j
= isl_basic_map_alloc_equality(bmap
);
911 isl_seq_cpy(bmap
->eq
[j
], hull
->eq
[i
],
912 1 + isl_basic_set_total_dim(hull
));
914 isl_vec_free(bmap
->sample
);
915 bmap
->sample
= isl_vec_copy(hull
->sample
);
916 isl_basic_set_free(hull
);
917 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
| ISL_BASIC_MAP_ALL_EQUALITIES
);
918 bmap
= isl_basic_map_simplify(bmap
);
919 return isl_basic_map_finalize(bmap
);
921 isl_basic_set_free(hull
);
922 isl_basic_map_free(bmap
);
926 __isl_give isl_basic_set
*isl_basic_set_detect_equalities(
927 __isl_take isl_basic_set
*bset
)
929 return bset_from_bmap(
930 isl_basic_map_detect_equalities(bset_to_bmap(bset
)));
933 __isl_give isl_map
*isl_map_detect_equalities(__isl_take isl_map
*map
)
935 return isl_map_inline_foreach_basic_map(map
,
936 &isl_basic_map_detect_equalities
);
939 __isl_give isl_set
*isl_set_detect_equalities(__isl_take isl_set
*set
)
941 return set_from_map(isl_map_detect_equalities(set_to_map(set
)));
944 /* Return the superset of "bmap" described by the equalities
945 * satisfied by "bmap" that are already known.
947 __isl_give isl_basic_map
*isl_basic_map_plain_affine_hull(
948 __isl_take isl_basic_map
*bmap
)
950 bmap
= isl_basic_map_cow(bmap
);
952 isl_basic_map_free_inequality(bmap
, bmap
->n_ineq
);
953 bmap
= isl_basic_map_finalize(bmap
);
957 /* Return the superset of "bset" described by the equalities
958 * satisfied by "bset" that are already known.
960 __isl_give isl_basic_set
*isl_basic_set_plain_affine_hull(
961 __isl_take isl_basic_set
*bset
)
963 return isl_basic_map_plain_affine_hull(bset
);
966 /* After computing the rational affine hull (by detecting the implicit
967 * equalities), we compute the additional equalities satisfied by
968 * the integer points (if any) and add the original equalities back in.
970 __isl_give isl_basic_map
*isl_basic_map_affine_hull(
971 __isl_take isl_basic_map
*bmap
)
973 bmap
= isl_basic_map_detect_equalities(bmap
);
974 bmap
= isl_basic_map_plain_affine_hull(bmap
);
978 struct isl_basic_set
*isl_basic_set_affine_hull(struct isl_basic_set
*bset
)
980 return bset_from_bmap(isl_basic_map_affine_hull(bset_to_bmap(bset
)));
983 /* Given a rational affine matrix "M", add stride constraints to "bmap"
988 * is an integer vector. The variables x include all the variables
989 * of "bmap" except the unknown divs.
991 * If d is the common denominator of M, then we need to impose that
997 * exists alpha : d M(x) = d alpha
999 * This function is similar to add_strides in isl_morph.c
1001 static __isl_give isl_basic_map
*add_strides(__isl_take isl_basic_map
*bmap
,
1002 __isl_keep isl_mat
*M
, int n_known
)
1007 if (isl_int_is_one(M
->row
[0][0]))
1010 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
),
1011 M
->n_row
- 1, M
->n_row
- 1, 0);
1014 for (i
= 1; i
< M
->n_row
; ++i
) {
1015 isl_seq_gcd(M
->row
[i
], M
->n_col
, &gcd
);
1016 if (isl_int_is_divisible_by(gcd
, M
->row
[0][0]))
1018 div
= isl_basic_map_alloc_div(bmap
);
1021 isl_int_set_si(bmap
->div
[div
][0], 0);
1022 k
= isl_basic_map_alloc_equality(bmap
);
1025 isl_seq_cpy(bmap
->eq
[k
], M
->row
[i
], M
->n_col
);
1026 isl_seq_clr(bmap
->eq
[k
] + M
->n_col
, bmap
->n_div
- n_known
);
1027 isl_int_set(bmap
->eq
[k
][M
->n_col
- n_known
+ div
],
1035 isl_basic_map_free(bmap
);
1039 /* If there are any equalities that involve (multiple) unknown divs,
1040 * then extract the stride information encoded by those equalities
1041 * and make it explicitly available in "bmap".
1043 * We first sort the divs so that the unknown divs appear last and
1044 * then we count how many equalities involve these divs.
1046 * Let these equalities be of the form
1050 * where y represents the unknown divs and x the remaining variables.
1051 * Let [H 0] be the Hermite Normal Form of B, i.e.,
1055 * Then x is a solution of the equalities iff
1057 * H^-1 A(x) (= - [I 0] Q y)
1059 * is an integer vector. Let d be the common denominator of H^-1.
1062 * d H^-1 A(x) = d alpha
1064 * in add_strides, with alpha fresh existentially quantified variables.
1066 static __isl_give isl_basic_map
*isl_basic_map_make_strides_explicit(
1067 __isl_take isl_basic_map
*bmap
)
1076 known
= isl_basic_map_divs_known(bmap
);
1078 return isl_basic_map_free(bmap
);
1081 bmap
= isl_basic_map_sort_divs(bmap
);
1082 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1086 for (n_known
= 0; n_known
< bmap
->n_div
; ++n_known
)
1087 if (isl_int_is_zero(bmap
->div
[n_known
][0]))
1089 ctx
= isl_basic_map_get_ctx(bmap
);
1090 total
= isl_space_dim(bmap
->dim
, isl_dim_all
);
1091 for (n
= 0; n
< bmap
->n_eq
; ++n
)
1092 if (isl_seq_first_non_zero(bmap
->eq
[n
] + 1 + total
+ n_known
,
1093 bmap
->n_div
- n_known
) == -1)
1097 B
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 0, 1 + total
+ n_known
);
1098 n_col
= bmap
->n_div
- n_known
;
1099 A
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 1 + total
+ n_known
, n_col
);
1100 A
= isl_mat_left_hermite(A
, 0, NULL
, NULL
);
1101 A
= isl_mat_drop_cols(A
, n
, n_col
- n
);
1102 A
= isl_mat_lin_to_aff(A
);
1103 A
= isl_mat_right_inverse(A
);
1104 B
= isl_mat_insert_zero_rows(B
, 0, 1);
1105 B
= isl_mat_set_element_si(B
, 0, 0, 1);
1106 M
= isl_mat_product(A
, B
);
1108 return isl_basic_map_free(bmap
);
1109 bmap
= add_strides(bmap
, M
, n_known
);
1110 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1116 /* Compute the affine hull of each basic map in "map" separately
1117 * and make all stride information explicit so that we can remove
1118 * all unknown divs without losing this information.
1119 * The result is also guaranteed to be gaussed.
1121 * In simple cases where a div is determined by an equality,
1122 * calling isl_basic_map_gauss is enough to make the stride information
1123 * explicit, as it will derive an explicit representation for the div
1124 * from the equality. If, however, the stride information
1125 * is encoded through multiple unknown divs then we need to make
1126 * some extra effort in isl_basic_map_make_strides_explicit.
1128 static __isl_give isl_map
*isl_map_local_affine_hull(__isl_take isl_map
*map
)
1132 map
= isl_map_cow(map
);
1136 for (i
= 0; i
< map
->n
; ++i
) {
1137 map
->p
[i
] = isl_basic_map_affine_hull(map
->p
[i
]);
1138 map
->p
[i
] = isl_basic_map_gauss(map
->p
[i
], NULL
);
1139 map
->p
[i
] = isl_basic_map_make_strides_explicit(map
->p
[i
]);
1141 return isl_map_free(map
);
1147 static __isl_give isl_set
*isl_set_local_affine_hull(__isl_take isl_set
*set
)
1149 return isl_map_local_affine_hull(set
);
1152 /* Return an empty basic map living in the same space as "map".
1154 static __isl_give isl_basic_map
*replace_map_by_empty_basic_map(
1155 __isl_take isl_map
*map
)
1159 space
= isl_map_get_space(map
);
1161 return isl_basic_map_empty(space
);
1164 /* Compute the affine hull of "map".
1166 * We first compute the affine hull of each basic map separately.
1167 * Then we align the divs and recompute the affine hulls of the basic
1168 * maps since some of them may now have extra divs.
1169 * In order to avoid performing parametric integer programming to
1170 * compute explicit expressions for the divs, possible leading to
1171 * an explosion in the number of basic maps, we first drop all unknown
1172 * divs before aligning the divs. Note that isl_map_local_affine_hull tries
1173 * to make sure that all stride information is explicitly available
1174 * in terms of known divs. This involves calling isl_basic_set_gauss,
1175 * which is also needed because affine_hull assumes its input has been gaussed,
1176 * while isl_map_affine_hull may be called on input that has not been gaussed,
1177 * in particular from initial_facet_constraint.
1178 * Similarly, align_divs may reorder some divs so that we need to
1179 * gauss the result again.
1180 * Finally, we combine the individual affine hulls into a single
1183 __isl_give isl_basic_map
*isl_map_affine_hull(__isl_take isl_map
*map
)
1185 struct isl_basic_map
*model
= NULL
;
1186 struct isl_basic_map
*hull
= NULL
;
1187 struct isl_set
*set
;
1188 isl_basic_set
*bset
;
1190 map
= isl_map_detect_equalities(map
);
1191 map
= isl_map_local_affine_hull(map
);
1192 map
= isl_map_remove_empty_parts(map
);
1193 map
= isl_map_remove_unknown_divs(map
);
1194 map
= isl_map_align_divs_internal(map
);
1200 return replace_map_by_empty_basic_map(map
);
1202 model
= isl_basic_map_copy(map
->p
[0]);
1203 set
= isl_map_underlying_set(map
);
1204 set
= isl_set_cow(set
);
1205 set
= isl_set_local_affine_hull(set
);
1210 set
->p
[0] = affine_hull(set
->p
[0], set
->p
[--set
->n
]);
1212 bset
= isl_basic_set_copy(set
->p
[0]);
1213 hull
= isl_basic_map_overlying_set(bset
, model
);
1215 hull
= isl_basic_map_simplify(hull
);
1216 return isl_basic_map_finalize(hull
);
1218 isl_basic_map_free(model
);
1223 struct isl_basic_set
*isl_set_affine_hull(struct isl_set
*set
)
1225 return bset_from_bmap(isl_map_affine_hull(set_to_map(set
)));