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 tab
= isl_tab_add_sample(tab
, isl_vec_copy(sample
));
481 hull
= add_adjacent_points(hull
, isl_vec_copy(sample
),
483 point
= isl_basic_set_from_vec(sample
);
484 hull
= affine_hull(hull
, point
);
491 isl_basic_set_free(hull
);
495 /* Drop all constraints in bmap that involve any of the dimensions
496 * first to first+n-1.
498 static __isl_give isl_basic_map
*isl_basic_map_drop_constraints_involving(
499 __isl_take isl_basic_map
*bmap
, unsigned first
, unsigned n
)
506 bmap
= isl_basic_map_cow(bmap
);
511 for (i
= bmap
->n_eq
- 1; i
>= 0; --i
) {
512 if (isl_seq_first_non_zero(bmap
->eq
[i
] + 1 + first
, n
) == -1)
514 isl_basic_map_drop_equality(bmap
, i
);
517 for (i
= bmap
->n_ineq
- 1; i
>= 0; --i
) {
518 if (isl_seq_first_non_zero(bmap
->ineq
[i
] + 1 + first
, n
) == -1)
520 isl_basic_map_drop_inequality(bmap
, i
);
526 /* Drop all constraints in bset that involve any of the dimensions
527 * first to first+n-1.
529 __isl_give isl_basic_set
*isl_basic_set_drop_constraints_involving(
530 __isl_take isl_basic_set
*bset
, unsigned first
, unsigned n
)
532 return isl_basic_map_drop_constraints_involving(bset
, first
, n
);
535 /* Drop all constraints in bmap that do not involve any of the dimensions
536 * first to first + n - 1 of the given type.
538 __isl_give isl_basic_map
*isl_basic_map_drop_constraints_not_involving_dims(
539 __isl_take isl_basic_map
*bmap
,
540 enum isl_dim_type type
, unsigned first
, unsigned n
)
546 return isl_basic_map_set_to_empty(bmap
);
547 bmap
= isl_basic_map_cow(bmap
);
551 dim
= isl_basic_map_dim(bmap
, type
);
552 if (first
+ n
> dim
|| first
+ n
< first
)
553 isl_die(isl_basic_map_get_ctx(bmap
), isl_error_invalid
,
554 "index out of bounds", return isl_basic_map_free(bmap
));
556 first
+= isl_basic_map_offset(bmap
, type
) - 1;
558 for (i
= bmap
->n_eq
- 1; i
>= 0; --i
) {
559 if (isl_seq_first_non_zero(bmap
->eq
[i
] + 1 + first
, n
) != -1)
561 isl_basic_map_drop_equality(bmap
, i
);
564 for (i
= bmap
->n_ineq
- 1; i
>= 0; --i
) {
565 if (isl_seq_first_non_zero(bmap
->ineq
[i
] + 1 + first
, n
) != -1)
567 isl_basic_map_drop_inequality(bmap
, i
);
573 /* Drop all constraints in bset that do not involve any of the dimensions
574 * first to first + n - 1 of the given type.
576 __isl_give isl_basic_set
*isl_basic_set_drop_constraints_not_involving_dims(
577 __isl_take isl_basic_set
*bset
,
578 enum isl_dim_type type
, unsigned first
, unsigned n
)
580 return isl_basic_map_drop_constraints_not_involving_dims(bset
,
584 /* Drop all constraints in bmap that involve any of the dimensions
585 * first to first + n - 1 of the given type.
587 __isl_give isl_basic_map
*isl_basic_map_drop_constraints_involving_dims(
588 __isl_take isl_basic_map
*bmap
,
589 enum isl_dim_type type
, unsigned first
, unsigned n
)
598 dim
= isl_basic_map_dim(bmap
, type
);
599 if (first
+ n
> dim
|| first
+ n
< first
)
600 isl_die(isl_basic_map_get_ctx(bmap
), isl_error_invalid
,
601 "index out of bounds", return isl_basic_map_free(bmap
));
603 bmap
= isl_basic_map_remove_divs_involving_dims(bmap
, type
, first
, n
);
604 first
+= isl_basic_map_offset(bmap
, type
) - 1;
605 return isl_basic_map_drop_constraints_involving(bmap
, first
, n
);
608 /* Drop all constraints in bset that involve any of the dimensions
609 * first to first + n - 1 of the given type.
611 __isl_give isl_basic_set
*isl_basic_set_drop_constraints_involving_dims(
612 __isl_take isl_basic_set
*bset
,
613 enum isl_dim_type type
, unsigned first
, unsigned n
)
615 return isl_basic_map_drop_constraints_involving_dims(bset
,
619 /* Drop all constraints in map that involve any of the dimensions
620 * first to first + n - 1 of the given type.
622 __isl_give isl_map
*isl_map_drop_constraints_involving_dims(
623 __isl_take isl_map
*map
,
624 enum isl_dim_type type
, unsigned first
, unsigned n
)
634 dim
= isl_map_dim(map
, type
);
635 if (first
+ n
> dim
|| first
+ n
< first
)
636 isl_die(isl_map_get_ctx(map
), isl_error_invalid
,
637 "index out of bounds", return isl_map_free(map
));
639 map
= isl_map_cow(map
);
643 for (i
= 0; i
< map
->n
; ++i
) {
644 map
->p
[i
] = isl_basic_map_drop_constraints_involving_dims(
645 map
->p
[i
], type
, first
, n
);
647 return isl_map_free(map
);
653 /* Drop all constraints in set that involve any of the dimensions
654 * first to first + n - 1 of the given type.
656 __isl_give isl_set
*isl_set_drop_constraints_involving_dims(
657 __isl_take isl_set
*set
,
658 enum isl_dim_type type
, unsigned first
, unsigned n
)
660 return isl_map_drop_constraints_involving_dims(set
, type
, first
, n
);
663 /* Construct an initial underapproximatino of the hull of "bset"
664 * from "sample" and any of its adjacent points that also belong to "bset".
666 static __isl_give isl_basic_set
*initialize_hull(__isl_keep isl_basic_set
*bset
,
667 __isl_take isl_vec
*sample
)
671 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
672 hull
= add_adjacent_points(hull
, sample
, bset
);
677 /* Look for all equalities satisfied by the integer points in bset,
678 * which is assumed to be bounded.
680 * The equalities are obtained by successively looking for
681 * a point that is affinely independent of the points found so far.
682 * In particular, for each equality satisfied by the points so far,
683 * we check if there is any point on a hyperplane parallel to the
684 * corresponding hyperplane shifted by at least one (in either direction).
686 static struct isl_basic_set
*uset_affine_hull_bounded(struct isl_basic_set
*bset
)
688 struct isl_vec
*sample
= NULL
;
689 struct isl_basic_set
*hull
;
690 struct isl_tab
*tab
= NULL
;
693 if (isl_basic_set_plain_is_empty(bset
))
696 dim
= isl_basic_set_n_dim(bset
);
698 if (bset
->sample
&& bset
->sample
->size
== 1 + dim
) {
699 int contains
= isl_basic_set_contains(bset
, bset
->sample
);
705 sample
= isl_vec_copy(bset
->sample
);
707 isl_vec_free(bset
->sample
);
712 tab
= isl_tab_from_basic_set(bset
, 1);
717 isl_vec_free(sample
);
718 return isl_basic_set_set_to_empty(bset
);
722 struct isl_tab_undo
*snap
;
723 snap
= isl_tab_snap(tab
);
724 sample
= isl_tab_sample(tab
);
725 if (isl_tab_rollback(tab
, snap
) < 0)
727 isl_vec_free(tab
->bmap
->sample
);
728 tab
->bmap
->sample
= isl_vec_copy(sample
);
733 if (sample
->size
== 0) {
735 isl_vec_free(sample
);
736 return isl_basic_set_set_to_empty(bset
);
739 hull
= initialize_hull(bset
, sample
);
741 hull
= extend_affine_hull(tab
, hull
, bset
);
742 isl_basic_set_free(bset
);
747 isl_vec_free(sample
);
749 isl_basic_set_free(bset
);
753 /* Given an unbounded tableau and an integer point satisfying the tableau,
754 * construct an initial affine hull containing the recession cone
755 * shifted to the given point.
757 * The unbounded directions are taken from the last rows of the basis,
758 * which is assumed to have been initialized appropriately.
760 static __isl_give isl_basic_set
*initial_hull(struct isl_tab
*tab
,
761 __isl_take isl_vec
*vec
)
765 struct isl_basic_set
*bset
= NULL
;
772 isl_assert(ctx
, vec
->size
!= 0, goto error
);
774 bset
= isl_basic_set_alloc(ctx
, 0, vec
->size
- 1, 0, vec
->size
- 1, 0);
777 dim
= isl_basic_set_n_dim(bset
) - tab
->n_unbounded
;
778 for (i
= 0; i
< dim
; ++i
) {
779 k
= isl_basic_set_alloc_equality(bset
);
782 isl_seq_cpy(bset
->eq
[k
] + 1, tab
->basis
->row
[1 + i
] + 1,
784 isl_seq_inner_product(bset
->eq
[k
] + 1, vec
->el
+1,
785 vec
->size
- 1, &bset
->eq
[k
][0]);
786 isl_int_neg(bset
->eq
[k
][0], bset
->eq
[k
][0]);
789 bset
= isl_basic_set_gauss(bset
, NULL
);
793 isl_basic_set_free(bset
);
798 /* Given a tableau of a set and a tableau of the corresponding
799 * recession cone, detect and add all equalities to the tableau.
800 * If the tableau is bounded, then we can simply keep the
801 * tableau in its state after the return from extend_affine_hull.
802 * However, if the tableau is unbounded, then
803 * isl_tab_set_initial_basis_with_cone will add some additional
804 * constraints to the tableau that have to be removed again.
805 * In this case, we therefore rollback to the state before
806 * any constraints were added and then add the equalities back in.
808 struct isl_tab
*isl_tab_detect_equalities(struct isl_tab
*tab
,
809 struct isl_tab
*tab_cone
)
812 struct isl_vec
*sample
;
813 struct isl_basic_set
*hull
= NULL
;
814 struct isl_tab_undo
*snap
;
816 if (!tab
|| !tab_cone
)
819 snap
= isl_tab_snap(tab
);
821 isl_mat_free(tab
->basis
);
824 isl_assert(tab
->mat
->ctx
, tab
->bmap
, goto error
);
825 isl_assert(tab
->mat
->ctx
, tab
->samples
, goto error
);
826 isl_assert(tab
->mat
->ctx
, tab
->samples
->n_col
== 1 + tab
->n_var
, goto error
);
827 isl_assert(tab
->mat
->ctx
, tab
->n_sample
> tab
->n_outside
, goto error
);
829 if (isl_tab_set_initial_basis_with_cone(tab
, tab_cone
) < 0)
832 sample
= isl_vec_alloc(tab
->mat
->ctx
, 1 + tab
->n_var
);
836 isl_seq_cpy(sample
->el
, tab
->samples
->row
[tab
->n_outside
], sample
->size
);
838 isl_vec_free(tab
->bmap
->sample
);
839 tab
->bmap
->sample
= isl_vec_copy(sample
);
841 if (tab
->n_unbounded
== 0)
842 hull
= isl_basic_set_from_vec(isl_vec_copy(sample
));
844 hull
= initial_hull(tab
, isl_vec_copy(sample
));
846 for (j
= tab
->n_outside
+ 1; j
< tab
->n_sample
; ++j
) {
847 isl_seq_cpy(sample
->el
, tab
->samples
->row
[j
], sample
->size
);
848 hull
= affine_hull(hull
,
849 isl_basic_set_from_vec(isl_vec_copy(sample
)));
852 isl_vec_free(sample
);
854 hull
= extend_affine_hull(tab
, hull
, NULL
);
858 if (tab
->n_unbounded
== 0) {
859 isl_basic_set_free(hull
);
863 if (isl_tab_rollback(tab
, snap
) < 0)
866 if (hull
->n_eq
> tab
->n_zero
) {
867 for (j
= 0; j
< hull
->n_eq
; ++j
) {
868 isl_seq_normalize(tab
->mat
->ctx
, hull
->eq
[j
], 1 + tab
->n_var
);
869 if (isl_tab_add_eq(tab
, hull
->eq
[j
]) < 0)
874 isl_basic_set_free(hull
);
878 isl_basic_set_free(hull
);
883 /* Compute the affine hull of "bset", where "cone" is the recession cone
886 * We first compute a unimodular transformation that puts the unbounded
887 * directions in the last dimensions. In particular, we take a transformation
888 * that maps all equalities to equalities (in HNF) on the first dimensions.
889 * Let x be the original dimensions and y the transformed, with y_1 bounded
892 * [ y_1 ] [ y_1 ] [ Q_1 ]
893 * x = U [ y_2 ] [ y_2 ] = [ Q_2 ] x
895 * Let's call the input basic set S. We compute S' = preimage(S, U)
896 * and drop the final dimensions including any constraints involving them.
897 * This results in set S''.
898 * Then we compute the affine hull A'' of S''.
899 * Let F y_1 >= g be the constraint system of A''. In the transformed
900 * space the y_2 are unbounded, so we can add them back without any constraints,
904 * [ F 0 ] [ y_2 ] >= g
907 * [ F 0 ] [ Q_2 ] x >= g
911 * The affine hull in the original space is then obtained as
912 * A = preimage(A'', Q_1).
914 static struct isl_basic_set
*affine_hull_with_cone(struct isl_basic_set
*bset
,
915 struct isl_basic_set
*cone
)
919 struct isl_basic_set
*hull
;
920 struct isl_mat
*M
, *U
, *Q
;
925 total
= isl_basic_set_total_dim(cone
);
926 cone_dim
= total
- cone
->n_eq
;
928 M
= isl_mat_sub_alloc6(bset
->ctx
, cone
->eq
, 0, cone
->n_eq
, 1, total
);
929 M
= isl_mat_left_hermite(M
, 0, &U
, &Q
);
934 U
= isl_mat_lin_to_aff(U
);
935 bset
= isl_basic_set_preimage(bset
, isl_mat_copy(U
));
937 bset
= isl_basic_set_drop_constraints_involving(bset
, total
- cone_dim
,
939 bset
= isl_basic_set_drop_dims(bset
, total
- cone_dim
, cone_dim
);
941 Q
= isl_mat_lin_to_aff(Q
);
942 Q
= isl_mat_drop_rows(Q
, 1 + total
- cone_dim
, cone_dim
);
944 if (bset
&& bset
->sample
&& bset
->sample
->size
== 1 + total
)
945 bset
->sample
= isl_mat_vec_product(isl_mat_copy(Q
), bset
->sample
);
947 hull
= uset_affine_hull_bounded(bset
);
953 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
954 U
= isl_mat_drop_cols(U
, 1 + total
- cone_dim
, cone_dim
);
955 if (sample
&& sample
->size
> 0)
956 sample
= isl_mat_vec_product(U
, sample
);
959 hull
= isl_basic_set_preimage(hull
, Q
);
961 isl_vec_free(hull
->sample
);
962 hull
->sample
= sample
;
964 isl_vec_free(sample
);
967 isl_basic_set_free(cone
);
971 isl_basic_set_free(bset
);
972 isl_basic_set_free(cone
);
976 /* Look for all equalities satisfied by the integer points in bset,
977 * which is assumed not to have any explicit equalities.
979 * The equalities are obtained by successively looking for
980 * a point that is affinely independent of the points found so far.
981 * In particular, for each equality satisfied by the points so far,
982 * we check if there is any point on a hyperplane parallel to the
983 * corresponding hyperplane shifted by at least one (in either direction).
985 * Before looking for any outside points, we first compute the recession
986 * cone. The directions of this recession cone will always be part
987 * of the affine hull, so there is no need for looking for any points
988 * in these directions.
989 * In particular, if the recession cone is full-dimensional, then
990 * the affine hull is simply the whole universe.
992 static struct isl_basic_set
*uset_affine_hull(struct isl_basic_set
*bset
)
994 struct isl_basic_set
*cone
;
996 if (isl_basic_set_plain_is_empty(bset
))
999 cone
= isl_basic_set_recession_cone(isl_basic_set_copy(bset
));
1002 if (cone
->n_eq
== 0) {
1003 struct isl_basic_set
*hull
;
1004 isl_basic_set_free(cone
);
1005 hull
= isl_basic_set_universe_like(bset
);
1006 isl_basic_set_free(bset
);
1010 if (cone
->n_eq
< isl_basic_set_total_dim(cone
))
1011 return affine_hull_with_cone(bset
, cone
);
1013 isl_basic_set_free(cone
);
1014 return uset_affine_hull_bounded(bset
);
1016 isl_basic_set_free(bset
);
1020 /* Look for all equalities satisfied by the integer points in bmap
1021 * that are independent of the equalities already explicitly available
1024 * We first remove all equalities already explicitly available,
1025 * then look for additional equalities in the reduced space
1026 * and then transform the result to the original space.
1027 * The original equalities are _not_ added to this set. This is
1028 * the responsibility of the calling function.
1029 * The resulting basic set has all meaning about the dimensions removed.
1030 * In particular, dimensions that correspond to existential variables
1031 * in bmap and that are found to be fixed are not removed.
1033 static struct isl_basic_set
*equalities_in_underlying_set(
1034 struct isl_basic_map
*bmap
)
1036 struct isl_mat
*T1
= NULL
;
1037 struct isl_mat
*T2
= NULL
;
1038 struct isl_basic_set
*bset
= NULL
;
1039 struct isl_basic_set
*hull
= NULL
;
1041 bset
= isl_basic_map_underlying_set(bmap
);
1045 bset
= isl_basic_set_remove_equalities(bset
, &T1
, &T2
);
1049 hull
= uset_affine_hull(bset
);
1057 struct isl_vec
*sample
= isl_vec_copy(hull
->sample
);
1058 if (sample
&& sample
->size
> 0)
1059 sample
= isl_mat_vec_product(T1
, sample
);
1062 hull
= isl_basic_set_preimage(hull
, T2
);
1064 isl_vec_free(hull
->sample
);
1065 hull
->sample
= sample
;
1067 isl_vec_free(sample
);
1074 isl_basic_set_free(bset
);
1075 isl_basic_set_free(hull
);
1079 /* Detect and make explicit all equalities satisfied by the (integer)
1082 struct isl_basic_map
*isl_basic_map_detect_equalities(
1083 struct isl_basic_map
*bmap
)
1086 struct isl_basic_set
*hull
= NULL
;
1090 if (bmap
->n_ineq
== 0)
1092 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_EMPTY
))
1094 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_ALL_EQUALITIES
))
1096 if (ISL_F_ISSET(bmap
, ISL_BASIC_MAP_RATIONAL
))
1097 return isl_basic_map_implicit_equalities(bmap
);
1099 hull
= equalities_in_underlying_set(isl_basic_map_copy(bmap
));
1102 if (ISL_F_ISSET(hull
, ISL_BASIC_SET_EMPTY
)) {
1103 isl_basic_set_free(hull
);
1104 return isl_basic_map_set_to_empty(bmap
);
1106 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
), 0,
1108 for (i
= 0; i
< hull
->n_eq
; ++i
) {
1109 j
= isl_basic_map_alloc_equality(bmap
);
1112 isl_seq_cpy(bmap
->eq
[j
], hull
->eq
[i
],
1113 1 + isl_basic_set_total_dim(hull
));
1115 isl_vec_free(bmap
->sample
);
1116 bmap
->sample
= isl_vec_copy(hull
->sample
);
1117 isl_basic_set_free(hull
);
1118 ISL_F_SET(bmap
, ISL_BASIC_MAP_NO_IMPLICIT
| ISL_BASIC_MAP_ALL_EQUALITIES
);
1119 bmap
= isl_basic_map_simplify(bmap
);
1120 return isl_basic_map_finalize(bmap
);
1122 isl_basic_set_free(hull
);
1123 isl_basic_map_free(bmap
);
1127 __isl_give isl_basic_set
*isl_basic_set_detect_equalities(
1128 __isl_take isl_basic_set
*bset
)
1130 return (isl_basic_set
*)
1131 isl_basic_map_detect_equalities((isl_basic_map
*)bset
);
1134 __isl_give isl_map
*isl_map_detect_equalities(__isl_take isl_map
*map
)
1136 return isl_map_inline_foreach_basic_map(map
,
1137 &isl_basic_map_detect_equalities
);
1140 __isl_give isl_set
*isl_set_detect_equalities(__isl_take isl_set
*set
)
1142 return (isl_set
*)isl_map_detect_equalities((isl_map
*)set
);
1145 /* After computing the rational affine hull (by detecting the implicit
1146 * equalities), we compute the additional equalities satisfied by
1147 * the integer points (if any) and add the original equalities back in.
1149 struct isl_basic_map
*isl_basic_map_affine_hull(struct isl_basic_map
*bmap
)
1151 bmap
= isl_basic_map_detect_equalities(bmap
);
1152 bmap
= isl_basic_map_cow(bmap
);
1154 isl_basic_map_free_inequality(bmap
, bmap
->n_ineq
);
1155 bmap
= isl_basic_map_finalize(bmap
);
1159 struct isl_basic_set
*isl_basic_set_affine_hull(struct isl_basic_set
*bset
)
1161 return (struct isl_basic_set
*)
1162 isl_basic_map_affine_hull((struct isl_basic_map
*)bset
);
1165 /* Given a rational affine matrix "M", add stride constraints to "bmap"
1170 * is an integer vector. The variables x include all the variables
1171 * of "bmap" except the unknown divs.
1173 * If d is the common denominator of M, then we need to impose that
1179 * exists alpha : d M(x) = d alpha
1181 * This function is similar to add_strides in isl_morph.c
1183 static __isl_give isl_basic_map
*add_strides(__isl_take isl_basic_map
*bmap
,
1184 __isl_keep isl_mat
*M
, int n_known
)
1189 if (isl_int_is_one(M
->row
[0][0]))
1192 bmap
= isl_basic_map_extend_space(bmap
, isl_space_copy(bmap
->dim
),
1193 M
->n_row
- 1, M
->n_row
- 1, 0);
1196 for (i
= 1; i
< M
->n_row
; ++i
) {
1197 isl_seq_gcd(M
->row
[i
], M
->n_col
, &gcd
);
1198 if (isl_int_is_divisible_by(gcd
, M
->row
[0][0]))
1200 div
= isl_basic_map_alloc_div(bmap
);
1203 isl_int_set_si(bmap
->div
[div
][0], 0);
1204 k
= isl_basic_map_alloc_equality(bmap
);
1207 isl_seq_cpy(bmap
->eq
[k
], M
->row
[i
], M
->n_col
);
1208 isl_seq_clr(bmap
->eq
[k
] + M
->n_col
, bmap
->n_div
- n_known
);
1209 isl_int_set(bmap
->eq
[k
][M
->n_col
- n_known
+ div
],
1217 isl_basic_map_free(bmap
);
1221 /* If there are any equalities that involve (multiple) unknown divs,
1222 * then extract the stride information encoded by those equalities
1223 * and make it explicitly available in "bmap".
1225 * We first sort the divs so that the unknown divs appear last and
1226 * then we count how many equalities involve these divs.
1228 * Let these equalities be of the form
1232 * where y represents the unknown divs and x the remaining variables.
1233 * Let [H 0] be the Hermite Normal Form of B, i.e.,
1237 * Then x is a solution of the equalities iff
1239 * H^-1 A(x) (= - [I 0] Q y)
1241 * is an integer vector. Let d be the common denominator of H^-1.
1244 * d H^-1 A(x) = d alpha
1246 * in add_strides, with alpha fresh existentially quantified variables.
1248 static __isl_give isl_basic_map
*isl_basic_map_make_strides_explicit(
1249 __isl_take isl_basic_map
*bmap
)
1258 known
= isl_basic_map_divs_known(bmap
);
1260 return isl_basic_map_free(bmap
);
1263 bmap
= isl_basic_map_sort_divs(bmap
);
1264 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1268 for (n_known
= 0; n_known
< bmap
->n_div
; ++n_known
)
1269 if (isl_int_is_zero(bmap
->div
[n_known
][0]))
1271 ctx
= isl_basic_map_get_ctx(bmap
);
1272 total
= isl_space_dim(bmap
->dim
, isl_dim_all
);
1273 for (n
= 0; n
< bmap
->n_eq
; ++n
)
1274 if (isl_seq_first_non_zero(bmap
->eq
[n
] + 1 + total
+ n_known
,
1275 bmap
->n_div
- n_known
) == -1)
1279 B
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 0, 1 + total
+ n_known
);
1280 n_col
= bmap
->n_div
- n_known
;
1281 A
= isl_mat_sub_alloc6(ctx
, bmap
->eq
, 0, n
, 1 + total
+ n_known
, n_col
);
1282 A
= isl_mat_left_hermite(A
, 0, NULL
, NULL
);
1283 A
= isl_mat_drop_cols(A
, n
, n_col
- n
);
1284 A
= isl_mat_lin_to_aff(A
);
1285 A
= isl_mat_right_inverse(A
);
1286 B
= isl_mat_insert_zero_rows(B
, 0, 1);
1287 B
= isl_mat_set_element_si(B
, 0, 0, 1);
1288 M
= isl_mat_product(A
, B
);
1290 return isl_basic_map_free(bmap
);
1291 bmap
= add_strides(bmap
, M
, n_known
);
1292 bmap
= isl_basic_map_gauss(bmap
, NULL
);
1298 /* Compute the affine hull of each basic map in "map" separately
1299 * and make all stride information explicit so that we can remove
1300 * all unknown divs without losing this information.
1301 * The result is also guaranteed to be gaussed.
1303 * In simple cases where a div is determined by an equality,
1304 * calling isl_basic_map_gauss is enough to make the stride information
1305 * explicit, as it will derive an explicit representation for the div
1306 * from the equality. If, however, the stride information
1307 * is encoded through multiple unknown divs then we need to make
1308 * some extra effort in isl_basic_map_make_strides_explicit.
1310 static __isl_give isl_map
*isl_map_local_affine_hull(__isl_take isl_map
*map
)
1314 map
= isl_map_cow(map
);
1318 for (i
= 0; i
< map
->n
; ++i
) {
1319 map
->p
[i
] = isl_basic_map_affine_hull(map
->p
[i
]);
1320 map
->p
[i
] = isl_basic_map_gauss(map
->p
[i
], NULL
);
1321 map
->p
[i
] = isl_basic_map_make_strides_explicit(map
->p
[i
]);
1323 return isl_map_free(map
);
1329 static __isl_give isl_set
*isl_set_local_affine_hull(__isl_take isl_set
*set
)
1331 return isl_map_local_affine_hull(set
);
1334 /* Compute the affine hull of "map".
1336 * We first compute the affine hull of each basic map separately.
1337 * Then we align the divs and recompute the affine hulls of the basic
1338 * maps since some of them may now have extra divs.
1339 * In order to avoid performing parametric integer programming to
1340 * compute explicit expressions for the divs, possible leading to
1341 * an explosion in the number of basic maps, we first drop all unknown
1342 * divs before aligning the divs. Note that isl_map_local_affine_hull tries
1343 * to make sure that all stride information is explicitly available
1344 * in terms of known divs. This involves calling isl_basic_set_gauss,
1345 * which is also needed because affine_hull assumes its input has been gaussed,
1346 * while isl_map_affine_hull may be called on input that has not been gaussed,
1347 * in particular from initial_facet_constraint.
1348 * Similarly, align_divs may reorder some divs so that we need to
1349 * gauss the result again.
1350 * Finally, we combine the individual affine hulls into a single
1353 __isl_give isl_basic_map
*isl_map_affine_hull(__isl_take isl_map
*map
)
1355 struct isl_basic_map
*model
= NULL
;
1356 struct isl_basic_map
*hull
= NULL
;
1357 struct isl_set
*set
;
1358 isl_basic_set
*bset
;
1360 map
= isl_map_detect_equalities(map
);
1361 map
= isl_map_local_affine_hull(map
);
1362 map
= isl_map_remove_empty_parts(map
);
1363 map
= isl_map_remove_unknown_divs(map
);
1364 map
= isl_map_align_divs(map
);
1370 hull
= isl_basic_map_empty_like_map(map
);
1375 model
= isl_basic_map_copy(map
->p
[0]);
1376 set
= isl_map_underlying_set(map
);
1377 set
= isl_set_cow(set
);
1378 set
= isl_set_local_affine_hull(set
);
1383 set
->p
[0] = affine_hull(set
->p
[0], set
->p
[--set
->n
]);
1385 bset
= isl_basic_set_copy(set
->p
[0]);
1386 hull
= isl_basic_map_overlying_set(bset
, model
);
1388 hull
= isl_basic_map_simplify(hull
);
1389 return isl_basic_map_finalize(hull
);
1391 isl_basic_map_free(model
);
1396 struct isl_basic_set
*isl_set_affine_hull(struct isl_set
*set
)
1398 return (struct isl_basic_set
*)
1399 isl_map_affine_hull((struct isl_map
*)set
);