isl_convex_hull.c: wrap_facet: allow unbounded facets again.
[platform/upstream/isl.git] / isl_convex_hull.c
index 1980cca..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"
@@ -50,7 +59,8 @@ int isl_basic_map_constraint_is_redundant(struct isl_basic_map **bmap,
        if (i < total)
                return 0;
 
-       res = isl_solve_lp(*bmap, 0, c, (*bmap)->ctx->one, opt_n, opt_d);
+       res = isl_basic_map_solve_lp(*bmap, 0, c, (*bmap)->ctx->one,
+                                       opt_n, opt_d, NULL);
        if (res == isl_lp_unbounded)
                return 0;
        if (res == isl_lp_error)
@@ -93,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)
@@ -128,8 +143,8 @@ static int uset_is_bound(struct isl_set *set, isl_int *c, unsigned len)
                if (ISL_F_ISSET(set->p[j], ISL_BASIC_SET_EMPTY))
                        continue;
 
-               res = isl_solve_lp((struct isl_basic_map*)set->p[j],
-                               0, c, set->ctx->one, &opt, &opt_denom);
+               res = isl_basic_set_solve_lp(set->p[j],
+                               0, c, set->ctx->one, &opt, &opt_denom, NULL);
                if (res == isl_lp_unbounded)
                        break;
                if (res == isl_lp_error)
@@ -140,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);
@@ -191,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];
@@ -281,14 +298,13 @@ static struct isl_basic_set *isl_basic_set_add_equality(
        struct isl_basic_set *bset, isl_int *c)
 {
        int i;
-       unsigned total;
        unsigned dim;
 
        if (ISL_F_ISSET(bset, ISL_BASIC_SET_EMPTY))
                return bset;
 
-       isl_assert(ctx, isl_basic_set_n_param(bset) == 0, goto error);
-       isl_assert(ctx, bset->n_div == 0, goto error);
+       isl_assert(bset->ctx, isl_basic_set_n_param(bset) == 0, goto error);
+       isl_assert(bset->ctx, bset->n_div == 0, goto error);
        dim = isl_basic_set_n_dim(bset);
        bset = isl_basic_set_cow(bset);
        bset = isl_basic_set_extend(bset, 0, dim, 0, 1, 0);
@@ -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)
 {
@@ -481,8 +498,8 @@ static isl_int *wrap_facet(struct isl_set *set, isl_int *facet, isl_int *ridge)
        }
        isl_int_init(num);
        isl_int_init(den);
-       res = isl_solve_lp((struct isl_basic_map *)lp, 0,
-                                   obj->block.data, set->ctx->one, &num, &den);
+       res = isl_basic_set_solve_lp(lp, 0,
+                           obj->block.data, set->ctx->one, &num, &den, NULL);
        if (res == isl_lp_ok) {
                isl_int_neg(num, num);
                isl_seq_combine(facet, num, facet, den, ridge, dim);
@@ -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".
  *
@@ -522,12 +589,12 @@ static struct isl_mat *initial_facet_constraint(struct isl_set *set,
        int i;
        unsigned dim = isl_set_n_dim(set);
 
-       isl_assert(ctx, set->n > 0, goto error);
-       isl_assert(ctx, bounds->n_row == dim, goto error);
+       isl_assert(set->ctx, set->n > 0, goto error);
+       isl_assert(set->ctx, bounds->n_row == dim, goto error);
 
        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(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;
@@ -666,9 +730,11 @@ static struct isl_basic_set *extend(struct isl_basic_set *hull,
        int k;
        struct isl_basic_set *facet = NULL;
        struct isl_basic_set *hull_facet = NULL;
-       unsigned total;
        unsigned dim;
 
+       if (!hull)
+               return NULL;
+
        isl_assert(set->ctx, set->n > 0, goto error);
 
        dim = isl_set_n_dim(set);
@@ -938,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;
@@ -1206,9 +1272,9 @@ static struct isl_vec *valid_direction(
                                bset1->ctx->one, dir->block.data,
                                sample->block.data[n++], bset1->ineq[i], 1 + d);
        isl_vec_free(sample);
+       isl_seq_normalize(bset1->ctx, dir->block.data + 1, dir->size - 1);
        isl_basic_set_free(bset1);
        isl_basic_set_free(bset2);
-       isl_seq_normalize(dir->block.data + 1, dir->size - 1);
        return dir;
 error:
        isl_vec_free(sample);
@@ -1333,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);
         
@@ -1374,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);
@@ -1452,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);
@@ -1494,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);
@@ -1568,7 +1634,7 @@ static void update_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
        struct max_constraint *c;
        uint32_t c_hash;
 
-       c_hash = isl_seq_hash(con + 1, len, isl_hash_init());
+       c_hash = isl_seq_get_hash(con + 1, len);
        entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
                        con + 1, 0);
        if (!entry)
@@ -1601,7 +1667,7 @@ static int has_constraint(struct isl_ctx *ctx, struct isl_hash_table *table,
        struct max_constraint *c;
        uint32_t c_hash;
 
-       c_hash = isl_seq_hash(con + 1, len, isl_hash_init());
+       c_hash = isl_seq_get_hash(con + 1, len);
        entry = isl_hash_table_find(ctx, table, c_hash, max_constraint_equal,
                        con + 1, 0);
        if (!entry)
@@ -1668,8 +1734,7 @@ static struct isl_basic_set *common_constraints(struct isl_basic_set *hull,
        for (i = 0; i < min_constraints; ++i) {
                struct isl_hash_table_entry *entry;
                uint32_t c_hash;
-               c_hash = isl_seq_hash(constraints[i].c->row[0] + 1, total,
-                                       isl_hash_init());
+               c_hash = isl_seq_get_hash(constraints[i].c->row[0] + 1, total);
                entry = isl_hash_table_find(hull->ctx, table, c_hash,
                        max_constraint_equal, constraints[i].c->row[0] + 1, 1);
                if (!entry)
@@ -1786,7 +1851,6 @@ static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set)
  */
 static struct isl_basic_set *uset_convex_hull(struct isl_set *set)
 {
-       int i;
        struct isl_basic_set *convex_hull = NULL;
        struct isl_basic_set *lin;
 
@@ -1835,7 +1899,6 @@ error:
  */
 static struct isl_basic_set *uset_convex_hull_wrap_bounded(struct isl_set *set)
 {
-       int i;
        struct isl_basic_set *convex_hull = NULL;
 
        if (isl_set_n_dim(set) == 0) {
@@ -1972,7 +2035,7 @@ struct sh_data {
        struct isl_ctx          *ctx;
        unsigned                n;
        struct isl_hash_table   *hull_table;
-       struct sh_data_entry    p[0];
+       struct sh_data_entry    p[1];
 };
 
 static void sh_data_free(struct sh_data *data)
@@ -2012,7 +2075,7 @@ static int hash_ineq(struct isl_ctx *ctx, struct isl_hash_table *table,
 
        v.len = len;
        v.p = ineq;
-       c_hash = isl_seq_hash(ineq + 1, len, isl_hash_init());
+       c_hash = isl_seq_get_hash(ineq + 1, len);
        entry = isl_hash_table_find(ctx, table, c_hash, has_ineq, &v, 1);
        if (!entry)
                return - 1;
@@ -2050,7 +2113,8 @@ static struct sh_data *sh_data_alloc(struct isl_set *set, unsigned n_ineq)
        int i;
 
        data = isl_calloc(set->ctx, struct sh_data,
-               sizeof(struct sh_data) + set->n * sizeof(struct sh_data_entry));
+               sizeof(struct sh_data) +
+               (set->n - 1) * sizeof(struct sh_data_entry));
        if (!data)
                return NULL;
        data->ctx = set->ctx;
@@ -2135,7 +2199,7 @@ static struct isl_basic_set *add_bound(struct isl_basic_set *hull,
 
        v.len = isl_basic_set_total_dim(hull);
        v.p = ineq;
-       c_hash = isl_seq_hash(ineq + 1, v.len, isl_hash_init());
+       c_hash = isl_seq_get_hash(ineq + 1, v.len);
 
        entry = isl_hash_table_find(hull->ctx, data->hull_table, c_hash,
                                        has_ineq, &v, 0);
@@ -2238,7 +2302,7 @@ static struct isl_basic_set *uset_simple_hull(struct isl_set *set)
        struct sh_data *data = NULL;
        struct isl_basic_set *hull = NULL;
        unsigned n_ineq;
-       int i, j;
+       int i;
 
        if (!set)
                return NULL;