isl_convex_hull.c: wrap_facet: allow unbounded facets again.
[platform/upstream/isl.git] / isl_convex_hull.c
index bbe06bd..e91bf12 100644 (file)
@@ -1,3 +1,12 @@
+/*
+ * 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"
@@ -94,13 +103,18 @@ struct isl_basic_map *isl_basic_map_convex_hull(struct isl_basic_map *bmap)
                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)
@@ -141,10 +155,11 @@ static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len)
                                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);
@@ -192,6 +207,7 @@ static int is_independent_bound(struct isl_set *set, isl_int *c,
        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];
@@ -302,7 +318,7 @@ error:
        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;
 
@@ -441,8 +457,9 @@ static struct isl_basic_set *wrap_constraints(struct isl_set *set)
  * 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)
 {
@@ -492,7 +509,8 @@ 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);
@@ -501,6 +519,55 @@ error:
        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".
  *
@@ -527,7 +594,7 @@ static struct isl_mat *initial_facet_constraint(struct isl_set *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;
@@ -549,12 +616,9 @@ static struct isl_mat *initial_facet_constraint(struct isl_set *set,
                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;
@@ -668,6 +732,9 @@ static struct isl_basic_set *extend(struct isl_basic_set *hull,
        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);
@@ -937,7 +1004,7 @@ static int isl_basic_set_is_bounded(struct isl_basic_set *bset)
        struct isl_tab *tab;
        int bounded;
 
-       tab = isl_tab_from_recession_cone((struct isl_basic_map *)bset);
+       tab = isl_tab_from_recession_cone(bset);
        bounded = isl_tab_cone_is_bounded(tab);
        isl_tab_free(tab);
        return bounded;
@@ -1332,8 +1399,8 @@ static struct isl_basic_set *convex_hull_pair_pointed(
        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);
         
@@ -1373,8 +1440,8 @@ static struct isl_basic_set *convex_hull_pair(struct isl_basic_set *bset1,
        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);
@@ -1451,7 +1518,7 @@ static struct isl_basic_set *uset_combined_lineality_space(struct isl_set *set)
 
        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);
@@ -1493,7 +1560,7 @@ static struct isl_basic_set *uset_convex_hull_unbounded(struct isl_set *set)
                        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);