From b266652bfdd1b5e9075ed5febb2036f1a5c62d84 Mon Sep 17 00:00:00 2001
From: Sven Verdoolaege <skimo@kotnet.org>
Date: Fri, 22 Aug 2008 15:29:16 +0200
Subject: [PATCH] isl_basic_map_affine_hull: compute integer affine hull

---
 include/isl_set.h           |   2 +
 isl_affine_hull.c           | 174 ++++++++++++++++++++++++++++++++++++++------
 isl_map.c                   |  14 +++-
 isl_map_private.h           |   2 +
 isl_test.c                  |  10 ++-
 test_inputs/affine2.polylib |   9 +++
 6 files changed, 186 insertions(+), 25 deletions(-)
 create mode 100644 test_inputs/affine2.polylib

diff --git a/include/isl_set.h b/include/isl_set.h
index 48d474c..0539b1c 100644
--- a/include/isl_set.h
+++ b/include/isl_set.h
@@ -69,6 +69,8 @@ struct isl_basic_set *isl_basic_set_finalize(struct isl_ctx *ctx,
 void isl_basic_set_free(struct isl_ctx *ctx, struct isl_basic_set *bset);
 struct isl_basic_set *isl_basic_set_copy(struct isl_ctx *ctx,
 					struct isl_basic_set *bset);
+struct isl_basic_set *isl_basic_set_empty(struct isl_ctx *ctx,
+		unsigned nparam, unsigned dim);
 struct isl_basic_set *isl_basic_set_universe(struct isl_ctx *ctx,
 		unsigned nparam, unsigned dim);
 void isl_basic_set_dump(struct isl_ctx *ctx, struct isl_basic_set *bset,
diff --git a/isl_affine_hull.c b/isl_affine_hull.c
index 4dae661..43e71c4 100644
--- a/isl_affine_hull.c
+++ b/isl_affine_hull.c
@@ -4,6 +4,8 @@
 #include "isl_lp.h"
 #include "isl_map.h"
 #include "isl_map_private.h"
+#include "isl_equalities.h"
+#include "isl_sample.h"
 
 struct isl_basic_map *isl_basic_map_implicit_equalities(struct isl_ctx *ctx,
 						struct isl_basic_map *bmap)
@@ -47,28 +49,6 @@ error:
 	return NULL;
 }
 
-struct isl_basic_map *isl_basic_map_affine_hull(struct isl_ctx *ctx,
-						struct isl_basic_map *bmap)
-{
-	bmap = isl_basic_map_implicit_equalities(ctx, bmap);
-	if (!bmap)
-		return NULL;
-	if (bmap->n_ineq == 0)
-		return bmap;
-	bmap = isl_basic_map_cow(ctx, bmap);
-	if (!bmap)
-		return NULL;
-	isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
-	return isl_basic_map_finalize(ctx, bmap);
-}
-
-struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
-						struct isl_basic_set *bset)
-{
-	return (struct isl_basic_set *)
-		isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
-}
-
 /* Make eq[row][col] of both bmaps equal so we can add the row
  * add the column to the common matrix.
  * Note that because of the echelon form, the columns of row row
@@ -223,7 +203,157 @@ static struct isl_basic_set *affine_hull(struct isl_ctx *ctx,
 		}
 	}
 	isl_basic_set_free(ctx, bset2);
+	isl_assert(ctx, row == bset1->n_eq, goto error);
 	return bset1;
+error:
+	isl_basic_set_free(ctx, bset1);
+	return NULL;
+}
+
+static struct isl_basic_set *isl_basic_set_from_vec(struct isl_ctx *ctx,
+	struct isl_vec *vec)
+{
+	int i;
+	int k;
+	struct isl_basic_set *bset = NULL;
+
+	if (!vec)
+		return NULL;
+	isl_assert(ctx, vec->size != 0, goto error);
+
+	bset = isl_basic_set_alloc(ctx, 0, vec->size - 1, 0, vec->size - 1, 0);
+	for (i = bset->dim - 1; i >= 0; --i) {
+		k = isl_basic_set_alloc_equality(ctx, bset);
+		if (k < 0)
+			goto error;
+		isl_seq_clr(bset->eq[k], 1 + bset->dim);
+		isl_int_neg(bset->eq[k][0], vec->block.data[1 + i]);
+		isl_int_set(bset->eq[k][1 + i], vec->block.data[0]);
+	}
+	isl_vec_free(ctx, vec);
+
+	return bset;
+error:
+	isl_basic_set_free(ctx, bset);
+	isl_vec_free(ctx, vec);
+	return NULL;
+}
+
+static struct isl_basic_set *outside_point(struct isl_ctx *ctx,
+	struct isl_basic_set *bset, isl_int *eq, int up)
+{
+	struct isl_basic_set *slice = NULL;
+	struct isl_vec *sample;
+	struct isl_basic_set *point;
+	int k;
+
+	slice = isl_basic_set_copy(ctx, bset);
+	if (!slice)
+		goto error;
+	slice = isl_basic_set_extend(ctx, slice, 0, slice->dim, 0, 1, 0);
+	k = isl_basic_set_alloc_equality(ctx, slice);
+	if (k < 0)
+		goto error;
+	isl_seq_cpy(slice->eq[k], eq, 1 + slice->dim);
+	if (up)
+		isl_int_add_ui(slice->eq[k][0], slice->eq[k][0], 1);
+	else
+		isl_int_sub_ui(slice->eq[k][0], slice->eq[k][0], 1);
+
+	sample = isl_basic_set_sample(ctx, slice);
+	if (!sample)
+		goto error;
+	if (sample->size == 0) {
+		isl_vec_free(ctx, sample);
+		point = isl_basic_set_empty(ctx, 0, bset->dim);
+	} else
+		point = isl_basic_set_from_vec(ctx, sample);
+
+	return point;
+error:
+	isl_basic_set_free(ctx, slice);
+	return NULL;
+}
+
+/* After computing the rational affine hull (by detecting the implicit
+ * equalities), we remove all equalities found so far, compute
+ * the integer affine hull of what is left, and then add the original
+ * equalities back in.
+ *
+ * The integer affine hull is constructed by successively looking
+ * a point that is affinely independent of the points found so far.
+ * In particular, for each equality satisfied by the points so far,
+ * we check if there is any point on the corresponding hyperplane
+ * shifted by one (in either direction).
+ */
+struct isl_basic_map *isl_basic_map_affine_hull(struct isl_ctx *ctx,
+						struct isl_basic_map *bmap)
+{
+	int i, j;
+	struct isl_mat *T2 = NULL;
+	struct isl_basic_set *bset = NULL;
+	struct isl_basic_set *hull = NULL;
+	struct isl_vec *sample;
+
+	bmap = isl_basic_map_implicit_equalities(ctx, bmap);
+	if (!bmap)
+		return NULL;
+	if (bmap->n_ineq == 0)
+		return bmap;
+
+	bset = isl_basic_map_underlying_set(ctx, isl_basic_map_copy(ctx, bmap));
+	bset = isl_basic_set_remove_equalities(ctx, bset, NULL, &T2);
+
+	sample = isl_basic_set_sample(ctx, isl_basic_set_copy(ctx, bset));
+	hull = isl_basic_set_from_vec(ctx, sample);
+
+	for (i = 0; i < bset->dim; ++i) {
+		struct isl_basic_set *point;
+		for (j = 0; j < hull->n_eq; ++j) {
+			point = outside_point(ctx, bset, hull->eq[j], 1);
+			if (!point)
+				goto error;
+			if (!F_ISSET(point, ISL_BASIC_SET_EMPTY))
+				break;
+			isl_basic_set_free(ctx, point);
+			point = outside_point(ctx, bset, hull->eq[j], 0);
+			if (!point)
+				goto error;
+			if (!F_ISSET(point, ISL_BASIC_SET_EMPTY))
+				break;
+			isl_basic_set_free(ctx, point);
+		}
+		if (j == hull->n_eq)
+			break;
+		hull = affine_hull(ctx, hull, point);
+	}
+
+	isl_basic_set_free(ctx, bset);
+	bset = NULL;
+	bmap = isl_basic_map_cow(ctx, bmap);
+	if (!bmap)
+		goto error;
+	isl_basic_map_free_inequality(ctx, bmap, bmap->n_ineq);
+	if (T2)
+		hull = isl_basic_set_preimage(ctx, hull, T2);
+	bmap = isl_basic_map_intersect(ctx, bmap,
+					isl_basic_map_overlying_set(ctx, hull,
+						isl_basic_map_copy(ctx, bmap)));
+
+	return isl_basic_map_finalize(ctx, bmap);
+error:
+	isl_mat_free(ctx, T2);
+	isl_basic_set_free(ctx, bset);
+	isl_basic_set_free(ctx, hull);
+	isl_basic_map_free(ctx, bmap);
+	return NULL;
+}
+
+struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
+						struct isl_basic_set *bset)
+{
+	return (struct isl_basic_set *)
+		isl_basic_map_affine_hull(ctx, (struct isl_basic_map *)bset);
 }
 
 struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
diff --git a/isl_map.c b/isl_map.c
index 67d1320..3c3b2e2 100644
--- a/isl_map.c
+++ b/isl_map.c
@@ -1996,7 +1996,7 @@ static int add_div_constraints(struct isl_ctx *ctx,
 	return j;
 }
 
-static struct isl_basic_set *isl_basic_map_underlying_set(
+struct isl_basic_set *isl_basic_map_underlying_set(
 		struct isl_ctx *ctx, struct isl_basic_map *bmap)
 {
 	if (!bmap)
@@ -2207,6 +2207,15 @@ struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
 	return bmap;
 }
 
+struct isl_basic_set *isl_basic_set_empty(struct isl_ctx *ctx,
+		unsigned nparam, unsigned dim)
+{
+	struct isl_basic_set *bset;
+	bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 1, 0);
+	bset = isl_basic_set_set_to_empty(ctx, bset);
+	return bset;
+}
+
 struct isl_basic_map *isl_basic_map_universe(struct isl_ctx *ctx,
 		unsigned nparam, unsigned in, unsigned out)
 {
@@ -2760,6 +2769,9 @@ static struct isl_basic_map *remove_redundant_divs(struct isl_ctx *ctx,
 {
 	int i;
 
+	if (!bmap)
+		return NULL;
+
 	for (i = bmap->n_div-1; i >= 0; --i) {
 		if (!div_is_redundant(ctx, bmap, i))
 			continue;
diff --git a/isl_map_private.h b/isl_map_private.h
index f3cef4e..ba73e14 100644
--- a/isl_map_private.h
+++ b/isl_map_private.h
@@ -52,6 +52,8 @@ struct isl_basic_set *isl_basic_set_gauss(struct isl_ctx *ctx,
 	struct isl_basic_set *bset, int *progress);
 struct isl_basic_map *isl_basic_map_implicit_equalities(struct isl_ctx *ctx,
 						struct isl_basic_map *bmap);
+struct isl_basic_set *isl_basic_map_underlying_set(
+		struct isl_ctx *ctx, struct isl_basic_map *bmap);
 struct isl_set *isl_map_underlying_set(struct isl_ctx *ctx, struct isl_map *map);
 struct isl_basic_map *isl_basic_map_overlying_set(
 	struct isl_ctx *ctx, struct isl_basic_set *bset,
diff --git a/isl_test.c b/isl_test.c
index 5616892..f54ce93 100644
--- a/isl_test.c
+++ b/isl_test.c
@@ -6,7 +6,7 @@
 
 static char *srcdir;
 
-void test_affine_hull(struct isl_ctx *ctx)
+void test_affine_hull_case(struct isl_ctx *ctx, const char *name)
 {
 	char filename[PATH_MAX];
 	FILE *input;
@@ -14,7 +14,7 @@ void test_affine_hull(struct isl_ctx *ctx)
 	struct isl_basic_set *bset1, *bset2;
 
 	n = snprintf(filename, sizeof(filename),
-			"%s/test_inputs/affine.polylib", srcdir);
+			"%s/test_inputs/%s.polylib", srcdir, name);
 	assert(n < sizeof(filename));
 	input = fopen(filename, "r");
 	assert(input);
@@ -32,6 +32,12 @@ void test_affine_hull(struct isl_ctx *ctx)
 	fclose(input);
 }
 
+void test_affine_hull(struct isl_ctx *ctx)
+{
+	test_affine_hull_case(ctx, "affine2");
+	test_affine_hull_case(ctx, "affine");
+}
+
 void test_convex_hull_case(struct isl_ctx *ctx, const char *name)
 {
 	char filename[PATH_MAX];
diff --git a/test_inputs/affine2.polylib b/test_inputs/affine2.polylib
new file mode 100644
index 0000000..c67db77
--- /dev/null
+++ b/test_inputs/affine2.polylib
@@ -0,0 +1,9 @@
+5 5
+1 -2 0 1 0
+1 2 0 -1 1
+1 0 -2 1 0
+1 0 2 -1 1
+1 0 0 1 -1
+
+1 5
+0 1 -1 0 0
-- 
2.7.4