From ad73f0c3a7ba2dc71049d5601bd65423742c25fe Mon Sep 17 00:00:00 2001
From: Sven Verdoolaege <skimo@kotnet.org>
Date: Sun, 17 Aug 2008 09:45:12 +0200
Subject: [PATCH] isl_map_affine_hull: perform affine hull on underlying set

---
 isl_affine_hull.c | 153 ++++++++++++++++++++++++++++++------------------------
 isl_map.c         | 151 +++++++++++++++++++++++++++++++++++++++++------------
 isl_map_private.h |   5 ++
 3 files changed, 209 insertions(+), 100 deletions(-)

diff --git a/isl_affine_hull.c b/isl_affine_hull.c
index ad9cbaf..e8df67c 100644
--- a/isl_affine_hull.c
+++ b/isl_affine_hull.c
@@ -75,42 +75,42 @@ struct isl_basic_set *isl_basic_set_affine_hull(struct isl_ctx *ctx,
  * after column col are zero.
  */
 static void set_common_multiple(
-	struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
+	struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 	unsigned row, unsigned col)
 {
 	isl_int m, c;
 
-	if (isl_int_eq(bmap1->eq[row][col], bmap2->eq[row][col]))
+	if (isl_int_eq(bset1->eq[row][col], bset2->eq[row][col]))
 		return;
 
 	isl_int_init(c);
 	isl_int_init(m);
-	isl_int_lcm(m, bmap1->eq[row][col], bmap2->eq[row][col]);
-	isl_int_divexact(c, m, bmap1->eq[row][col]);
-	isl_seq_scale(bmap1->eq[row], bmap1->eq[row], c, col+1);
-	isl_int_divexact(c, m, bmap2->eq[row][col]);
-	isl_seq_scale(bmap2->eq[row], bmap2->eq[row], c, col+1);
+	isl_int_lcm(m, bset1->eq[row][col], bset2->eq[row][col]);
+	isl_int_divexact(c, m, bset1->eq[row][col]);
+	isl_seq_scale(bset1->eq[row], bset1->eq[row], c, col+1);
+	isl_int_divexact(c, m, bset2->eq[row][col]);
+	isl_seq_scale(bset2->eq[row], bset2->eq[row], c, col+1);
 	isl_int_clear(c);
 	isl_int_clear(m);
 }
 
 /* Delete a given equality, moving all the following equalities one up.
  */
-static void delete_row(struct isl_basic_map *bmap, unsigned row)
+static void delete_row(struct isl_basic_set *bset, unsigned row)
 {
 	isl_int *t;
 	int r;
 
-	t = bmap->eq[row];
-	bmap->n_eq--;
-	for (r = row; r < bmap->n_eq; ++r)
-		bmap->eq[r] = bmap->eq[r+1];
-	bmap->eq[bmap->n_eq] = t;
+	t = bset->eq[row];
+	bset->n_eq--;
+	for (r = row; r < bset->n_eq; ++r)
+		bset->eq[r] = bset->eq[r+1];
+	bset->eq[bset->n_eq] = t;
 }
 
-/* Make first row entries in column col of bmap1 identical to
- * those of bmap2, using the fact that entry bmap1->eq[row][col]=a
- * is non-zero.  Initially, these elements of bmap1 are all zero.
+/* Make first row entries in column col of bset1 identical to
+ * those of bset2, using the fact that entry bset1->eq[row][col]=a
+ * is non-zero.  Initially, these elements of bset1 are all zero.
  * For each row i < row, we set
  *		A[i] = a * A[i] + B[i][col] * A[row]
  *		B[i] = a * B[i]
@@ -118,7 +118,7 @@ static void delete_row(struct isl_basic_map *bmap, unsigned row)
  *		A[i][col] = B[i][col] = a * old(B[i][col])
  */
 static void construct_column(
-	struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
+	struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 	unsigned row, unsigned col)
 {
 	int r;
@@ -128,24 +128,24 @@ static void construct_column(
 
 	isl_int_init(a);
 	isl_int_init(b);
-	total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
+	total = 1 + bset1->dim;
 	for (r = 0; r < row; ++r) {
-		if (isl_int_is_zero(bmap2->eq[r][col]))
+		if (isl_int_is_zero(bset2->eq[r][col]))
 			continue;
-		isl_int_gcd(b, bmap2->eq[r][col], bmap1->eq[row][col]);
-		isl_int_divexact(a, bmap1->eq[row][col], b);
-		isl_int_divexact(b, bmap2->eq[r][col], b);
-		isl_seq_combine(bmap1->eq[r], a, bmap1->eq[r],
-					      b, bmap1->eq[row], total);
-		isl_seq_scale(bmap2->eq[r], bmap2->eq[r], a, total);
+		isl_int_gcd(b, bset2->eq[r][col], bset1->eq[row][col]);
+		isl_int_divexact(a, bset1->eq[row][col], b);
+		isl_int_divexact(b, bset2->eq[r][col], b);
+		isl_seq_combine(bset1->eq[r], a, bset1->eq[r],
+					      b, bset1->eq[row], total);
+		isl_seq_scale(bset2->eq[r], bset2->eq[r], a, total);
 	}
 	isl_int_clear(a);
 	isl_int_clear(b);
-	delete_row(bmap1, row);
+	delete_row(bset1, row);
 }
 
-/* Make first row entries in column col of bmap1 identical to
- * those of bmap2, using only these entries of the two matrices.
+/* Make first row entries in column col of bset1 identical to
+ * those of bset2, using only these entries of the two matrices.
  * Let t be the last row with different entries.
  * For each row i < t, we set
  *	A[i] = (A[t][col]-B[t][col]) * A[i] + (B[i][col]-A[i][col) * A[t]
@@ -154,7 +154,7 @@ static void construct_column(
  *	A[i][col] = B[i][col] = old(A[t][col]*B[i][col]-A[i][col]*B[t][col])
  */
 static int transform_column(
-	struct isl_basic_map *bmap1, struct isl_basic_map *bmap2,
+	struct isl_basic_set *bset1, struct isl_basic_set *bset2,
 	unsigned row, unsigned col)
 {
 	int i, t;
@@ -162,31 +162,31 @@ static int transform_column(
 	unsigned total;
 
 	for (t = row-1; t >= 0; --t)
-		if (isl_int_ne(bmap1->eq[t][col], bmap2->eq[t][col]))
+		if (isl_int_ne(bset1->eq[t][col], bset2->eq[t][col]))
 			break;
 	if (t < 0)
 		return 0;
 
-	total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
+	total = 1 + bset1->dim;
 	isl_int_init(a);
 	isl_int_init(b);
 	isl_int_init(g);
-	isl_int_sub(b, bmap1->eq[t][col], bmap2->eq[t][col]);
+	isl_int_sub(b, bset1->eq[t][col], bset2->eq[t][col]);
 	for (i = 0; i < t; ++i) {
-		isl_int_sub(a, bmap2->eq[i][col], bmap1->eq[i][col]);
+		isl_int_sub(a, bset2->eq[i][col], bset1->eq[i][col]);
 		isl_int_gcd(g, a, b);
 		isl_int_divexact(a, a, g);
 		isl_int_divexact(g, b, g);
-		isl_seq_combine(bmap1->eq[i], g, bmap1->eq[i], a, bmap1->eq[t],
+		isl_seq_combine(bset1->eq[i], g, bset1->eq[i], a, bset1->eq[t],
 				total);
-		isl_seq_combine(bmap2->eq[i], g, bmap2->eq[i], a, bmap2->eq[t],
+		isl_seq_combine(bset2->eq[i], g, bset2->eq[i], a, bset2->eq[t],
 				total);
 	}
 	isl_int_clear(a);
 	isl_int_clear(b);
 	isl_int_clear(g);
-	delete_row(bmap1, t);
-	delete_row(bmap2, t);
+	delete_row(bset1, t);
+	delete_row(bset2, t);
 	return 1;
 }
 
@@ -195,74 +195,91 @@ static int transform_column(
  * except that the echelon form we use starts from the last column
  * and that we are dealing with integer coefficients.
  */
-static struct isl_basic_map *affine_hull(struct isl_ctx *ctx,
-	struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
+static struct isl_basic_set *affine_hull(struct isl_ctx *ctx,
+	struct isl_basic_set *bset1, struct isl_basic_set *bset2)
 {
 	unsigned total;
 	int col;
 	int row;
 
-	total = 1 + bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
+	total = 1 + bset1->dim;
 
 	row = 0;
 	for (col = total-1; col >= 0; --col) {
-		int is_zero1 = row >= bmap1->n_eq ||
-			isl_int_is_zero(bmap1->eq[row][col]);
-		int is_zero2 = row >= bmap2->n_eq ||
-			isl_int_is_zero(bmap2->eq[row][col]);
+		int is_zero1 = row >= bset1->n_eq ||
+			isl_int_is_zero(bset1->eq[row][col]);
+		int is_zero2 = row >= bset2->n_eq ||
+			isl_int_is_zero(bset2->eq[row][col]);
 		if (!is_zero1 && !is_zero2) {
-			set_common_multiple(bmap1, bmap2, row, col);
+			set_common_multiple(bset1, bset2, row, col);
 			++row;
 		} else if (!is_zero1 && is_zero2) {
-			construct_column(bmap1, bmap2, row, col);
+			construct_column(bset1, bset2, row, col);
 		} else if (is_zero1 && !is_zero2) {
-			construct_column(bmap2, bmap1, row, col);
+			construct_column(bset2, bset1, row, col);
 		} else {
-			if (transform_column(bmap1, bmap2, row, col))
+			if (transform_column(bset1, bset2, row, col))
 				--row;
 		}
 	}
-	isl_basic_map_free(ctx, bmap2);
-	return bmap1;
+	isl_basic_set_free(ctx, bset2);
+	return bset1;
 }
 
 struct isl_basic_map *isl_map_affine_hull(struct isl_ctx *ctx,
 						struct isl_map *map)
 {
 	int i;
-	struct isl_basic_map *bmap;
+	struct isl_basic_map *model = NULL;
+	struct isl_basic_map *hull = NULL;
+	struct isl_set *set;
 
 	if (!map)
 		return NULL;
 
 	if (map->n == 0) {
-		bmap = isl_basic_map_empty(ctx,
+		hull = isl_basic_map_empty(ctx,
 					    map->nparam, map->n_in, map->n_out);
 		isl_map_free(ctx, map);
-		return bmap;
+		return hull;
 	}
 
 	map = isl_map_align_divs(ctx, map);
-	map = isl_map_cow(ctx, map);
+	model = isl_basic_map_copy(ctx, map->p[0]);
+	set = isl_map_underlying_set(ctx, map);
+	set = isl_set_cow(ctx, set);
+	if (!set)
+		goto error;
 
-	for (i = 0; i < map->n; ++i) {
-		map->p[i] = isl_basic_map_cow(ctx, map->p[i]);
-		map->p[i] = isl_basic_map_affine_hull(ctx, map->p[i]);
-		map->p[i] = isl_basic_map_gauss(ctx, map->p[i], NULL);
-		if (!map->p[i])
+	for (i = 0; i < set->n; ++i) {
+		set->p[i] = isl_basic_set_cow(ctx, set->p[i]);
+		set->p[i] = isl_basic_set_affine_hull(ctx, set->p[i]);
+		set->p[i] = isl_basic_set_gauss(ctx, set->p[i], NULL);
+		if (!set->p[i])
 			goto error;
 	}
-	while (map->n > 1) {
-		map->p[0] = affine_hull(ctx, map->p[0], map->p[--map->n]);
-		if (!map->p[0])
-			goto error;
+	set = isl_set_remove_empty_parts(ctx, set);
+	if (set->n == 0) {
+		hull = isl_basic_map_empty(ctx,
+				    model->nparam, model->n_in, model->n_out);
+		isl_basic_map_free(ctx, model);
+	} else {
+		struct isl_basic_set *bset;
+		while (set->n > 1) {
+			set->p[0] = affine_hull(ctx, set->p[0],
+							set->p[--set->n]);
+			if (!set->p[0])
+				goto error;
+		}
+		bset = isl_basic_set_copy(ctx, set->p[0]);
+		hull = isl_basic_map_overlying_set(ctx, bset, model);
 	}
-	bmap = isl_basic_map_copy(ctx, map->p[0]);
-	isl_map_free(ctx, map);
-	bmap = isl_basic_map_simplify(ctx, bmap);
-	return isl_basic_map_finalize(ctx, bmap);
+	isl_set_free(ctx, set);
+	hull = isl_basic_map_simplify(ctx, hull);
+	return isl_basic_map_finalize(ctx, hull);
 error:
-	isl_map_free(ctx, map);
+	isl_basic_map_free(ctx, model);
+	isl_set_free(ctx, set);
 	return NULL;
 }
 
diff --git a/isl_map.c b/isl_map.c
index 033a1d3..fc001e9 100644
--- a/isl_map.c
+++ b/isl_map.c
@@ -489,7 +489,7 @@ struct isl_basic_map *isl_basic_map_cow(struct isl_ctx *ctx,
 	return bmap;
 }
 
-static struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set)
+struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set)
 {
 	if (!set)
 		return NULL;
@@ -1776,6 +1776,37 @@ error:
 	return NULL;
 }
 
+/* For a div d = floor(f/m), add the constraints
+ *
+ *		f - m d >= 0
+ *		-(f-(n-1)) + m d >= 0
+ *
+ * Note that the second constraint is the negation of
+ *
+ *		f - m d >= n
+ */
+static int add_div_constraints(struct isl_ctx *ctx,
+		struct isl_basic_map *bmap, unsigned div)
+{
+	int i, j;
+	unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
+	unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
+
+	i = isl_basic_map_alloc_inequality(ctx, bmap);
+	if (i < 0)
+		return -1;
+	isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
+	isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
+
+	j = isl_basic_map_alloc_inequality(ctx, bmap);
+	if (j < 0)
+		return -1;
+	isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
+	isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
+	isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
+	return j;
+}
+
 static struct isl_basic_set *isl_basic_map_underlying_set(
 		struct isl_ctx *ctx, struct isl_basic_map *bmap)
 {
@@ -1796,6 +1827,93 @@ error:
 	return NULL;
 }
 
+struct isl_basic_map *isl_basic_map_overlying_set(
+	struct isl_ctx *ctx, struct isl_basic_set *bset,
+	struct isl_basic_map *like)
+{
+	struct isl_basic_map *bmap;
+	unsigned total;
+	int i, k;
+
+	if (!bset || !like)
+		goto error;
+	isl_assert(ctx, bset->dim ==
+			like->nparam + like->n_in + like->n_out + like->n_div,
+			goto error);
+	if (like->nparam == 0 && like->n_in == 0 && like->n_div == 0) {
+		isl_basic_map_free(ctx, like);
+		return (struct isl_basic_map *)bset;
+	}
+	bset = isl_basic_set_cow(ctx, bset);
+	total = bset->dim + bset->extra;
+	if (!bset)
+		goto error;
+	bmap = (struct isl_basic_map *)bset;
+	bmap->nparam = like->nparam;
+	bmap->n_in = like->n_in;
+	bmap->n_out = like->n_out;
+	bmap->extra += like->n_div;
+	if (bmap->extra) {
+		unsigned ltotal;
+		ltotal = total - bmap->extra + like->extra;
+		if (ltotal > total)
+			ltotal = total;
+		bmap->block2 = isl_blk_extend(ctx, bmap->block2,
+					bmap->extra * (1 + 1 + total));
+		if (isl_blk_is_error(bmap->block2))
+			goto error;
+		bmap->div = isl_realloc_array(ctx, bmap->div, isl_int *,
+						bmap->extra);
+		if (!bmap->div)
+			goto error;
+		bmap = isl_basic_map_extend(ctx, bmap, bmap->nparam,
+			bmap->n_in, bmap->n_out, 0, 0, 2 * like->n_div);
+		for (i = 0; i < like->n_div; ++i) {
+			k = isl_basic_map_alloc_div(ctx, bmap);
+			if (k < 0)
+				goto error;
+			isl_seq_cpy(bmap->div[k], like->div[i], 1 + 1 + ltotal);
+			isl_seq_clr(bmap->div[k]+1+1+ltotal, total - ltotal);
+			if (add_div_constraints(ctx, bmap, k) < 0)
+				goto error;
+		}
+	}
+	isl_basic_map_free(ctx, like);
+	bmap = isl_basic_map_finalize(ctx, bmap);
+	return bmap;
+error:
+	isl_basic_map_free(ctx, like);
+	isl_basic_set_free(ctx, bset);
+	return NULL;
+}
+
+struct isl_set *isl_map_underlying_set(struct isl_ctx *ctx, struct isl_map *map)
+{
+	int i;
+
+	map = isl_map_align_divs(ctx, map);
+	map = isl_map_cow(ctx, map);
+	if (!map)
+		return NULL;
+
+	for (i = 0; i < map->n; ++i) {
+		map->p[i] = (struct isl_basic_map *)
+				isl_basic_map_underlying_set(ctx, map->p[i]);
+		if (!map->p[i])
+			goto error;
+	}
+	if (map->n == 0)
+		map->n_out += map->nparam + map->n_in;
+	else
+		map->n_out = map->p[0]->n_out;
+	map->nparam = 0;
+	map->n_in = 0;
+	return (struct isl_set *)map;
+error:
+	isl_map_free(ctx, map);
+	return NULL;
+}
+
 struct isl_basic_set *isl_basic_map_domain(struct isl_ctx *ctx,
 		struct isl_basic_map *bmap)
 {
@@ -2793,37 +2911,6 @@ static int find_div(struct isl_basic_map *dst,
 	return -1;
 }
 
-/* For a div d = floor(f/m), add the constraints
- *
- *		f - m d >= 0
- *		-(f-(n-1)) + m d >= 0
- *
- * Note that the second constraint is the negation of
- *
- *		f - m d >= n
- */
-static int add_div_constraints(struct isl_ctx *ctx,
-		struct isl_basic_map *bmap, unsigned div)
-{
-	int i, j;
-	unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
-	unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
-
-	i = isl_basic_map_alloc_inequality(ctx, bmap);
-	if (i < 0)
-		return -1;
-	isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
-	isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
-
-	j = isl_basic_map_alloc_inequality(ctx, bmap);
-	if (j < 0)
-		return -1;
-	isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
-	isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
-	isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
-	return j;
-}
-
 struct isl_basic_map *isl_basic_map_align_divs(struct isl_ctx *ctx,
 		struct isl_basic_map *dst, struct isl_basic_map *src)
 {
diff --git a/isl_map_private.h b/isl_map_private.h
index fd713b9..7f980e0 100644
--- a/isl_map_private.h
+++ b/isl_map_private.h
@@ -27,6 +27,7 @@ struct isl_basic_set *isl_basic_set_cow(struct isl_ctx *ctx,
 		struct isl_basic_set *bset);
 struct isl_basic_map *isl_basic_map_cow(struct isl_ctx *ctx,
 		struct isl_basic_map *bmap);
+struct isl_set *isl_set_cow(struct isl_ctx *ctx, struct isl_set *set);
 struct isl_map *isl_map_cow(struct isl_ctx *ctx, struct isl_map *map);
 
 struct isl_basic_map *isl_basic_map_set_to_empty(
@@ -43,6 +44,10 @@ struct isl_basic_map *isl_basic_map_gauss(struct isl_ctx *ctx,
 	struct isl_basic_map *bmap, int *progress);
 struct isl_basic_set *isl_basic_set_gauss(struct isl_ctx *ctx,
 	struct isl_basic_set *bset, int *progress);
+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,
+	struct isl_basic_map *like);
 
 struct isl_map *isl_map_remove_empty_parts(struct isl_ctx *ctx,
 	struct isl_map *map);
-- 
2.7.4