+/*
+ * Copyright 2008-2009 Katholieke Universiteit Leuven
+ *
+ * Use of this software is governed by the GNU LGPLv2.1 license
+ *
+ * Written by Sven Verdoolaege, K.U.Leuven, Departement
+ * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
+ */
+
#include "isl_lp.h"
#include "isl_map.h"
#include "isl_map_private.h"
return bmap;
tab = isl_tab_from_basic_map(bmap);
- tab = isl_tab_detect_equalities(tab);
- tab = isl_tab_detect_redundant(tab);
+ tab = isl_tab_detect_implicit_equalities(tab);
+ if (isl_tab_detect_redundant(tab) < 0)
+ goto error;
bmap = isl_basic_map_update_from_tab(bmap, tab);
isl_tab_free(tab);
ISL_F_SET(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
ISL_F_SET(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
return bmap;
+error:
+ isl_tab_free(tab);
+ isl_basic_map_free(bmap);
+ return NULL;
}
struct isl_basic_set *isl_basic_set_convex_hull(struct isl_basic_set *bset)
goto error;
continue;
}
- if (!isl_int_is_one(opt_denom))
- isl_seq_scale(c, c, opt_denom, len);
- if (first || isl_int_is_neg(opt))
+ if (first || isl_int_is_neg(opt)) {
+ if (!isl_int_is_one(opt_denom))
+ isl_seq_scale(c, c, opt_denom, len);
isl_int_sub(c[0], c[0], opt);
+ }
first = 0;
}
isl_int_clear(opt);
is_bound = uset_is_bound(set, dirs->row[n], dirs->n_col);
if (is_bound != 1)
return is_bound;
+ isl_seq_normalize(set->ctx, dirs->row[n], dirs->n_col);
if (i < n) {
int k;
isl_int *t = dirs->row[n];
return NULL;
}
-static struct isl_set *isl_set_add_equality(struct isl_set *set, isl_int *c)
+static struct isl_set *isl_set_add_basic_set_equality(struct isl_set *set, isl_int *c)
{
int i;
* In the original space, we need to take the same combination of the
* corresponding constraints "facet" and "ridge".
*
- * Note that a is always finite, since we only apply the wrapping
- * technique to a union of polytopes.
+ * If a = -infty = "-1/0", then we just return the original facet constraint.
+ * This means that the facet is unbounded, but has a bounded intersection
+ * with the union of sets.
*/
static isl_int *wrap_facet(struct isl_set *set, isl_int *facet, isl_int *ridge)
{
isl_vec_free(obj);
isl_basic_set_free(lp);
isl_set_free(set);
- isl_assert(set->ctx, res == isl_lp_ok, return NULL);
+ isl_assert(set->ctx, res == isl_lp_ok || res == isl_lp_unbounded,
+ return NULL);
return facet;
error:
isl_basic_set_free(lp);
return NULL;
}
+/* Drop rows in "rows" that are redundant with respect to earlier rows,
+ * assuming that "rows" is of full column rank.
+ *
+ * We compute the column echelon form. The non-redundant rows are
+ * those that are the first to contain a non-zero entry in a column.
+ * All the other rows can be removed.
+ */
+static __isl_give isl_mat *drop_redundant_rows(__isl_take isl_mat *rows)
+{
+ struct isl_mat *H = NULL;
+ int col;
+ int row;
+ int last_row;
+
+ if (!rows)
+ return NULL;
+
+ isl_assert(rows->ctx, rows->n_row >= rows->n_col, goto error);
+
+ if (rows->n_row == rows->n_col)
+ return rows;
+
+ H = isl_mat_left_hermite(isl_mat_copy(rows), 0, NULL, NULL);
+ if (!H)
+ goto error;
+
+ last_row = rows->n_row;
+ for (col = rows->n_col - 1; col >= 0; --col) {
+ for (row = col; row < last_row; ++row)
+ if (!isl_int_is_zero(H->row[row][col]))
+ break;
+ isl_assert(rows->ctx, row < last_row, goto error);
+ if (row + 1 < last_row) {
+ rows = isl_mat_drop_rows(rows, row + 1, last_row - (row + 1));
+ if (rows->n_row == rows->n_col)
+ break;
+ }
+ last_row = row;
+ }
+
+ isl_mat_free(H);
+
+ return rows;
+error:
+ isl_mat_free(H);
+ isl_mat_free(rows);
+ return NULL;
+}
+
/* Given a set of d linearly independent bounding constraints of the
* convex hull of "set", compute the constraint of a facet of "set".
*
while (bounds->n_row > 1) {
slice = isl_set_copy(set);
- slice = isl_set_add_equality(slice, bounds->row[0]);
+ slice = isl_set_add_basic_set_equality(slice, bounds->row[0]);
face = isl_set_affine_hull(slice);
if (!face)
goto error;
U = isl_mat_drop_cols(U, 0, 1);
Q = isl_mat_drop_rows(Q, 0, 1);
bounds = isl_mat_product(bounds, U);
+ bounds = drop_redundant_rows(bounds);
bounds = isl_mat_product(bounds, Q);
- while (isl_seq_first_non_zero(bounds->row[bounds->n_row-1],
- bounds->n_col) == -1) {
- bounds->n_row--;
- isl_assert(set->ctx, bounds->n_row > 1, goto error);
- }
+ isl_assert(set->ctx, bounds->n_row > 1, goto error);
if (!wrap_facet(set, bounds->row[0],
bounds->row[bounds->n_row-1]))
goto error;
struct isl_basic_set *hull_facet = NULL;
unsigned dim;
+ if (!hull)
+ return NULL;
+
isl_assert(set->ctx, set->n > 0, goto error);
dim = isl_set_n_dim(set);
bset1 = homogeneous_map(bset1, isl_mat_copy(T2));
bset2 = homogeneous_map(bset2, T2);
set = isl_set_alloc_dim(isl_basic_set_get_dim(bset1), 2, 0);
- set = isl_set_add(set, bset1);
- set = isl_set_add(set, bset2);
+ set = isl_set_add_basic_set(set, bset1);
+ set = isl_set_add_basic_set(set, bset2);
hull = uset_convex_hull(set);
hull = isl_basic_set_preimage(hull, T);
if (lin->n_eq < isl_basic_set_total_dim(lin)) {
struct isl_set *set;
set = isl_set_alloc_dim(isl_basic_set_get_dim(bset1), 2, 0);
- set = isl_set_add(set, bset1);
- set = isl_set_add(set, bset2);
+ set = isl_set_add_basic_set(set, bset1);
+ set = isl_set_add_basic_set(set, bset2);
return modulo_lineality(set, lin);
}
isl_basic_set_free(lin);
lin = isl_set_alloc_dim(isl_set_get_dim(set), set->n, 0);
for (i = 0; i < set->n; ++i)
- lin = isl_set_add(lin,
+ lin = isl_set_add_basic_set(lin,
isl_basic_set_lineality_space(isl_basic_set_copy(set->p[i])));
isl_set_free(set);
return isl_set_affine_hull(lin);
break;
}
if (t->n_eq < isl_basic_set_total_dim(t)) {
- set = isl_set_add(set, convex_hull);
+ set = isl_set_add_basic_set(set, convex_hull);
return modulo_lineality(set, t);
}
isl_basic_set_free(t);