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 struct isl_basic_map
*isl_basic_map_implicit_equalities(
33 struct 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 struct isl_map
*isl_map_implicit_equalities(struct isl_map
*map
)
76 for (i
= 0; i
< map
->n
; ++i
) {
77 map
->p
[i
] = isl_basic_map_implicit_equalities(map
->p
[i
]);
88 /* Make eq[row][col] of both bmaps equal so we can add the row
89 * add the column to the common matrix.
90 * Note that because of the echelon form, the columns of row row
91 * after column col are zero.
93 static void set_common_multiple(
94 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
95 unsigned row
, unsigned col
)
99 if (isl_int_eq(bset1
->eq
[row
][col
], bset2
->eq
[row
][col
]))
104 isl_int_lcm(m
, bset1
->eq
[row
][col
], bset2
->eq
[row
][col
]);
105 isl_int_divexact(c
, m
, bset1
->eq
[row
][col
]);
106 isl_seq_scale(bset1
->eq
[row
], bset1
->eq
[row
], c
, col
+1);
107 isl_int_divexact(c
, m
, bset2
->eq
[row
][col
]);
108 isl_seq_scale(bset2
->eq
[row
], bset2
->eq
[row
], c
, col
+1);
113 /* Delete a given equality, moving all the following equalities one up.
115 static void delete_row(struct isl_basic_set
*bset
, unsigned row
)
122 for (r
= row
; r
< bset
->n_eq
; ++r
)
123 bset
->eq
[r
] = bset
->eq
[r
+1];
124 bset
->eq
[bset
->n_eq
] = t
;
127 /* Make first row entries in column col of bset1 identical to
128 * those of bset2, using the fact that entry bset1->eq[row][col]=a
129 * is non-zero. Initially, these elements of bset1 are all zero.
130 * For each row i < row, we set
131 * A[i] = a * A[i] + B[i][col] * A[row]
134 * A[i][col] = B[i][col] = a * old(B[i][col])
136 static void construct_column(
137 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
138 unsigned row
, unsigned col
)
147 total
= 1 + isl_basic_set_n_dim(bset1
);
148 for (r
= 0; r
< row
; ++r
) {
149 if (isl_int_is_zero(bset2
->eq
[r
][col
]))
151 isl_int_gcd(b
, bset2
->eq
[r
][col
], bset1
->eq
[row
][col
]);
152 isl_int_divexact(a
, bset1
->eq
[row
][col
], b
);
153 isl_int_divexact(b
, bset2
->eq
[r
][col
], b
);
154 isl_seq_combine(bset1
->eq
[r
], a
, bset1
->eq
[r
],
155 b
, bset1
->eq
[row
], total
);
156 isl_seq_scale(bset2
->eq
[r
], bset2
->eq
[r
], a
, total
);
160 delete_row(bset1
, row
);
163 /* Make first row entries in column col of bset1 identical to
164 * those of bset2, using only these entries of the two matrices.
165 * Let t be the last row with different entries.
166 * For each row i < t, we set
167 * A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
168 * B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
170 * A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
172 static int transform_column(
173 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
174 unsigned row
, unsigned col
)
180 for (t
= row
-1; t
>= 0; --t
)
181 if (isl_int_ne(bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]))
186 total
= 1 + isl_basic_set_n_dim(bset1
);
190 isl_int_sub(b
, bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]);
191 for (i
= 0; i
< t
; ++i
) {
192 isl_int_sub(a
, bset2
->eq
[i
][col
], bset1
->eq
[i
][col
]);
193 isl_int_gcd(g
, a
, b
);
194 isl_int_divexact(a
, a
, g
);
195 isl_int_divexact(g
, b
, g
);
196 isl_seq_combine(bset1
->eq
[i
], g
, bset1
->eq
[i
], a
, bset1
->eq
[t
],
198 isl_seq_combine(bset2
->eq
[i
], g
, bset2
->eq
[i
], a
, bset2
->eq
[t
],
204 delete_row(bset1
, t
);
205 delete_row(bset2
, t
);
209 /* The implementation is based on Section 5.2 of Michael Karr,
210 * "Affine Relationships Among Variables of a Program",
211 * except that the echelon form we use starts from the last column
212 * and that we are dealing with integer coefficients.
214 static struct isl_basic_set
*affine_hull(
215 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
)
221 if (!bset1
|| !bset2
)
224 total
= 1 + isl_basic_set_n_dim(bset1
);
227 for (col
= total
-1; col
>= 0; --col
) {
228 int is_zero1
= row
>= bset1
->n_eq
||
229 isl_int_is_zero(bset1
->eq
[row
][col
]);
230 int is_zero2
= row
>= bset2
->n_eq
||
231 isl_int_is_zero(bset2
->eq
[row
][col
]);
232 if (!is_zero1
&& !is_zero2
) {
233 set_common_multiple(bset1
, bset2
, row
, col
);
235 } else if (!is_zero1
&& is_zero2
) {
236 construct_column(bset1
, bset2
, row
, col
);
237 } else if (is_zero1
&& !is_zero2
) {
238 construct_column(bset2
, bset1
, row
, col
);
240 if (transform_column(bset1
, bset2
, row
, col
))
244 isl_assert(bset1
->ctx
, row
== bset1
->n_eq
, goto error
);
245 isl_basic_set_free(bset2
);
246 bset1
= isl_basic_set_normalize_constraints(bset1
);
249 isl_basic_set_free(bset1
);
250 isl_basic_set_free(bset2
);
254 /* Find an integer point in the set represented by "tab"
255 * that lies outside of the equality "eq" e(x) = 0.
256 * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
257 * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
258 * The point, if found, is returned.
259 * If no point can be found, a zero-length vector is returned.
261 * Before solving an ILP problem, we first check if simply
262 * adding the normal of the constraint to one of the known
263 * integer points in the basic set represented by "tab"
264 * yields another point inside the basic set.
266 * The caller of this function ensures that the tableau is bounded or
267 * that tab->basis and tab->n_unbounded have been set appropriately.
269 static struct isl_vec
*outside_point(struct isl_tab
*tab
, isl_int
*eq
, int up
)
272 struct isl_vec
*sample
= NULL
;
273 struct isl_tab_undo
*snap
;
281 sample
= isl_vec_alloc(ctx
, 1 + dim
);
284 isl_int_set_si(sample
->el
[0], 1);
285 isl_seq_combine(sample
->el
+ 1,
286 ctx
->one
, tab
->bmap
->sample
->el
+ 1,
287 up
? ctx
->one
: ctx
->negone
, eq
+ 1, dim
);
288 if (isl_basic_map_contains(tab
->bmap
, sample
))
290 isl_vec_free(sample
);
293 snap
= isl_tab_snap(tab
);
296 isl_seq_neg(eq
, eq
, 1 + dim
);
297 isl_int_sub_ui(eq
[0], eq
[0], 1);
299 if (isl_tab_extend_cons(tab
, 1) < 0)
301 if (isl_tab_add_ineq(tab
, eq
) < 0)
304 sample
= isl_tab_sample(tab
);
306 isl_int_add_ui(eq
[0], eq
[0], 1);
308 isl_seq_neg(eq
, eq
, 1 + dim
);
310 if (sample
&& isl_tab_rollback(tab
, snap
) < 0)
315 isl_vec_free(sample
);
319 struct isl_basic_set
*isl_basic_set_recession_cone(struct isl_basic_set
*bset
)
323 bset
= isl_basic_set_cow(bset
);
326 isl_assert(bset
->ctx
, bset
->n_div
== 0, goto error
);
328 for (i
= 0; i
< bset
->n_eq
; ++i
)
329 isl_int_set_si(bset
->eq
[i
][0], 0);
331 for (i
= 0; i
< bset
->n_ineq
; ++i
)
332 isl_int_set_si(bset
->ineq
[i
][0], 0);
334 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
335 return isl_basic_set_implicit_equalities(bset
);
337 isl_basic_set_free(bset
);
341 __isl_give isl_set
*isl_set_recession_cone(__isl_take isl_set
*set
)
350 set
= isl_set_remove_divs(set
);
351 set
= isl_set_cow(set
);
355 for (i
= 0; i
< set
->n
; ++i
) {
356 set
->p
[i
] = isl_basic_set_recession_cone(set
->p
[i
]);
367 /* Move "sample" to a point that is one up (or down) from the original
368 * point in dimension "pos".
370 static void adjacent_point(__isl_keep isl_vec
*sample
, int pos
, int up
)
373 isl_int_add_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
375 isl_int_sub_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
378 /* Check if any points that are adjacent to "sample" also belong to "bset".
379 * If so, add them to "hull" and return the updated hull.
381 * Before checking whether and adjacent point belongs to "bset", we first
382 * check whether it already belongs to "hull" as this test is typically
385 static __isl_give isl_basic_set
*add_adjacent_points(
386 __isl_take isl_basic_set
*hull
, __isl_take isl_vec
*sample
,
387 __isl_keep isl_basic_set
*bset
)
395 dim
= isl_basic_set_dim(hull
, isl_dim_set
);
397 for (i
= 0; i
< dim
; ++i
) {
398 for (up
= 0; up
<= 1; ++up
) {
400 isl_basic_set
*point
;
402 adjacent_point(sample
, i
, up
);
403 contains
= isl_basic_set_contains(hull
, sample
);
407 adjacent_point(sample
, i
, !up
);
410 contains
= isl_basic_set_contains(bset
, sample
);
414 point
= isl_basic_set_from_vec(
415 isl_vec_copy(sample
));
416 hull
= affine_hull(hull
, point
);
418 adjacent_point(sample
, i
, !up
);
424 isl_vec_free(sample
);
428 isl_vec_free(sample
);
429 isl_basic_set_free(hull
);
433 /* Extend an initial (under-)approximation of the affine hull of basic
434 * set represented by the tableau "tab"
435 * by looking for points that do not satisfy one of the equalities
436 * in the current approximation and adding them to that approximation
437 * until no such points can be found any more.
439 * The caller of this function ensures that "tab" is bounded or
440 * that tab->basis and tab->n_unbounded have been set appropriately.
442 * "bset" may be either NULL or the basic set represented by "tab".
443 * If "bset" is not NULL, we check for any point we find if any
444 * of its adjacent points also belong to "bset".
446 static __isl_give isl_basic_set
*extend_affine_hull(struct isl_tab
*tab
,
447 __isl_take isl_basic_set
*hull
, __isl_keep isl_basic_set
*bset
)
457 if (isl_tab_extend_cons(tab
, 2 * dim
+ 1) < 0)
460 for (i
= 0; i
< dim
; ++i
) {
461 struct isl_vec
*sample
;
462 struct isl_basic_set
*point
;
463 for (j
= 0; j
< hull
->n_eq
; ++j
) {
464 sample
= outside_point(tab
, hull
->eq
[j
], 1);
467 if (sample
->size
> 0)
469 isl_vec_free(sample
);
470 sample
= outside_point(tab
, hull
->eq
[j
], 0);
473 if (sample
->size
> 0)
475 isl_vec_free(sample
);
477 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
483 isl_tab_add_sample(tab
, isl_vec_copy(sample
)) < 0)
484 hull
= isl_basic_set_free(hull
);
486 hull
= add_adjacent_points(hull
, isl_vec_copy(sample
),
488 point
= isl_basic_set_from_vec(sample
);
489 hull
= affine_hull(hull
, point
);
496 isl_basic_set_free(hull
);
500 /* Construct an initial underapproximation of the hull of "bset"
501 * from "sample" and any of its adjacent points that also belong to "bset".
503 static __isl_give isl_basic_set
*initialize_hull(__isl_keep isl_basic_set
*bset
,
504 __isl_take isl_vec
*sample
)
508 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
509 hull
= add_adjacent_points(hull
, sample
, bset
);
514 /* Look for all equalities satisfied by the integer points in bset,
515 * which is assumed to be bounded.
517 * The equalities are obtained by successively looking for
518 * a point that is affinely independent of the points found so far.
519 * In particular, for each equality satisfied by the points so far,
520 * we check if there is any point on a hyperplane parallel to the
521 * corresponding hyperplane shifted by at least one (in either direction).
523 static struct isl_basic_set
*uset_affine_hull_bounded(struct isl_basic_set
*bset
)
525 struct isl_vec
*sample
= NULL
;
526 struct isl_basic_set
*hull
;
527 struct isl_tab
*tab
= NULL
;
530 if (isl_basic_set_plain_is_empty(bset
))
533 dim
= isl_basic_set_n_dim(bset
);
535 if (bset
->sample
&& bset
->sample
->size
== 1 + dim
) {
536 int contains
= isl_basic_set_contains(bset
, bset
->sample
);
542 sample
= isl_vec_copy(bset
->sample
);
544 isl_vec_free(bset
->sample
);
549 tab
= isl_tab_from_basic_set(bset
, 1);
554 isl_vec_free(sample
);
555 return isl_basic_set_set_to_empty(bset
);
559 struct isl_tab_undo
*snap
;
560 snap
= isl_tab_snap(tab
);
561 sample
= isl_tab_sample(tab
);
562 if (isl_tab_rollback(tab
, snap
) < 0)
564 isl_vec_free(tab
->bmap
->sample
);
565 tab
->bmap
->sample
= isl_vec_copy(sample
);
570 if (sample
->size
== 0) {
572 isl_vec_free(sample
);
573 return isl_basic_set_set_to_empty(bset
);
576 hull
= initialize_hull(bset
, sample
);
578 hull
= extend_affine_hull(tab
, hull
, bset
);
579 isl_basic_set_free(bset
);
584 isl_vec_free(sample
);
586 isl_basic_set_free(bset
);
590 /* Given an unbounded tableau and an integer point satisfying the tableau,
591 * construct an initial affine hull containing the recession cone
592 * shifted to the given point.
594 * The unbounded directions are taken from the last rows of the basis,
595 * which is assumed to have been initialized appropriately.
597 static __isl_give isl_basic_set
*initial_hull(struct isl_tab
*tab
,
598 __isl_take isl_vec
*vec
)
602 struct isl_basic_set
*bset
= NULL
;
609 isl_assert(ctx
, vec
->size
!= 0, goto error
);
611 bset
= isl_basic_set_alloc(ctx
, 0, vec
->size
- 1, 0, vec
->size
- 1, 0);
614 dim
= isl_basic_set_n_dim(bset
) - tab
->n_unbounded
;
615 for (i
= 0; i
< dim
; ++i
) {
616 k
= isl_basic_set_alloc_equality(bset
);
619 isl_seq_cpy(bset
->eq
[k
] + 1, tab
->basis
->row
[1 + i
] + 1,
621 isl_seq_inner_product(bset
->eq
[k
] + 1, vec
->el
+1,
622 vec
->size
- 1, &bset
->eq
[k
][0]);
623 isl_int_neg(bset
->eq
[k
][0], bset
->eq
[k
][0]);
626 bset
= isl_basic_set_gauss(bset
, NULL
);
630 isl_basic_set_free(bset
);
635 /* Given a tableau of a set and a tableau of the corresponding
636 * recession cone, detect and add all equalities to the tableau.
637 * If the tableau is bounded, then we can simply keep the
638 * tableau in its state after the return from extend_affine_hull.
639 * However, if the tableau is unbounded, then
640 * isl_tab_set_initial_basis_with_cone will add some additional
641 * constraints to the tableau that have to be removed again.
642 * In this case, we therefore rollback to the state before
643 * any constraints were added and then add the equalities back in.
645 struct isl_tab
*isl_tab_detect_equalities(struct isl_tab
*tab
,
646 struct isl_tab
*tab_cone
)
649 struct isl_vec
*sample
;
650 struct isl_basic_set
*hull
= NULL
;
651 struct isl_tab_undo
*snap
;
653 if (!tab
|| !tab_cone
)
656 snap
= isl_tab_snap(tab
);
658 isl_mat_free(tab
->basis
);
661 isl_assert(tab
->mat
->ctx
, tab
->bmap
, goto error
);
662 isl_assert(tab
->mat
->ctx
, tab
->samples
, goto error
);
663 isl_assert(tab
->mat
->ctx
, tab
->samples
->n_col
== 1 + tab
->n_var
, goto error
);
664 isl_assert(tab
->mat
->ctx
, tab
->n_sample
> tab
->n_outside
, goto error
);
666 if (isl_tab_set_initial_basis_with_cone(tab
, tab_cone
) < 0)
669 sample
= isl_vec_alloc(tab
->mat
->ctx
, 1 + tab
->n_var
);
673 isl_seq_cpy(sample
->el
, tab
->samples
->row
[tab
->n_outside
], sample
->size
);
675 isl_vec_free(tab
->bmap
->sample
);
676 tab
->bmap
->sample
= isl_vec_copy(sample
);
678 if (tab
->n_unbounded
== 0)
679 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
681 hull
= initial_hull(tab
, isl_vec_copy(sample
));
683 for (j
= tab
->n_outside
+ 1; j
< tab
->n_sample
; ++j
) {
684 isl_seq_cpy(sample
->el
, tab
->samples
->row
[j
], sample
->size
);
685 hull
= affine_hull(hull
,
686 isl_basic_set_from_vec(isl_vec_copy(sample
)));
689 isl_vec_free(sample
);
691 hull
= extend_affine_hull(tab
, hull
, NULL
);
695 if (tab
->n_unbounded
== 0) {
696 isl_basic_set_free(hull
);
700 if (isl_tab_rollback(tab
, snap
) < 0)
703 if (hull
->n_eq
> tab
->n_zero
) {
704 for (j
= 0; j
< hull
->n_eq
; ++j
) {
705 isl_seq_normalize(tab
->mat
->ctx
, hull
->eq
[j
], 1 + tab
->n_var
);
706 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
711 isl_basic_set_free(hull
);
715 isl_basic_set_free(hull
);
720 /* Compute the affine hull of "bset", where "cone" is the recession cone
723 * We first compute a unimodular transformation that puts the unbounded
724 * directions in the last dimensions. In particular, we take a transformation
725 * that maps all equalities to equalities (in HNF) on the first dimensions.
726 * Let x be the original dimensions and y the transformed, with y_1 bounded
729 * [ y_1 ] [ y_1 ] [ Q_1 ]
730 * x = U [ y_2 ] [ y_2 ] = [ Q_2 ] x
732 * Let's call the input basic set S. We compute S' = preimage(S, U)
733 * and drop the final dimensions including any constraints involving them.
734 * This results in set S''.
735 * Then we compute the affine hull A'' of S''.
736 * Let F y_1 >= g be the constraint system of A''. In the transformed
737 * space the y_2 are unbounded, so we can add them back without any constraints,
741 * [ F 0 ] [ y_2 ] >= g
744 * [ F 0 ] [ Q_2 ] x >= g
748 * The affine hull in the original space is then obtained as
749 * A = preimage(A'', Q_1).
751 static struct isl_basic_set
*affine_hull_with_cone(struct isl_basic_set
*bset
,
752 struct isl_basic_set
*cone
)
756 struct isl_basic_set
*hull
;
757 struct isl_mat
*M
, *U
, *Q
;
762 total
= isl_basic_set_total_dim(cone
);
763 cone_dim
= total
- cone
->n_eq
;
765 M
= isl_mat_sub_alloc6(bset
->ctx
, cone
->eq
, 0, cone
->n_eq
, 1, total
);
766 M
= isl_mat_left_hermite(M
, 0, &U
, &Q
);
771 U
= isl_mat_lin_to_aff(U
);
772 bset
= isl_basic_set_preimage(bset
, isl_mat_copy(U
));
774 bset
= isl_basic_set_drop_constraints_involving(bset
, total
- cone_dim
,
776 bset
= isl_basic_set_drop_dims(bset
, total
- cone_dim
, cone_dim
);
778 Q
= isl_mat_lin_to_aff(Q
);
779 Q
= isl_mat_drop_rows(Q
, 1 + total
- cone_dim
, cone_dim
);
781 if (bset
&& bset
->sample
&& bset
->sample
->size
== 1 + total
)
782 bset
->sample
= isl_mat_vec_product(isl_mat_copy(Q
), bset
->sample
);
784 hull
= uset_affine_hull_bounded(bset
);
790 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
791 U
= isl_mat_drop_cols(U
, 1 + total
- cone_dim
, cone_dim
);
792 if (sample
&& sample
->size
> 0)
793 sample
= isl_mat_vec_product(U
, sample
);
796 hull
= isl_basic_set_preimage(hull
, Q
);
798 isl_vec_free(hull
->sample
);
799 hull
->sample
= sample
;
801 isl_vec_free(sample
);
804 isl_basic_set_free(cone
);
808 isl_basic_set_free(bset
);
809 isl_basic_set_free(cone
);
813 /* Look for all equalities satisfied by the integer points in bset,
814 * which is assumed not to have any explicit equalities.
816 * The equalities are obtained by successively looking for
817 * a point that is affinely independent of the points found so far.
818 * In particular, for each equality satisfied by the points so far,
819 * we check if there is any point on a hyperplane parallel to the
820 * corresponding hyperplane shifted by at least one (in either direction).
822 * Before looking for any outside points, we first compute the recession
823 * cone. The directions of this recession cone will always be part
824 * of the affine hull, so there is no need for looking for any points
825 * in these directions.
826 * In particular, if the recession cone is full-dimensional, then
827 * the affine hull is simply the whole universe.
829 static struct isl_basic_set
*uset_affine_hull(struct isl_basic_set
*bset
)
831 struct isl_basic_set
*cone
;
833 if (isl_basic_set_plain_is_empty(bset
))
836 cone
= isl_basic_set_recession_cone(isl_basic_set_copy(bset
));
839 if (cone
->n_eq
== 0) {
841 space
= isl_basic_set_get_space(bset
);
842 isl_basic_set_free(cone
);
843 isl_basic_set_free(bset
);
844 return isl_basic_set_universe(space
);
847 if (cone
->n_eq
< isl_basic_set_total_dim(cone
))
848 return affine_hull_with_cone(bset
, cone
);
850 isl_basic_set_free(cone
);
851 return uset_affine_hull_bounded(bset
);
853 isl_basic_set_free(bset
);
857 /* Look for all equalities satisfied by the integer points in bmap
858 * that are independent of the equalities already explicitly available
861 * We first remove all equalities already explicitly available,
862 * then look for additional equalities in the reduced space
863 * and then transform the result to the original space.
864 * The original equalities are _not_ added to this set. This is
865 * the responsibility of the calling function.
866 * The resulting basic set has all meaning about the dimensions removed.
867 * In particular, dimensions that correspond to existential variables
868 * in bmap and that are found to be fixed are not removed.
870 static struct isl_basic_set
*equalities_in_underlying_set(
871 struct isl_basic_map
*bmap
)
873 struct isl_mat
*T1
= NULL
;
874 struct isl_mat
*T2
= NULL
;
875 struct isl_basic_set
*bset
= NULL
;
876 struct isl_basic_set
*hull
= NULL
;
878 bset
= isl_basic_map_underlying_set(bmap
);
882 bset
= isl_basic_set_remove_equalities(bset
, &T1
, &T2
);
886 hull
= uset_affine_hull(bset
);
894 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
895 if (sample
&& sample
->size
> 0)
896 sample
= isl_mat_vec_product(T1
, sample
);
899 hull
= isl_basic_set_preimage(hull
, T2
);
901 isl_vec_free(hull
->sample
);
902 hull
->sample
= sample
;
904 isl_vec_free(sample
);
911 isl_basic_set_free(bset
);
912 isl_basic_set_free(hull
);
916 /* Detect and make explicit all equalities satisfied by the (integer)
919 struct isl_basic_map
*isl_basic_map_detect_equalities(
920 struct isl_basic_map
*bmap
)
923 struct isl_basic_set
*hull
= NULL
;
927 if (bmap
->n_ineq
== 0)
929 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
931 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_ALL_EQUALITIES
))
933 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_RATIONAL
))
934 return isl_basic_map_implicit_equalities(bmap
);
936 hull
= equalities_in_underlying_set(isl_basic_map_copy(bmap
));
939 if (ISL_F_ISSET(hull
, ISL_BASIC_SET_EMPTY
)) {
940 isl_basic_set_free(hull
);
941 return isl_basic_map_set_to_empty(bmap
);
943 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
), 0,
945 for (i
= 0; i
< hull
->n_eq
; ++i
) {
946 j
= isl_basic_map_alloc_equality(bmap
);
949 isl_seq_cpy(bmap
->eq
[j
], hull
->eq
[i
],
950 1 + isl_basic_set_total_dim(hull
));
952 isl_vec_free(bmap
->sample
);
953 bmap
->sample
= isl_vec_copy(hull
->sample
);
954 isl_basic_set_free(hull
);
955 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
| ISL_BASIC_MAP_ALL_EQUALITIES
);
956 bmap
= isl_basic_map_simplify(bmap
);
957 return isl_basic_map_finalize(bmap
);
959 isl_basic_set_free(hull
);
960 isl_basic_map_free(bmap
);
964 __isl_give isl_basic_set
*isl_basic_set_detect_equalities(
965 __isl_take isl_basic_set
*bset
)
967 return bset_from_bmap(
968 isl_basic_map_detect_equalities(bset_to_bmap(bset
)));
971 __isl_give isl_map
*isl_map_detect_equalities(__isl_take isl_map
*map
)
973 return isl_map_inline_foreach_basic_map(map
,
974 &isl_basic_map_detect_equalities
);
977 __isl_give isl_set
*isl_set_detect_equalities(__isl_take isl_set
*set
)
979 return set_from_map(isl_map_detect_equalities(set_to_map(set
)));
982 /* Return the superset of "bmap" described by the equalities
983 * satisfied by "bmap" that are already known.
985 __isl_give isl_basic_map
*isl_basic_map_plain_affine_hull(
986 __isl_take isl_basic_map
*bmap
)
988 bmap
= isl_basic_map_cow(bmap
);
990 isl_basic_map_free_inequality(bmap
, bmap
->n_ineq
);
991 bmap
= isl_basic_map_finalize(bmap
);
995 /* Return the superset of "bset" described by the equalities
996 * satisfied by "bset" that are already known.
998 __isl_give isl_basic_set
*isl_basic_set_plain_affine_hull(
999 __isl_take isl_basic_set
*bset
)
1001 return isl_basic_map_plain_affine_hull(bset
);
1004 /* After computing the rational affine hull (by detecting the implicit
1005 * equalities), we compute the additional equalities satisfied by
1006 * the integer points (if any) and add the original equalities back in.
1008 struct isl_basic_map
*isl_basic_map_affine_hull(struct isl_basic_map
*bmap
)
1010 bmap
= isl_basic_map_detect_equalities(bmap
);
1011 bmap
= isl_basic_map_plain_affine_hull(bmap
);
1015 struct isl_basic_set
*isl_basic_set_affine_hull(struct isl_basic_set
*bset
)
1017 return bset_from_bmap(isl_basic_map_affine_hull(bset_to_bmap(bset
)));
1020 /* Given a rational affine matrix "M", add stride constraints to "bmap"
1025 * is an integer vector. The variables x include all the variables
1026 * of "bmap" except the unknown divs.
1028 * If d is the common denominator of M, then we need to impose that
1034 * exists alpha : d M(x) = d alpha
1036 * This function is similar to add_strides in isl_morph.c
1038 static __isl_give isl_basic_map
*add_strides(__isl_take isl_basic_map
*bmap
,
1039 __isl_keep isl_mat
*M
, int n_known
)
1044 if (isl_int_is_one(M
->row
[0][0]))
1047 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
),
1048 M
->n_row
- 1, M
->n_row
- 1, 0);
1051 for (i
= 1; i
< M
->n_row
; ++i
) {
1052 isl_seq_gcd(M
->row
[i
], M
->n_col
, &gcd
);
1053 if (isl_int_is_divisible_by(gcd
, M
->row
[0][0]))
1055 div
= isl_basic_map_alloc_div(bmap
);
1058 isl_int_set_si(bmap
->div
[div
][0], 0);
1059 k
= isl_basic_map_alloc_equality(bmap
);
1062 isl_seq_cpy(bmap
->eq
[k
], M
->row
[i
], M
->n_col
);
1063 isl_seq_clr(bmap
->eq
[k
] + M
->n_col
, bmap
->n_div
- n_known
);
1064 isl_int_set(bmap
->eq
[k
][M
->n_col
- n_known
+ div
],
1072 isl_basic_map_free(bmap
);
1076 /* If there are any equalities that involve (multiple) unknown divs,
1077 * then extract the stride information encoded by those equalities
1078 * and make it explicitly available in "bmap".
1080 * We first sort the divs so that the unknown divs appear last and
1081 * then we count how many equalities involve these divs.
1083 * Let these equalities be of the form
1087 * where y represents the unknown divs and x the remaining variables.
1088 * Let [H 0] be the Hermite Normal Form of B, i.e.,
1092 * Then x is a solution of the equalities iff
1094 * H^-1 A(x) (= - [I 0] Q y)
1096 * is an integer vector. Let d be the common denominator of H^-1.
1099 * d H^-1 A(x) = d alpha
1101 * in add_strides, with alpha fresh existentially quantified variables.
1103 static __isl_give isl_basic_map
*isl_basic_map_make_strides_explicit(
1104 __isl_take isl_basic_map
*bmap
)
1113 known
= isl_basic_map_divs_known(bmap
);
1115 return isl_basic_map_free(bmap
);
1118 bmap
= isl_basic_map_sort_divs(bmap
);
1119 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1123 for (n_known
= 0; n_known
< bmap
->n_div
; ++n_known
)
1124 if (isl_int_is_zero(bmap
->div
[n_known
][0]))
1126 ctx
= isl_basic_map_get_ctx(bmap
);
1127 total
= isl_space_dim(bmap
->dim
, isl_dim_all
);
1128 for (n
= 0; n
< bmap
->n_eq
; ++n
)
1129 if (isl_seq_first_non_zero(bmap
->eq
[n
] + 1 + total
+ n_known
,
1130 bmap
->n_div
- n_known
) == -1)
1134 B
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 0, 1 + total
+ n_known
);
1135 n_col
= bmap
->n_div
- n_known
;
1136 A
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 1 + total
+ n_known
, n_col
);
1137 A
= isl_mat_left_hermite(A
, 0, NULL
, NULL
);
1138 A
= isl_mat_drop_cols(A
, n
, n_col
- n
);
1139 A
= isl_mat_lin_to_aff(A
);
1140 A
= isl_mat_right_inverse(A
);
1141 B
= isl_mat_insert_zero_rows(B
, 0, 1);
1142 B
= isl_mat_set_element_si(B
, 0, 0, 1);
1143 M
= isl_mat_product(A
, B
);
1145 return isl_basic_map_free(bmap
);
1146 bmap
= add_strides(bmap
, M
, n_known
);
1147 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1153 /* Compute the affine hull of each basic map in "map" separately
1154 * and make all stride information explicit so that we can remove
1155 * all unknown divs without losing this information.
1156 * The result is also guaranteed to be gaussed.
1158 * In simple cases where a div is determined by an equality,
1159 * calling isl_basic_map_gauss is enough to make the stride information
1160 * explicit, as it will derive an explicit representation for the div
1161 * from the equality. If, however, the stride information
1162 * is encoded through multiple unknown divs then we need to make
1163 * some extra effort in isl_basic_map_make_strides_explicit.
1165 static __isl_give isl_map
*isl_map_local_affine_hull(__isl_take isl_map
*map
)
1169 map
= isl_map_cow(map
);
1173 for (i
= 0; i
< map
->n
; ++i
) {
1174 map
->p
[i
] = isl_basic_map_affine_hull(map
->p
[i
]);
1175 map
->p
[i
] = isl_basic_map_gauss(map
->p
[i
], NULL
);
1176 map
->p
[i
] = isl_basic_map_make_strides_explicit(map
->p
[i
]);
1178 return isl_map_free(map
);
1184 static __isl_give isl_set
*isl_set_local_affine_hull(__isl_take isl_set
*set
)
1186 return isl_map_local_affine_hull(set
);
1189 /* Return an empty basic map living in the same space as "map".
1191 static __isl_give isl_basic_map
*replace_map_by_empty_basic_map(
1192 __isl_take isl_map
*map
)
1196 space
= isl_map_get_space(map
);
1198 return isl_basic_map_empty(space
);
1201 /* Compute the affine hull of "map".
1203 * We first compute the affine hull of each basic map separately.
1204 * Then we align the divs and recompute the affine hulls of the basic
1205 * maps since some of them may now have extra divs.
1206 * In order to avoid performing parametric integer programming to
1207 * compute explicit expressions for the divs, possible leading to
1208 * an explosion in the number of basic maps, we first drop all unknown
1209 * divs before aligning the divs. Note that isl_map_local_affine_hull tries
1210 * to make sure that all stride information is explicitly available
1211 * in terms of known divs. This involves calling isl_basic_set_gauss,
1212 * which is also needed because affine_hull assumes its input has been gaussed,
1213 * while isl_map_affine_hull may be called on input that has not been gaussed,
1214 * in particular from initial_facet_constraint.
1215 * Similarly, align_divs may reorder some divs so that we need to
1216 * gauss the result again.
1217 * Finally, we combine the individual affine hulls into a single
1220 __isl_give isl_basic_map
*isl_map_affine_hull(__isl_take isl_map
*map
)
1222 struct isl_basic_map
*model
= NULL
;
1223 struct isl_basic_map
*hull
= NULL
;
1224 struct isl_set
*set
;
1225 isl_basic_set
*bset
;
1227 map
= isl_map_detect_equalities(map
);
1228 map
= isl_map_local_affine_hull(map
);
1229 map
= isl_map_remove_empty_parts(map
);
1230 map
= isl_map_remove_unknown_divs(map
);
1231 map
= isl_map_align_divs(map
);
1237 return replace_map_by_empty_basic_map(map
);
1239 model
= isl_basic_map_copy(map
->p
[0]);
1240 set
= isl_map_underlying_set(map
);
1241 set
= isl_set_cow(set
);
1242 set
= isl_set_local_affine_hull(set
);
1247 set
->p
[0] = affine_hull(set
->p
[0], set
->p
[--set
->n
]);
1249 bset
= isl_basic_set_copy(set
->p
[0]);
1250 hull
= isl_basic_map_overlying_set(bset
, model
);
1252 hull
= isl_basic_map_simplify(hull
);
1253 return isl_basic_map_finalize(hull
);
1255 isl_basic_map_free(model
);
1260 struct isl_basic_set
*isl_set_affine_hull(struct isl_set
*set
)
1262 return bset_from_bmap(isl_map_affine_hull(set_to_map(set
)));