privately export isl_set_is_bounded
[platform/upstream/isl.git] / isl_convex_hull.c
index b634353..17212c3 100644 (file)
@@ -155,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);
@@ -170,93 +171,6 @@ error:
        return -1;
 }
 
-/* Check if "c" is a direction that is independent of the previously found "n"
- * bounds in "dirs".
- * If so, add it to the list, with the negative of the lower bound
- * in the constant position, i.e., such that c corresponds to a bounding
- * hyperplane (but not necessarily a facet).
- * Assumes set "set" is bounded.
- */
-static int is_independent_bound(struct isl_set *set, isl_int *c,
-       struct isl_mat *dirs, int n)
-{
-       int is_bound;
-       int i = 0;
-
-       isl_seq_cpy(dirs->row[n]+1, c+1, dirs->n_col-1);
-       if (n != 0) {
-               int pos = isl_seq_first_non_zero(dirs->row[n]+1, dirs->n_col-1);
-               if (pos < 0)
-                       return 0;
-               for (i = 0; i < n; ++i) {
-                       int pos_i;
-                       pos_i = isl_seq_first_non_zero(dirs->row[i]+1, dirs->n_col-1);
-                       if (pos_i < pos)
-                               continue;
-                       if (pos_i > pos)
-                               break;
-                       isl_seq_elim(dirs->row[n]+1, dirs->row[i]+1, pos,
-                                       dirs->n_col-1, NULL);
-                       pos = isl_seq_first_non_zero(dirs->row[n]+1, dirs->n_col-1);
-                       if (pos < 0)
-                               return 0;
-               }
-       }
-
-       is_bound = uset_is_bound(set, dirs->row[n], dirs->n_col);
-       if (is_bound != 1)
-               return is_bound;
-       if (i < n) {
-               int k;
-               isl_int *t = dirs->row[n];
-               for (k = n; k > i; --k)
-                       dirs->row[k] = dirs->row[k-1];
-               dirs->row[i] = t;
-       }
-       return 1;
-}
-
-/* Compute and return a maximal set of linearly independent bounds
- * on the set "set", based on the constraints of the basic sets
- * in "set".
- */
-static struct isl_mat *independent_bounds(struct isl_set *set)
-{
-       int i, j, n;
-       struct isl_mat *dirs = NULL;
-       unsigned dim = isl_set_n_dim(set);
-
-       dirs = isl_mat_alloc(set->ctx, dim, 1+dim);
-       if (!dirs)
-               goto error;
-
-       n = 0;
-       for (i = 0; n < dim && i < set->n; ++i) {
-               int f;
-               struct isl_basic_set *bset = set->p[i];
-
-               for (j = 0; n < dim && j < bset->n_eq; ++j) {
-                       f = is_independent_bound(set, bset->eq[j], dirs, n);
-                       if (f < 0)
-                               goto error;
-                       if (f)
-                               ++n;
-               }
-               for (j = 0; n < dim && j < bset->n_ineq; ++j) {
-                       f = is_independent_bound(set, bset->ineq[j], dirs, n);
-                       if (f < 0)
-                               goto error;
-                       if (f)
-                               ++n;
-               }
-       }
-       dirs->n_row = n;
-       return dirs;
-error:
-       isl_mat_free(dirs);
-       return NULL;
-}
-
 struct isl_basic_set *isl_basic_set_set_rational(struct isl_basic_set *bset)
 {
        if (!bset)
@@ -316,7 +230,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;
 
@@ -455,10 +369,12 @@ 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)
+isl_int *isl_set_wrap_facet(__isl_keep isl_set *set,
+       isl_int *facet, isl_int *ridge)
 {
        int i;
        struct isl_mat *T = NULL;
@@ -469,6 +385,7 @@ static isl_int *wrap_facet(struct isl_set *set, isl_int *facet, isl_int *ridge)
        unsigned dim;
 
        set = isl_set_copy(set);
+       set = isl_set_set_rational(set);
 
        dim = 1 + isl_set_n_dim(set);
        T = isl_mat_alloc(set->ctx, 3, dim);
@@ -506,7 +423,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);
@@ -515,33 +433,44 @@ error:
        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".
- *
- * We first compute the intersection with the first bounding hyperplane
- * and remove the component corresponding to this hyperplane from
- * other bounds (in homogeneous space).
- * We then wrap around one of the remaining bounding constraints
- * and continue the process until all bounding constraints have been
- * taken into account.
- * The resulting linear combination of the bounding constraints will
- * correspond to a facet of the convex hull.
+/* Compute the constraint of a facet of "set".
+ *
+ * We first compute the intersection with a bounding constraint
+ * that is orthogonal to one of the coordinate axes.
+ * If the affine hull of this intersection has only one equality,
+ * we have found a facet.
+ * Otherwise, we wrap the current bounding constraint around
+ * one of the equalities of the face (one that is not equal to
+ * the current bounding constraint).
+ * This process continues until we have found a facet.
+ * The dimension of the intersection increases by at least
+ * one on each iteration, so termination is guaranteed.
  */
-static struct isl_mat *initial_facet_constraint(struct isl_set *set,
-       struct isl_mat *bounds)
+static __isl_give isl_mat *initial_facet_constraint(__isl_keep isl_set *set)
 {
        struct isl_set *slice = NULL;
        struct isl_basic_set *face = NULL;
        struct isl_mat *m, *U, *Q;
        int i;
        unsigned dim = isl_set_n_dim(set);
+       int is_bound;
+       isl_mat *bounds;
 
        isl_assert(set->ctx, set->n > 0, goto error);
-       isl_assert(set->ctx, bounds->n_row == dim, goto error);
+       bounds = isl_mat_alloc(set->ctx, 1, 1 + dim);
+       if (!bounds)
+               return NULL;
+
+       isl_seq_clr(bounds->row[0], dim);
+       isl_int_set_si(bounds->row[0][1 + dim - 1], 1);
+       is_bound = uset_is_bound(set, bounds->row[0], 1 + dim);
+       isl_assert(set->ctx, is_bound == 1, goto error);
+       isl_seq_normalize(set->ctx, bounds->row[0], 1 + dim);
+       bounds->n_row = 1;
 
-       while (bounds->n_row > 1) {
+       for (;;) {
                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,32 +478,18 @@ static struct isl_mat *initial_facet_constraint(struct isl_set *set,
                        isl_basic_set_free(face);
                        break;
                }
-               m = isl_mat_alloc(set->ctx, 1 + face->n_eq, 1 + dim);
-               if (!m)
-                       goto error;
-               isl_int_set_si(m->row[0][0], 1);
-               isl_seq_clr(m->row[0]+1, dim);
                for (i = 0; i < face->n_eq; ++i)
-                       isl_seq_cpy(m->row[1 + i], face->eq[i], 1 + dim);
-               U = isl_mat_right_inverse(m);
-               Q = isl_mat_right_inverse(isl_mat_copy(U));
-               U = isl_mat_drop_cols(U, 1 + face->n_eq, dim - face->n_eq);
-               Q = isl_mat_drop_rows(Q, 1 + face->n_eq, dim - face->n_eq);
-               U = isl_mat_drop_cols(U, 0, 1);
-               Q = isl_mat_drop_rows(Q, 0, 1);
-               bounds = isl_mat_product(bounds, U);
-               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);
-               }
-               if (!wrap_facet(set, bounds->row[0],
-                                         bounds->row[bounds->n_row-1]))
+                       if (!isl_seq_eq(bounds->row[0], face->eq[i], 1 + dim) &&
+                           !isl_seq_is_neg(bounds->row[0],
+                                               face->eq[i], 1 + dim))
+                               break;
+               isl_assert(set->ctx, i < face->n_eq, goto error);
+               if (!isl_set_wrap_facet(set, bounds->row[0], face->eq[i]))
                        goto error;
+               isl_seq_normalize(set->ctx, bounds->row[0], bounds->n_col);
                isl_basic_set_free(face);
-               bounds->n_row--;
        }
+
        return bounds;
 error:
        isl_basic_set_free(face);
@@ -682,6 +597,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);
@@ -711,7 +629,7 @@ static struct isl_basic_set *extend(struct isl_basic_set *hull,
                        if (k < 0)
                                goto error;
                        isl_seq_cpy(hull->ineq[k], hull->ineq[i], 1+dim);
-                       if (!wrap_facet(set, hull->ineq[k], facet->ineq[j]))
+                       if (!isl_set_wrap_facet(set, hull->ineq[k], facet->ineq[j]))
                                goto error;
                }
                isl_basic_set_free(hull_facet);
@@ -951,13 +869,18 @@ static int isl_basic_set_is_bounded(struct isl_basic_set *bset)
        struct isl_tab *tab;
        int bounded;
 
+       if (!bset)
+               return -1;
+       if (isl_basic_set_fast_is_empty(bset))
+               return 1;
+
        tab = isl_tab_from_recession_cone(bset);
        bounded = isl_tab_cone_is_bounded(tab);
        isl_tab_free(tab);
        return bounded;
 }
 
-static int isl_set_is_bounded(struct isl_set *set)
+int isl_set_is_bounded(__isl_keep isl_set *set)
 {
        int i;
 
@@ -1168,7 +1091,7 @@ error:
  * (including the "positivity constraint" 1 >= 0) and \alpha_{ij}
  * strictly positive numbers.  For simplicity we impose \alpha_{ij} >= 1.
  * We first set up an LP with as variables the \alpha{ij}.
- * In this formulateion, for each polyhedron i,
+ * In this formulation, for each polyhedron i,
  * the first constraint is the positivity constraint, followed by pairs
  * of variables for the equalities, followed by variables for the inequalities.
  * We then simply pick a feasible solution and compute s using (*).
@@ -1219,7 +1142,7 @@ 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_seq_normalize(bset1->ctx, dir->el, dir->size);
        isl_basic_set_free(bset1);
        isl_basic_set_free(bset2);
        return dir;
@@ -1346,8 +1269,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);
         
@@ -1361,18 +1284,47 @@ error:
        return NULL;
 }
 
+static struct isl_basic_set *uset_convex_hull_wrap(struct isl_set *set);
+static struct isl_basic_set *modulo_affine_hull(
+       struct isl_set *set, struct isl_basic_set *affine_hull);
+
 /* Compute the convex hull of a pair of basic sets without any parameters or
  * integer divisions.
  *
+ * This function is called from uset_convex_hull_unbounded, which
+ * means that the complete convex hull is unbounded.  Some pairs
+ * of basic sets may still be bounded, though.
+ * They may even lie inside a lower dimensional space, in which
+ * case they need to be handled inside their affine hull since
+ * the main algorithm assumes that the result is full-dimensional.
+ *
  * If the convex hull of the two basic sets would have a non-trivial
  * lineality space, we first project out this lineality space.
  */
 static struct isl_basic_set *convex_hull_pair(struct isl_basic_set *bset1,
        struct isl_basic_set *bset2)
 {
-       struct isl_basic_set *lin;
+       isl_basic_set *lin, *aff;
+       int bounded1, bounded2;
+
+       aff = isl_set_affine_hull(isl_basic_set_union(isl_basic_set_copy(bset1),
+                                                   isl_basic_set_copy(bset2)));
+       if (!aff)
+               goto error;
+       if (aff->n_eq != 0) 
+               return modulo_affine_hull(isl_basic_set_union(bset1, bset2), aff);
+       isl_basic_set_free(aff);
+
+       bounded1 = isl_basic_set_is_bounded(bset1);
+       bounded2 = isl_basic_set_is_bounded(bset2);
+
+       if (bounded1 < 0 || bounded2 < 0)
+               goto error;
 
-       if (isl_basic_set_is_bounded(bset1) || isl_basic_set_is_bounded(bset2))
+       if (bounded1 && bounded2)
+               uset_convex_hull_wrap(isl_basic_set_union(bset1, bset2));
+
+       if (bounded1 || bounded2)
                return convex_hull_pair_pointed(bset1, bset2);
 
        lin = induced_lineality_space(isl_basic_set_copy(bset1),
@@ -1387,8 +1339,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);
@@ -1465,7 +1417,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);
@@ -1507,7 +1459,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);
@@ -1521,13 +1473,8 @@ error:
 }
 
 /* Compute an initial hull for wrapping containing a single initial
- * facet by first computing bounds on the set and then using these
- * bounds to construct an initial facet.
- * This function is a remnant of an older implementation where the
- * bounds were also used to check whether the set was bounded.
- * Since this function will now only be called when we know the
- * set to be bounded, the initial facet should probably be constructed 
- * by simply using the coordinate directions instead.
+ * facet.
+ * This function assumes that the given set is bounded.
  */
 static struct isl_basic_set *initial_hull(struct isl_basic_set *hull,
        struct isl_set *set)
@@ -1538,11 +1485,7 @@ static struct isl_basic_set *initial_hull(struct isl_basic_set *hull,
 
        if (!hull)
                goto error;
-       bounds = independent_bounds(set);
-       if (!bounds)
-               goto error;
-       isl_assert(set->ctx, bounds->n_row == isl_set_n_dim(set), goto error);
-       bounds = initial_facet_constraint(set, bounds);
+       bounds = initial_facet_constraint(set);
        if (!bounds)
                goto error;
        k = isl_basic_set_alloc_inequality(hull);
@@ -1881,7 +1824,7 @@ error:
  * convex hull of the transformed set and then add the equalities back
  * (after performing the inverse transformation.
  */
-static struct isl_basic_set *modulo_affine_hull(struct isl_ctx *ctx,
+static struct isl_basic_set *modulo_affine_hull(
        struct isl_set *set, struct isl_basic_set *affine_hull)
 {
        struct isl_mat *T;
@@ -1940,7 +1883,7 @@ struct isl_basic_map *isl_map_convex_hull(struct isl_map *map)
        if (!affine_hull)
                goto error;
        if (affine_hull->n_eq != 0)
-               bset = modulo_affine_hull(ctx, set, affine_hull);
+               bset = modulo_affine_hull(set, affine_hull);
        else {
                isl_basic_set_free(affine_hull);
                bset = uset_convex_hull(set);
@@ -2112,7 +2055,7 @@ static int is_bound(struct sh_data *data, struct isl_set *set, int j,
 
        isl_int_clear(opt);
 
-       return res == isl_lp_ok ? 1 :
+       return (res == isl_lp_ok || res == isl_lp_empty) ? 1 :
               res == isl_lp_unbounded ? 0 : -1;
 }
 
@@ -2233,11 +2176,11 @@ static struct isl_basic_set *add_bounds(struct isl_basic_set *bset,
        for (j = 0; j < set->p[i]->n_eq; ++j) {
                for (k = 0; k < 2; ++k) {
                        isl_seq_neg(set->p[i]->eq[j], set->p[i]->eq[j], 1+dim);
-                       add_bound(bset, data, set, i, set->p[i]->eq[j]);
+                       bset = add_bound(bset, data, set, i, set->p[i]->eq[j]);
                }
        }
        for (j = 0; j < set->p[i]->n_ineq; ++j)
-               add_bound(bset, data, set, i, set->p[i]->ineq[j]);
+               bset = add_bound(bset, data, set, i, set->p[i]->ineq[j]);
        return bset;
 }