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>
26 struct isl_basic_map
*isl_basic_map_implicit_equalities(
27 struct isl_basic_map
*bmap
)
34 bmap
= isl_basic_map_gauss(bmap
, NULL
);
35 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
37 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
))
39 if (bmap
->n_ineq
<= 1)
42 tab
= isl_tab_from_basic_map(bmap
, 0);
43 if (isl_tab_detect_implicit_equalities(tab
) < 0)
45 bmap
= isl_basic_map_update_from_tab(bmap
, tab
);
47 bmap
= isl_basic_map_gauss(bmap
, NULL
);
48 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
);
52 isl_basic_map_free(bmap
);
56 struct isl_basic_set
*isl_basic_set_implicit_equalities(
57 struct isl_basic_set
*bset
)
59 return (struct isl_basic_set
*)
60 isl_basic_map_implicit_equalities((struct isl_basic_map
*)bset
);
63 struct isl_map
*isl_map_implicit_equalities(struct isl_map
*map
)
70 for (i
= 0; i
< map
->n
; ++i
) {
71 map
->p
[i
] = isl_basic_map_implicit_equalities(map
->p
[i
]);
82 /* Make eq[row][col] of both bmaps equal so we can add the row
83 * add the column to the common matrix.
84 * Note that because of the echelon form, the columns of row row
85 * after column col are zero.
87 static void set_common_multiple(
88 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
89 unsigned row
, unsigned col
)
93 if (isl_int_eq(bset1
->eq
[row
][col
], bset2
->eq
[row
][col
]))
98 isl_int_lcm(m
, bset1
->eq
[row
][col
], bset2
->eq
[row
][col
]);
99 isl_int_divexact(c
, m
, bset1
->eq
[row
][col
]);
100 isl_seq_scale(bset1
->eq
[row
], bset1
->eq
[row
], c
, col
+1);
101 isl_int_divexact(c
, m
, bset2
->eq
[row
][col
]);
102 isl_seq_scale(bset2
->eq
[row
], bset2
->eq
[row
], c
, col
+1);
107 /* Delete a given equality, moving all the following equalities one up.
109 static void delete_row(struct isl_basic_set
*bset
, unsigned row
)
116 for (r
= row
; r
< bset
->n_eq
; ++r
)
117 bset
->eq
[r
] = bset
->eq
[r
+1];
118 bset
->eq
[bset
->n_eq
] = t
;
121 /* Make first row entries in column col of bset1 identical to
122 * those of bset2, using the fact that entry bset1->eq[row][col]=a
123 * is non-zero. Initially, these elements of bset1 are all zero.
124 * For each row i < row, we set
125 * A[i] = a * A[i] + B[i][col] * A[row]
128 * A[i][col] = B[i][col] = a * old(B[i][col])
130 static void construct_column(
131 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
132 unsigned row
, unsigned col
)
141 total
= 1 + isl_basic_set_n_dim(bset1
);
142 for (r
= 0; r
< row
; ++r
) {
143 if (isl_int_is_zero(bset2
->eq
[r
][col
]))
145 isl_int_gcd(b
, bset2
->eq
[r
][col
], bset1
->eq
[row
][col
]);
146 isl_int_divexact(a
, bset1
->eq
[row
][col
], b
);
147 isl_int_divexact(b
, bset2
->eq
[r
][col
], b
);
148 isl_seq_combine(bset1
->eq
[r
], a
, bset1
->eq
[r
],
149 b
, bset1
->eq
[row
], total
);
150 isl_seq_scale(bset2
->eq
[r
], bset2
->eq
[r
], a
, total
);
154 delete_row(bset1
, row
);
157 /* Make first row entries in column col of bset1 identical to
158 * those of bset2, using only these entries of the two matrices.
159 * Let t be the last row with different entries.
160 * For each row i < t, we set
161 * A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
162 * B[i] = (A[t][col]-B[t][col]) * B[i] + (B[i][col]-A[i][col) * B[t]
164 * A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
166 static int transform_column(
167 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
,
168 unsigned row
, unsigned col
)
174 for (t
= row
-1; t
>= 0; --t
)
175 if (isl_int_ne(bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]))
180 total
= 1 + isl_basic_set_n_dim(bset1
);
184 isl_int_sub(b
, bset1
->eq
[t
][col
], bset2
->eq
[t
][col
]);
185 for (i
= 0; i
< t
; ++i
) {
186 isl_int_sub(a
, bset2
->eq
[i
][col
], bset1
->eq
[i
][col
]);
187 isl_int_gcd(g
, a
, b
);
188 isl_int_divexact(a
, a
, g
);
189 isl_int_divexact(g
, b
, g
);
190 isl_seq_combine(bset1
->eq
[i
], g
, bset1
->eq
[i
], a
, bset1
->eq
[t
],
192 isl_seq_combine(bset2
->eq
[i
], g
, bset2
->eq
[i
], a
, bset2
->eq
[t
],
198 delete_row(bset1
, t
);
199 delete_row(bset2
, t
);
203 /* The implementation is based on Section 5.2 of Michael Karr,
204 * "Affine Relationships Among Variables of a Program",
205 * except that the echelon form we use starts from the last column
206 * and that we are dealing with integer coefficients.
208 static struct isl_basic_set
*affine_hull(
209 struct isl_basic_set
*bset1
, struct isl_basic_set
*bset2
)
215 if (!bset1
|| !bset2
)
218 total
= 1 + isl_basic_set_n_dim(bset1
);
221 for (col
= total
-1; col
>= 0; --col
) {
222 int is_zero1
= row
>= bset1
->n_eq
||
223 isl_int_is_zero(bset1
->eq
[row
][col
]);
224 int is_zero2
= row
>= bset2
->n_eq
||
225 isl_int_is_zero(bset2
->eq
[row
][col
]);
226 if (!is_zero1
&& !is_zero2
) {
227 set_common_multiple(bset1
, bset2
, row
, col
);
229 } else if (!is_zero1
&& is_zero2
) {
230 construct_column(bset1
, bset2
, row
, col
);
231 } else if (is_zero1
&& !is_zero2
) {
232 construct_column(bset2
, bset1
, row
, col
);
234 if (transform_column(bset1
, bset2
, row
, col
))
238 isl_assert(bset1
->ctx
, row
== bset1
->n_eq
, goto error
);
239 isl_basic_set_free(bset2
);
240 bset1
= isl_basic_set_normalize_constraints(bset1
);
243 isl_basic_set_free(bset1
);
244 isl_basic_set_free(bset2
);
248 /* Find an integer point in the set represented by "tab"
249 * that lies outside of the equality "eq" e(x) = 0.
250 * If "up" is true, look for a point satisfying e(x) - 1 >= 0.
251 * Otherwise, look for a point satisfying -e(x) - 1 >= 0 (i.e., e(x) <= -1).
252 * The point, if found, is returned.
253 * If no point can be found, a zero-length vector is returned.
255 * Before solving an ILP problem, we first check if simply
256 * adding the normal of the constraint to one of the known
257 * integer points in the basic set represented by "tab"
258 * yields another point inside the basic set.
260 * The caller of this function ensures that the tableau is bounded or
261 * that tab->basis and tab->n_unbounded have been set appropriately.
263 static struct isl_vec
*outside_point(struct isl_tab
*tab
, isl_int
*eq
, int up
)
266 struct isl_vec
*sample
= NULL
;
267 struct isl_tab_undo
*snap
;
275 sample
= isl_vec_alloc(ctx
, 1 + dim
);
278 isl_int_set_si(sample
->el
[0], 1);
279 isl_seq_combine(sample
->el
+ 1,
280 ctx
->one
, tab
->bmap
->sample
->el
+ 1,
281 up
? ctx
->one
: ctx
->negone
, eq
+ 1, dim
);
282 if (isl_basic_map_contains(tab
->bmap
, sample
))
284 isl_vec_free(sample
);
287 snap
= isl_tab_snap(tab
);
290 isl_seq_neg(eq
, eq
, 1 + dim
);
291 isl_int_sub_ui(eq
[0], eq
[0], 1);
293 if (isl_tab_extend_cons(tab
, 1) < 0)
295 if (isl_tab_add_ineq(tab
, eq
) < 0)
298 sample
= isl_tab_sample(tab
);
300 isl_int_add_ui(eq
[0], eq
[0], 1);
302 isl_seq_neg(eq
, eq
, 1 + dim
);
304 if (sample
&& isl_tab_rollback(tab
, snap
) < 0)
309 isl_vec_free(sample
);
313 struct isl_basic_set
*isl_basic_set_recession_cone(struct isl_basic_set
*bset
)
317 bset
= isl_basic_set_cow(bset
);
320 isl_assert(bset
->ctx
, bset
->n_div
== 0, goto error
);
322 for (i
= 0; i
< bset
->n_eq
; ++i
)
323 isl_int_set_si(bset
->eq
[i
][0], 0);
325 for (i
= 0; i
< bset
->n_ineq
; ++i
)
326 isl_int_set_si(bset
->ineq
[i
][0], 0);
328 ISL_F_CLR(bset
, ISL_BASIC_SET_NO_IMPLICIT
);
329 return isl_basic_set_implicit_equalities(bset
);
331 isl_basic_set_free(bset
);
335 __isl_give isl_set
*isl_set_recession_cone(__isl_take isl_set
*set
)
344 set
= isl_set_remove_divs(set
);
345 set
= isl_set_cow(set
);
349 for (i
= 0; i
< set
->n
; ++i
) {
350 set
->p
[i
] = isl_basic_set_recession_cone(set
->p
[i
]);
361 /* Move "sample" to a point that is one up (or down) from the original
362 * point in dimension "pos".
364 static void adjacent_point(__isl_keep isl_vec
*sample
, int pos
, int up
)
367 isl_int_add_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
369 isl_int_sub_ui(sample
->el
[1 + pos
], sample
->el
[1 + pos
], 1);
372 /* Check if any points that are adjacent to "sample" also belong to "bset".
373 * If so, add them to "hull" and return the updated hull.
375 * Before checking whether and adjacent point belongs to "bset", we first
376 * check whether it already belongs to "hull" as this test is typically
379 static __isl_give isl_basic_set
*add_adjacent_points(
380 __isl_take isl_basic_set
*hull
, __isl_take isl_vec
*sample
,
381 __isl_keep isl_basic_set
*bset
)
389 dim
= isl_basic_set_dim(hull
, isl_dim_set
);
391 for (i
= 0; i
< dim
; ++i
) {
392 for (up
= 0; up
<= 1; ++up
) {
394 isl_basic_set
*point
;
396 adjacent_point(sample
, i
, up
);
397 contains
= isl_basic_set_contains(hull
, sample
);
401 adjacent_point(sample
, i
, !up
);
404 contains
= isl_basic_set_contains(bset
, sample
);
408 point
= isl_basic_set_from_vec(
409 isl_vec_copy(sample
));
410 hull
= affine_hull(hull
, point
);
412 adjacent_point(sample
, i
, !up
);
418 isl_vec_free(sample
);
422 isl_vec_free(sample
);
423 isl_basic_set_free(hull
);
427 /* Extend an initial (under-)approximation of the affine hull of basic
428 * set represented by the tableau "tab"
429 * by looking for points that do not satisfy one of the equalities
430 * in the current approximation and adding them to that approximation
431 * until no such points can be found any more.
433 * The caller of this function ensures that "tab" is bounded or
434 * that tab->basis and tab->n_unbounded have been set appropriately.
436 * "bset" may be either NULL or the basic set represented by "tab".
437 * If "bset" is not NULL, we check for any point we find if any
438 * of its adjacent points also belong to "bset".
440 static __isl_give isl_basic_set
*extend_affine_hull(struct isl_tab
*tab
,
441 __isl_take isl_basic_set
*hull
, __isl_keep isl_basic_set
*bset
)
451 if (isl_tab_extend_cons(tab
, 2 * dim
+ 1) < 0)
454 for (i
= 0; i
< dim
; ++i
) {
455 struct isl_vec
*sample
;
456 struct isl_basic_set
*point
;
457 for (j
= 0; j
< hull
->n_eq
; ++j
) {
458 sample
= outside_point(tab
, hull
->eq
[j
], 1);
461 if (sample
->size
> 0)
463 isl_vec_free(sample
);
464 sample
= outside_point(tab
, hull
->eq
[j
], 0);
467 if (sample
->size
> 0)
469 isl_vec_free(sample
);
471 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
477 isl_tab_add_sample(tab
, isl_vec_copy(sample
)) < 0)
478 hull
= isl_basic_set_free(hull
);
480 hull
= add_adjacent_points(hull
, isl_vec_copy(sample
),
482 point
= isl_basic_set_from_vec(sample
);
483 hull
= affine_hull(hull
, point
);
490 isl_basic_set_free(hull
);
494 /* Drop all constraints in bmap that involve any of the dimensions
495 * first to first+n-1.
497 static __isl_give isl_basic_map
*isl_basic_map_drop_constraints_involving(
498 __isl_take isl_basic_map
*bmap
, unsigned first
, unsigned n
)
505 bmap
= isl_basic_map_cow(bmap
);
510 for (i
= bmap
->n_eq
- 1; i
>= 0; --i
) {
511 if (isl_seq_first_non_zero(bmap
->eq
[i
] + 1 + first
, n
) == -1)
513 isl_basic_map_drop_equality(bmap
, i
);
516 for (i
= bmap
->n_ineq
- 1; i
>= 0; --i
) {
517 if (isl_seq_first_non_zero(bmap
->ineq
[i
] + 1 + first
, n
) == -1)
519 isl_basic_map_drop_inequality(bmap
, i
);
525 /* Drop all constraints in bset that involve any of the dimensions
526 * first to first+n-1.
528 __isl_give isl_basic_set
*isl_basic_set_drop_constraints_involving(
529 __isl_take isl_basic_set
*bset
, unsigned first
, unsigned n
)
531 return isl_basic_map_drop_constraints_involving(bset
, first
, n
);
534 /* Drop all constraints in bmap that do not involve any of the dimensions
535 * first to first + n - 1 of the given type.
537 __isl_give isl_basic_map
*isl_basic_map_drop_constraints_not_involving_dims(
538 __isl_take isl_basic_map
*bmap
,
539 enum isl_dim_type type
, unsigned first
, unsigned n
)
545 return isl_basic_map_set_to_empty(bmap
);
546 bmap
= isl_basic_map_cow(bmap
);
550 dim
= isl_basic_map_dim(bmap
, type
);
551 if (first
+ n
> dim
|| first
+ n
< first
)
552 isl_die(isl_basic_map_get_ctx(bmap
), isl_error_invalid
,
553 "index out of bounds", return isl_basic_map_free(bmap
));
555 first
+= isl_basic_map_offset(bmap
, type
) - 1;
557 for (i
= bmap
->n_eq
- 1; i
>= 0; --i
) {
558 if (isl_seq_first_non_zero(bmap
->eq
[i
] + 1 + first
, n
) != -1)
560 isl_basic_map_drop_equality(bmap
, i
);
563 for (i
= bmap
->n_ineq
- 1; i
>= 0; --i
) {
564 if (isl_seq_first_non_zero(bmap
->ineq
[i
] + 1 + first
, n
) != -1)
566 isl_basic_map_drop_inequality(bmap
, i
);
572 /* Drop all constraints in bset that do not involve any of the dimensions
573 * first to first + n - 1 of the given type.
575 __isl_give isl_basic_set
*isl_basic_set_drop_constraints_not_involving_dims(
576 __isl_take isl_basic_set
*bset
,
577 enum isl_dim_type type
, unsigned first
, unsigned n
)
579 return isl_basic_map_drop_constraints_not_involving_dims(bset
,
583 /* Drop all constraints in bmap that involve any of the dimensions
584 * first to first + n - 1 of the given type.
586 __isl_give isl_basic_map
*isl_basic_map_drop_constraints_involving_dims(
587 __isl_take isl_basic_map
*bmap
,
588 enum isl_dim_type type
, unsigned first
, unsigned n
)
597 dim
= isl_basic_map_dim(bmap
, type
);
598 if (first
+ n
> dim
|| first
+ n
< first
)
599 isl_die(isl_basic_map_get_ctx(bmap
), isl_error_invalid
,
600 "index out of bounds", return isl_basic_map_free(bmap
));
602 bmap
= isl_basic_map_remove_divs_involving_dims(bmap
, type
, first
, n
);
603 first
+= isl_basic_map_offset(bmap
, type
) - 1;
604 return isl_basic_map_drop_constraints_involving(bmap
, first
, n
);
607 /* Drop all constraints in bset that involve any of the dimensions
608 * first to first + n - 1 of the given type.
610 __isl_give isl_basic_set
*isl_basic_set_drop_constraints_involving_dims(
611 __isl_take isl_basic_set
*bset
,
612 enum isl_dim_type type
, unsigned first
, unsigned n
)
614 return isl_basic_map_drop_constraints_involving_dims(bset
,
618 /* Drop all constraints in map that involve any of the dimensions
619 * first to first + n - 1 of the given type.
621 __isl_give isl_map
*isl_map_drop_constraints_involving_dims(
622 __isl_take isl_map
*map
,
623 enum isl_dim_type type
, unsigned first
, unsigned n
)
633 dim
= isl_map_dim(map
, type
);
634 if (first
+ n
> dim
|| first
+ n
< first
)
635 isl_die(isl_map_get_ctx(map
), isl_error_invalid
,
636 "index out of bounds", return isl_map_free(map
));
638 map
= isl_map_cow(map
);
642 for (i
= 0; i
< map
->n
; ++i
) {
643 map
->p
[i
] = isl_basic_map_drop_constraints_involving_dims(
644 map
->p
[i
], type
, first
, n
);
646 return isl_map_free(map
);
652 /* Drop all constraints in set that involve any of the dimensions
653 * first to first + n - 1 of the given type.
655 __isl_give isl_set
*isl_set_drop_constraints_involving_dims(
656 __isl_take isl_set
*set
,
657 enum isl_dim_type type
, unsigned first
, unsigned n
)
659 return isl_map_drop_constraints_involving_dims(set
, type
, first
, n
);
662 /* Construct an initial underapproximatino of the hull of "bset"
663 * from "sample" and any of its adjacent points that also belong to "bset".
665 static __isl_give isl_basic_set
*initialize_hull(__isl_keep isl_basic_set
*bset
,
666 __isl_take isl_vec
*sample
)
670 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
671 hull
= add_adjacent_points(hull
, sample
, bset
);
676 /* Look for all equalities satisfied by the integer points in bset,
677 * which is assumed to be bounded.
679 * The equalities are obtained by successively looking for
680 * a point that is affinely independent of the points found so far.
681 * In particular, for each equality satisfied by the points so far,
682 * we check if there is any point on a hyperplane parallel to the
683 * corresponding hyperplane shifted by at least one (in either direction).
685 static struct isl_basic_set
*uset_affine_hull_bounded(struct isl_basic_set
*bset
)
687 struct isl_vec
*sample
= NULL
;
688 struct isl_basic_set
*hull
;
689 struct isl_tab
*tab
= NULL
;
692 if (isl_basic_set_plain_is_empty(bset
))
695 dim
= isl_basic_set_n_dim(bset
);
697 if (bset
->sample
&& bset
->sample
->size
== 1 + dim
) {
698 int contains
= isl_basic_set_contains(bset
, bset
->sample
);
704 sample
= isl_vec_copy(bset
->sample
);
706 isl_vec_free(bset
->sample
);
711 tab
= isl_tab_from_basic_set(bset
, 1);
716 isl_vec_free(sample
);
717 return isl_basic_set_set_to_empty(bset
);
721 struct isl_tab_undo
*snap
;
722 snap
= isl_tab_snap(tab
);
723 sample
= isl_tab_sample(tab
);
724 if (isl_tab_rollback(tab
, snap
) < 0)
726 isl_vec_free(tab
->bmap
->sample
);
727 tab
->bmap
->sample
= isl_vec_copy(sample
);
732 if (sample
->size
== 0) {
734 isl_vec_free(sample
);
735 return isl_basic_set_set_to_empty(bset
);
738 hull
= initialize_hull(bset
, sample
);
740 hull
= extend_affine_hull(tab
, hull
, bset
);
741 isl_basic_set_free(bset
);
746 isl_vec_free(sample
);
748 isl_basic_set_free(bset
);
752 /* Given an unbounded tableau and an integer point satisfying the tableau,
753 * construct an initial affine hull containing the recession cone
754 * shifted to the given point.
756 * The unbounded directions are taken from the last rows of the basis,
757 * which is assumed to have been initialized appropriately.
759 static __isl_give isl_basic_set
*initial_hull(struct isl_tab
*tab
,
760 __isl_take isl_vec
*vec
)
764 struct isl_basic_set
*bset
= NULL
;
771 isl_assert(ctx
, vec
->size
!= 0, goto error
);
773 bset
= isl_basic_set_alloc(ctx
, 0, vec
->size
- 1, 0, vec
->size
- 1, 0);
776 dim
= isl_basic_set_n_dim(bset
) - tab
->n_unbounded
;
777 for (i
= 0; i
< dim
; ++i
) {
778 k
= isl_basic_set_alloc_equality(bset
);
781 isl_seq_cpy(bset
->eq
[k
] + 1, tab
->basis
->row
[1 + i
] + 1,
783 isl_seq_inner_product(bset
->eq
[k
] + 1, vec
->el
+1,
784 vec
->size
- 1, &bset
->eq
[k
][0]);
785 isl_int_neg(bset
->eq
[k
][0], bset
->eq
[k
][0]);
788 bset
= isl_basic_set_gauss(bset
, NULL
);
792 isl_basic_set_free(bset
);
797 /* Given a tableau of a set and a tableau of the corresponding
798 * recession cone, detect and add all equalities to the tableau.
799 * If the tableau is bounded, then we can simply keep the
800 * tableau in its state after the return from extend_affine_hull.
801 * However, if the tableau is unbounded, then
802 * isl_tab_set_initial_basis_with_cone will add some additional
803 * constraints to the tableau that have to be removed again.
804 * In this case, we therefore rollback to the state before
805 * any constraints were added and then add the equalities back in.
807 struct isl_tab
*isl_tab_detect_equalities(struct isl_tab
*tab
,
808 struct isl_tab
*tab_cone
)
811 struct isl_vec
*sample
;
812 struct isl_basic_set
*hull
= NULL
;
813 struct isl_tab_undo
*snap
;
815 if (!tab
|| !tab_cone
)
818 snap
= isl_tab_snap(tab
);
820 isl_mat_free(tab
->basis
);
823 isl_assert(tab
->mat
->ctx
, tab
->bmap
, goto error
);
824 isl_assert(tab
->mat
->ctx
, tab
->samples
, goto error
);
825 isl_assert(tab
->mat
->ctx
, tab
->samples
->n_col
== 1 + tab
->n_var
, goto error
);
826 isl_assert(tab
->mat
->ctx
, tab
->n_sample
> tab
->n_outside
, goto error
);
828 if (isl_tab_set_initial_basis_with_cone(tab
, tab_cone
) < 0)
831 sample
= isl_vec_alloc(tab
->mat
->ctx
, 1 + tab
->n_var
);
835 isl_seq_cpy(sample
->el
, tab
->samples
->row
[tab
->n_outside
], sample
->size
);
837 isl_vec_free(tab
->bmap
->sample
);
838 tab
->bmap
->sample
= isl_vec_copy(sample
);
840 if (tab
->n_unbounded
== 0)
841 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
843 hull
= initial_hull(tab
, isl_vec_copy(sample
));
845 for (j
= tab
->n_outside
+ 1; j
< tab
->n_sample
; ++j
) {
846 isl_seq_cpy(sample
->el
, tab
->samples
->row
[j
], sample
->size
);
847 hull
= affine_hull(hull
,
848 isl_basic_set_from_vec(isl_vec_copy(sample
)));
851 isl_vec_free(sample
);
853 hull
= extend_affine_hull(tab
, hull
, NULL
);
857 if (tab
->n_unbounded
== 0) {
858 isl_basic_set_free(hull
);
862 if (isl_tab_rollback(tab
, snap
) < 0)
865 if (hull
->n_eq
> tab
->n_zero
) {
866 for (j
= 0; j
< hull
->n_eq
; ++j
) {
867 isl_seq_normalize(tab
->mat
->ctx
, hull
->eq
[j
], 1 + tab
->n_var
);
868 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
873 isl_basic_set_free(hull
);
877 isl_basic_set_free(hull
);
882 /* Compute the affine hull of "bset", where "cone" is the recession cone
885 * We first compute a unimodular transformation that puts the unbounded
886 * directions in the last dimensions. In particular, we take a transformation
887 * that maps all equalities to equalities (in HNF) on the first dimensions.
888 * Let x be the original dimensions and y the transformed, with y_1 bounded
891 * [ y_1 ] [ y_1 ] [ Q_1 ]
892 * x = U [ y_2 ] [ y_2 ] = [ Q_2 ] x
894 * Let's call the input basic set S. We compute S' = preimage(S, U)
895 * and drop the final dimensions including any constraints involving them.
896 * This results in set S''.
897 * Then we compute the affine hull A'' of S''.
898 * Let F y_1 >= g be the constraint system of A''. In the transformed
899 * space the y_2 are unbounded, so we can add them back without any constraints,
903 * [ F 0 ] [ y_2 ] >= g
906 * [ F 0 ] [ Q_2 ] x >= g
910 * The affine hull in the original space is then obtained as
911 * A = preimage(A'', Q_1).
913 static struct isl_basic_set
*affine_hull_with_cone(struct isl_basic_set
*bset
,
914 struct isl_basic_set
*cone
)
918 struct isl_basic_set
*hull
;
919 struct isl_mat
*M
, *U
, *Q
;
924 total
= isl_basic_set_total_dim(cone
);
925 cone_dim
= total
- cone
->n_eq
;
927 M
= isl_mat_sub_alloc6(bset
->ctx
, cone
->eq
, 0, cone
->n_eq
, 1, total
);
928 M
= isl_mat_left_hermite(M
, 0, &U
, &Q
);
933 U
= isl_mat_lin_to_aff(U
);
934 bset
= isl_basic_set_preimage(bset
, isl_mat_copy(U
));
936 bset
= isl_basic_set_drop_constraints_involving(bset
, total
- cone_dim
,
938 bset
= isl_basic_set_drop_dims(bset
, total
- cone_dim
, cone_dim
);
940 Q
= isl_mat_lin_to_aff(Q
);
941 Q
= isl_mat_drop_rows(Q
, 1 + total
- cone_dim
, cone_dim
);
943 if (bset
&& bset
->sample
&& bset
->sample
->size
== 1 + total
)
944 bset
->sample
= isl_mat_vec_product(isl_mat_copy(Q
), bset
->sample
);
946 hull
= uset_affine_hull_bounded(bset
);
952 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
953 U
= isl_mat_drop_cols(U
, 1 + total
- cone_dim
, cone_dim
);
954 if (sample
&& sample
->size
> 0)
955 sample
= isl_mat_vec_product(U
, sample
);
958 hull
= isl_basic_set_preimage(hull
, Q
);
960 isl_vec_free(hull
->sample
);
961 hull
->sample
= sample
;
963 isl_vec_free(sample
);
966 isl_basic_set_free(cone
);
970 isl_basic_set_free(bset
);
971 isl_basic_set_free(cone
);
975 /* Look for all equalities satisfied by the integer points in bset,
976 * which is assumed not to have any explicit equalities.
978 * The equalities are obtained by successively looking for
979 * a point that is affinely independent of the points found so far.
980 * In particular, for each equality satisfied by the points so far,
981 * we check if there is any point on a hyperplane parallel to the
982 * corresponding hyperplane shifted by at least one (in either direction).
984 * Before looking for any outside points, we first compute the recession
985 * cone. The directions of this recession cone will always be part
986 * of the affine hull, so there is no need for looking for any points
987 * in these directions.
988 * In particular, if the recession cone is full-dimensional, then
989 * the affine hull is simply the whole universe.
991 static struct isl_basic_set
*uset_affine_hull(struct isl_basic_set
*bset
)
993 struct isl_basic_set
*cone
;
995 if (isl_basic_set_plain_is_empty(bset
))
998 cone
= isl_basic_set_recession_cone(isl_basic_set_copy(bset
));
1001 if (cone
->n_eq
== 0) {
1002 struct isl_basic_set
*hull
;
1003 isl_basic_set_free(cone
);
1004 hull
= isl_basic_set_universe_like(bset
);
1005 isl_basic_set_free(bset
);
1009 if (cone
->n_eq
< isl_basic_set_total_dim(cone
))
1010 return affine_hull_with_cone(bset
, cone
);
1012 isl_basic_set_free(cone
);
1013 return uset_affine_hull_bounded(bset
);
1015 isl_basic_set_free(bset
);
1019 /* Look for all equalities satisfied by the integer points in bmap
1020 * that are independent of the equalities already explicitly available
1023 * We first remove all equalities already explicitly available,
1024 * then look for additional equalities in the reduced space
1025 * and then transform the result to the original space.
1026 * The original equalities are _not_ added to this set. This is
1027 * the responsibility of the calling function.
1028 * The resulting basic set has all meaning about the dimensions removed.
1029 * In particular, dimensions that correspond to existential variables
1030 * in bmap and that are found to be fixed are not removed.
1032 static struct isl_basic_set
*equalities_in_underlying_set(
1033 struct isl_basic_map
*bmap
)
1035 struct isl_mat
*T1
= NULL
;
1036 struct isl_mat
*T2
= NULL
;
1037 struct isl_basic_set
*bset
= NULL
;
1038 struct isl_basic_set
*hull
= NULL
;
1040 bset
= isl_basic_map_underlying_set(bmap
);
1044 bset
= isl_basic_set_remove_equalities(bset
, &T1
, &T2
);
1048 hull
= uset_affine_hull(bset
);
1056 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
1057 if (sample
&& sample
->size
> 0)
1058 sample
= isl_mat_vec_product(T1
, sample
);
1061 hull
= isl_basic_set_preimage(hull
, T2
);
1063 isl_vec_free(hull
->sample
);
1064 hull
->sample
= sample
;
1066 isl_vec_free(sample
);
1073 isl_basic_set_free(bset
);
1074 isl_basic_set_free(hull
);
1078 /* Detect and make explicit all equalities satisfied by the (integer)
1081 struct isl_basic_map
*isl_basic_map_detect_equalities(
1082 struct isl_basic_map
*bmap
)
1085 struct isl_basic_set
*hull
= NULL
;
1089 if (bmap
->n_ineq
== 0)
1091 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
1093 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_ALL_EQUALITIES
))
1095 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_RATIONAL
))
1096 return isl_basic_map_implicit_equalities(bmap
);
1098 hull
= equalities_in_underlying_set(isl_basic_map_copy(bmap
));
1101 if (ISL_F_ISSET(hull
, ISL_BASIC_SET_EMPTY
)) {
1102 isl_basic_set_free(hull
);
1103 return isl_basic_map_set_to_empty(bmap
);
1105 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
), 0,
1107 for (i
= 0; i
< hull
->n_eq
; ++i
) {
1108 j
= isl_basic_map_alloc_equality(bmap
);
1111 isl_seq_cpy(bmap
->eq
[j
], hull
->eq
[i
],
1112 1 + isl_basic_set_total_dim(hull
));
1114 isl_vec_free(bmap
->sample
);
1115 bmap
->sample
= isl_vec_copy(hull
->sample
);
1116 isl_basic_set_free(hull
);
1117 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
| ISL_BASIC_MAP_ALL_EQUALITIES
);
1118 bmap
= isl_basic_map_simplify(bmap
);
1119 return isl_basic_map_finalize(bmap
);
1121 isl_basic_set_free(hull
);
1122 isl_basic_map_free(bmap
);
1126 __isl_give isl_basic_set
*isl_basic_set_detect_equalities(
1127 __isl_take isl_basic_set
*bset
)
1129 return (isl_basic_set
*)
1130 isl_basic_map_detect_equalities((isl_basic_map
*)bset
);
1133 __isl_give isl_map
*isl_map_detect_equalities(__isl_take isl_map
*map
)
1135 return isl_map_inline_foreach_basic_map(map
,
1136 &isl_basic_map_detect_equalities
);
1139 __isl_give isl_set
*isl_set_detect_equalities(__isl_take isl_set
*set
)
1141 return (isl_set
*)isl_map_detect_equalities((isl_map
*)set
);
1144 /* After computing the rational affine hull (by detecting the implicit
1145 * equalities), we compute the additional equalities satisfied by
1146 * the integer points (if any) and add the original equalities back in.
1148 struct isl_basic_map
*isl_basic_map_affine_hull(struct isl_basic_map
*bmap
)
1150 bmap
= isl_basic_map_detect_equalities(bmap
);
1151 bmap
= isl_basic_map_cow(bmap
);
1153 isl_basic_map_free_inequality(bmap
, bmap
->n_ineq
);
1154 bmap
= isl_basic_map_finalize(bmap
);
1158 struct isl_basic_set
*isl_basic_set_affine_hull(struct isl_basic_set
*bset
)
1160 return (struct isl_basic_set
*)
1161 isl_basic_map_affine_hull((struct isl_basic_map
*)bset
);
1164 /* Given a rational affine matrix "M", add stride constraints to "bmap"
1169 * is an integer vector. The variables x include all the variables
1170 * of "bmap" except the unknown divs.
1172 * If d is the common denominator of M, then we need to impose that
1178 * exists alpha : d M(x) = d alpha
1180 * This function is similar to add_strides in isl_morph.c
1182 static __isl_give isl_basic_map
*add_strides(__isl_take isl_basic_map
*bmap
,
1183 __isl_keep isl_mat
*M
, int n_known
)
1188 if (isl_int_is_one(M
->row
[0][0]))
1191 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
),
1192 M
->n_row
- 1, M
->n_row
- 1, 0);
1195 for (i
= 1; i
< M
->n_row
; ++i
) {
1196 isl_seq_gcd(M
->row
[i
], M
->n_col
, &gcd
);
1197 if (isl_int_is_divisible_by(gcd
, M
->row
[0][0]))
1199 div
= isl_basic_map_alloc_div(bmap
);
1202 isl_int_set_si(bmap
->div
[div
][0], 0);
1203 k
= isl_basic_map_alloc_equality(bmap
);
1206 isl_seq_cpy(bmap
->eq
[k
], M
->row
[i
], M
->n_col
);
1207 isl_seq_clr(bmap
->eq
[k
] + M
->n_col
, bmap
->n_div
- n_known
);
1208 isl_int_set(bmap
->eq
[k
][M
->n_col
- n_known
+ div
],
1216 isl_basic_map_free(bmap
);
1220 /* If there are any equalities that involve (multiple) unknown divs,
1221 * then extract the stride information encoded by those equalities
1222 * and make it explicitly available in "bmap".
1224 * We first sort the divs so that the unknown divs appear last and
1225 * then we count how many equalities involve these divs.
1227 * Let these equalities be of the form
1231 * where y represents the unknown divs and x the remaining variables.
1232 * Let [H 0] be the Hermite Normal Form of B, i.e.,
1236 * Then x is a solution of the equalities iff
1238 * H^-1 A(x) (= - [I 0] Q y)
1240 * is an integer vector. Let d be the common denominator of H^-1.
1243 * d H^-1 A(x) = d alpha
1245 * in add_strides, with alpha fresh existentially quantified variables.
1247 static __isl_give isl_basic_map
*isl_basic_map_make_strides_explicit(
1248 __isl_take isl_basic_map
*bmap
)
1257 known
= isl_basic_map_divs_known(bmap
);
1259 return isl_basic_map_free(bmap
);
1262 bmap
= isl_basic_map_sort_divs(bmap
);
1263 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1267 for (n_known
= 0; n_known
< bmap
->n_div
; ++n_known
)
1268 if (isl_int_is_zero(bmap
->div
[n_known
][0]))
1270 ctx
= isl_basic_map_get_ctx(bmap
);
1271 total
= isl_space_dim(bmap
->dim
, isl_dim_all
);
1272 for (n
= 0; n
< bmap
->n_eq
; ++n
)
1273 if (isl_seq_first_non_zero(bmap
->eq
[n
] + 1 + total
+ n_known
,
1274 bmap
->n_div
- n_known
) == -1)
1278 B
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 0, 1 + total
+ n_known
);
1279 n_col
= bmap
->n_div
- n_known
;
1280 A
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 1 + total
+ n_known
, n_col
);
1281 A
= isl_mat_left_hermite(A
, 0, NULL
, NULL
);
1282 A
= isl_mat_drop_cols(A
, n
, n_col
- n
);
1283 A
= isl_mat_lin_to_aff(A
);
1284 A
= isl_mat_right_inverse(A
);
1285 B
= isl_mat_insert_zero_rows(B
, 0, 1);
1286 B
= isl_mat_set_element_si(B
, 0, 0, 1);
1287 M
= isl_mat_product(A
, B
);
1289 return isl_basic_map_free(bmap
);
1290 bmap
= add_strides(bmap
, M
, n_known
);
1291 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1297 /* Compute the affine hull of each basic map in "map" separately
1298 * and make all stride information explicit so that we can remove
1299 * all unknown divs without losing this information.
1300 * The result is also guaranteed to be gaussed.
1302 * In simple cases where a div is determined by an equality,
1303 * calling isl_basic_map_gauss is enough to make the stride information
1304 * explicit, as it will derive an explicit representation for the div
1305 * from the equality. If, however, the stride information
1306 * is encoded through multiple unknown divs then we need to make
1307 * some extra effort in isl_basic_map_make_strides_explicit.
1309 static __isl_give isl_map
*isl_map_local_affine_hull(__isl_take isl_map
*map
)
1313 map
= isl_map_cow(map
);
1317 for (i
= 0; i
< map
->n
; ++i
) {
1318 map
->p
[i
] = isl_basic_map_affine_hull(map
->p
[i
]);
1319 map
->p
[i
] = isl_basic_map_gauss(map
->p
[i
], NULL
);
1320 map
->p
[i
] = isl_basic_map_make_strides_explicit(map
->p
[i
]);
1322 return isl_map_free(map
);
1328 static __isl_give isl_set
*isl_set_local_affine_hull(__isl_take isl_set
*set
)
1330 return isl_map_local_affine_hull(set
);
1333 /* Compute the affine hull of "map".
1335 * We first compute the affine hull of each basic map separately.
1336 * Then we align the divs and recompute the affine hulls of the basic
1337 * maps since some of them may now have extra divs.
1338 * In order to avoid performing parametric integer programming to
1339 * compute explicit expressions for the divs, possible leading to
1340 * an explosion in the number of basic maps, we first drop all unknown
1341 * divs before aligning the divs. Note that isl_map_local_affine_hull tries
1342 * to make sure that all stride information is explicitly available
1343 * in terms of known divs. This involves calling isl_basic_set_gauss,
1344 * which is also needed because affine_hull assumes its input has been gaussed,
1345 * while isl_map_affine_hull may be called on input that has not been gaussed,
1346 * in particular from initial_facet_constraint.
1347 * Similarly, align_divs may reorder some divs so that we need to
1348 * gauss the result again.
1349 * Finally, we combine the individual affine hulls into a single
1352 __isl_give isl_basic_map
*isl_map_affine_hull(__isl_take isl_map
*map
)
1354 struct isl_basic_map
*model
= NULL
;
1355 struct isl_basic_map
*hull
= NULL
;
1356 struct isl_set
*set
;
1357 isl_basic_set
*bset
;
1359 map
= isl_map_detect_equalities(map
);
1360 map
= isl_map_local_affine_hull(map
);
1361 map
= isl_map_remove_empty_parts(map
);
1362 map
= isl_map_remove_unknown_divs(map
);
1363 map
= isl_map_align_divs(map
);
1369 hull
= isl_basic_map_empty_like_map(map
);
1374 model
= isl_basic_map_copy(map
->p
[0]);
1375 set
= isl_map_underlying_set(map
);
1376 set
= isl_set_cow(set
);
1377 set
= isl_set_local_affine_hull(set
);
1382 set
->p
[0] = affine_hull(set
->p
[0], set
->p
[--set
->n
]);
1384 bset
= isl_basic_set_copy(set
->p
[0]);
1385 hull
= isl_basic_map_overlying_set(bset
, model
);
1387 hull
= isl_basic_map_simplify(hull
);
1388 return isl_basic_map_finalize(hull
);
1390 isl_basic_map_free(model
);
1395 struct isl_basic_set
*isl_set_affine_hull(struct isl_set
*set
)
1397 return (struct isl_basic_set
*)
1398 isl_map_affine_hull((struct isl_map
*)set
);