add isl_map_power and isl_map_transitive_closure
authorSven Verdoolaege <skimo@kotnet.org>
Wed, 3 Feb 2010 17:29:52 +0000 (18:29 +0100)
committerSven Verdoolaege <skimo@kotnet.org>
Wed, 3 Feb 2010 17:41:55 +0000 (18:41 +0100)
Makefile.am
doc/user.pod
include/isl_map.h
isl_test.c
isl_transitive_closure.c [new file with mode: 0644]

index 0f3ca74..85815bb 100644 (file)
@@ -72,6 +72,7 @@ libisl_la_SOURCES = \
        isl_tab.c \
        isl_tab.h \
        isl_tab_pip.c \
+       isl_transitive_closure.c \
        isl_vec.c
 EXTRA_libisl_la_SOURCES = \
        isl_lp_piplib.c \
index 1bfb878..e976d9d 100644 (file)
@@ -767,6 +767,29 @@ variables, then the result of these operations is currently undefined.
        __isl_give isl_basic_map *isl_map_affine_hull(
                __isl_take isl_map *map);
 
+=item * Power
+
+       __isl_give isl_map *isl_map_power(__isl_take isl_map *map,
+               unsigned param, int *exact);
+
+Compute a parametric representation for all positive powers I<k> of C<map>.
+The power I<k> is equated to the parameter at position C<param>.
+The result may be an overapproximation.  If the result is exact,
+then C<*exact> is set to C<1>.
+The current implementation only produces exact results for particular
+cases of piecewise translations (i.e., piecewise uniform dependences).
+
+=item * Transitive closure
+
+       __isl_give isl_map *isl_map_transitive_closure(
+               __isl_take isl_map *map, int *exact);
+
+Compute the transitive closure of C<map>.
+The result may be an overapproximation.  If the result is exact,
+then C<*exact> is set to C<1>.
+The current implementation only produces exact results for particular
+cases of piecewise translations (i.e., piecewise uniform dependences).
+
 =back
 
 =head2 Binary Operations
index 7a50710..c78cfc3 100644 (file)
@@ -357,6 +357,11 @@ int isl_map_foreach_basic_map(__isl_keep isl_map *map,
 
 __isl_give isl_map *isl_set_lifting(__isl_take isl_set *set);
 
+__isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
+       int *exact);
+__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
+       int *exact);
+
 #if defined(__cplusplus)
 }
 #endif
index 8119f42..04f199e 100644 (file)
@@ -601,6 +601,114 @@ void test_coalesce(struct isl_ctx *ctx)
        isl_set_free(set);
 }
 
+void test_closure(struct isl_ctx *ctx)
+{
+       isl_set *dom;
+       isl_map *up, *right;
+       isl_map *map, *map2;
+       int exact;
+
+       // COCOA example 1
+       map = isl_map_read_from_str(ctx,
+               "[n,k] -> { [i,j] -> [i2,j2] : i2 = i + 1 and j2 = j + 1 and "
+                       "1 <= i and i < n and 1 <= j and j < n or "
+                       "i2 = i + 1 and j2 = j - 1 and "
+                       "1 <= i and i < n and 2 <= j and j <= n }", -1);
+       map = isl_map_power(map, 1, &exact);
+       assert(exact);
+       isl_map_free(map);
+
+       // COCOA example 1
+       map = isl_map_read_from_str(ctx,
+               "[n] -> { [i,j] -> [i2,j2] : i2 = i + 1 and j2 = j + 1 and "
+                       "1 <= i and i < n and 1 <= j and j < n or "
+                       "i2 = i + 1 and j2 = j - 1 and "
+                       "1 <= i and i < n and 2 <= j and j <= n }", -1);
+       map = isl_map_transitive_closure(map, &exact);
+       assert(exact);
+       map2 = isl_map_read_from_str(ctx,
+               "[n] -> { [i,j] -> [i2,j2] : exists (k1,k2,k : "
+                       "1 <= i and i < n and 1 <= j and j <= n and "
+                       "2 <= i2 and i2 <= n and 1 <= j2 and j2 <= n and "
+                       "i2 = i + k1 + k2 and j2 = j + k1 - k2 and "
+                       "k1 >= 0 and k2 >= 0 and k1 + k2 = k and k >= 1 )}", -1);
+       assert(isl_map_is_equal(map, map2));
+       isl_map_free(map2);
+       isl_map_free(map);
+
+       map = isl_map_read_from_str(ctx,
+               "[n] -> { [x] -> [y] : y = x + 1 and 0 <= x and x <= n and "
+                                    " 0 <= y and y <= n }", -1);
+       map = isl_map_transitive_closure(map, &exact);
+       map2 = isl_map_read_from_str(ctx,
+               "[n] -> { [x] -> [y] : y > x and 0 <= x and x <= n and "
+                                    " 0 <= y and y <= n }", -1);
+       assert(isl_map_is_equal(map, map2));
+       isl_map_free(map2);
+       isl_map_free(map);
+
+       // COCOA example 2
+       map = isl_map_read_from_str(ctx,
+               "[n] -> { [i,j] -> [i2,j2] : i2 = i + 2 and j2 = j + 2 and "
+                       "1 <= i and i < n - 1 and 1 <= j and j < n - 1 or "
+                       "i2 = i + 2 and j2 = j - 2 and "
+                       "1 <= i and i < n - 1 and 3 <= j and j <= n }", -1);
+       map = isl_map_transitive_closure(map, &exact);
+       assert(exact);
+       map2 = isl_map_read_from_str(ctx,
+               "[n] -> { [i,j] -> [i2,j2] : exists (k1,k2,k : "
+                       "1 <= i and i < n - 1 and 1 <= j and j <= n and "
+                       "3 <= i2 and i2 <= n and 1 <= j2 and j2 <= n and "
+                       "i2 = i + 2 k1 + 2 k2 and j2 = j + 2 k1 - 2 k2 and "
+                       "k1 >= 0 and k2 >= 0 and k1 + k2 = k and k >= 1) }", -1);
+       assert(isl_map_is_equal(map, map2));
+       isl_map_free(map);
+       isl_map_free(map2);
+
+       // COCOA Fig.2 left
+       map = isl_map_read_from_str(ctx,
+               "[n] -> { [i,j] -> [i2,j2] : i2 = i + 2 and j2 = j and "
+                       "i <= 2 j - 3 and i <= n - 2 and j <= 2 i - 1 and "
+                       "j <= n or "
+                       "i2 = i and j2 = j + 2 and i <= 2 j - 1 and i <= n and "
+                       "j <= 2 i - 3 and j <= n - 2 or "
+                       "i2 = i + 1 and j2 = j + 1 and i <= 2 j - 1 and "
+                       "i <= n - 1 and j <= 2 i - 1 and j <= n - 1 }", -1);
+       map = isl_map_transitive_closure(map, &exact);
+       assert(exact);
+       isl_map_free(map);
+
+       // COCOA Fig.2 right
+       map = isl_map_read_from_str(ctx,
+               "[n] -> { [i,j] -> [i2,j2] : i2 = i + 3 and j2 = j and "
+                       "i <= 2 j - 4 and i <= n - 3 and j <= 2 i - 1 and "
+                       "j <= n or "
+                       "i2 = i and j2 = j + 3 and i <= 2 j - 1 and i <= n and "
+                       "j <= 2 i - 4 and j <= n - 3 or "
+                       "i2 = i + 1 and j2 = j + 1 and i <= 2 j - 1 and "
+                       "i <= n - 1 and j <= 2 i - 1 and j <= n - 1 }", -1);
+       map = isl_map_transitive_closure(map, &exact);
+       assert(!exact);
+       isl_map_free(map);
+
+       // COCOA Fig.1 right
+       dom = isl_set_read_from_str(ctx,
+               "{ [x,y] : x >= 0 and -2 x + 3 y >= 0 and x <= 3 and "
+                       "2 x - 3 y + 3 >= 0 }", -1);
+       right = isl_map_read_from_str(ctx,
+               "{ [x,y] -> [x2,y2] : x2 = x + 1 and y2 = y }", -1);
+       up = isl_map_read_from_str(ctx,
+               "{ [x,y] -> [x2,y2] : x2 = x and y2 = y + 1 }", -1);
+       right = isl_map_intersect_domain(right, isl_set_copy(dom));
+       right = isl_map_intersect_range(right, isl_set_copy(dom));
+       up = isl_map_intersect_domain(up, isl_set_copy(dom));
+       up = isl_map_intersect_range(up, dom);
+       map = isl_map_union(up, right);
+       map = isl_map_transitive_closure(map, &exact);
+       assert(!exact);
+       isl_map_free(map);
+}
+
 int main()
 {
        struct isl_ctx *ctx;
@@ -618,6 +726,7 @@ int main()
        test_convex_hull(ctx);
        test_gist(ctx);
        test_coalesce(ctx);
+       test_closure(ctx);
        isl_ctx_free(ctx);
        return 0;
 }
diff --git a/isl_transitive_closure.c b/isl_transitive_closure.c
new file mode 100644 (file)
index 0000000..6095d6a
--- /dev/null
@@ -0,0 +1,304 @@
+/*
+ * Copyright 2010      INRIA Saclay
+ *
+ * Use of this software is governed by the GNU LGPLv2.1 license
+ *
+ * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
+ * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
+ * 91893 Orsay, France 
+ */
+
+#include "isl_map.h"
+#include "isl_map_private.h"
+/*
+ * The transitive closure implementation is based on the paper
+ * "Computing the Transitive Closure of a Union of Affine Integer
+ * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
+ * Albert Cohen.
+ */
+
+/* Given a union of translations (uniform dependences), return a matrix
+ * with as many rows as there are disjuncts in the union and as many
+ * columns as there are input dimensions (which should be equal to
+ * the number of output dimensions).
+ * Each row contains the translation performed by the corresponding disjunct.
+ * If "map" turns out not to be a union of translations, then the contents
+ * of the returned matrix are undefined and *ok is set to 0.
+ */
+static __isl_give isl_mat *extract_steps(__isl_keep isl_map *map, int *ok)
+{
+       int i, j;
+       struct isl_mat *steps;
+       unsigned dim = isl_map_dim(map, isl_dim_in);
+
+       *ok = 1;
+
+       steps = isl_mat_alloc(map->ctx, map->n, dim);
+       if (!steps)
+               return NULL;
+
+       for (i = 0; i < map->n; ++i) {
+               struct isl_basic_set *delta;
+
+               delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
+
+               for (j = 0; j < dim; ++j) {
+                       int fixed;
+
+                       fixed = isl_basic_set_fast_dim_is_fixed(delta, j,
+                                                           &steps->row[i][j]);
+                       if (fixed < 0) {
+                               isl_basic_set_free(delta);
+                               goto error;
+                       }
+                       if (!fixed)
+                               break;
+               }
+
+               isl_basic_set_free(delta);
+
+               if (j < dim)
+                       break;
+       }
+
+       if (i < map->n)
+               *ok = 0;
+
+       return steps;
+error:
+       isl_mat_free(steps);
+       return NULL;
+}
+
+/* Given a set of n offsets v_i (the rows of "steps"), construct a relation
+ * of the given dimension specification that maps a element x to any
+ * element that can be reached by taking a positive number of steps
+ * along any of the offsets, where the number of steps k is equal to
+ * parameter "param".  That is, construct
+ *
+ * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v_i, k = \sum_i k_i > 0 }
+ *
+ * If strict is non-negative, then at least one step should be taken
+ * along the given offset v_strict, i.e., k_strict > 0.
+ */
+static __isl_give isl_map *path_along_steps(__isl_take isl_dim *dim,
+       __isl_keep isl_mat *steps, unsigned param, int strict)
+{
+       int i, j, k;
+       struct isl_basic_map *path = NULL;
+       unsigned d;
+       unsigned n;
+       unsigned nparam;
+
+       if (!dim || !steps)
+               goto error;
+
+       d = isl_dim_size(dim, isl_dim_in);
+       n = steps->n_row;
+       nparam = isl_dim_size(dim, isl_dim_param);
+
+       path = isl_basic_map_alloc_dim(isl_dim_copy(dim), n, d + 1, n + 1);
+
+       for (i = 0; i < n; ++i) {
+               k = isl_basic_map_alloc_div(path);
+               if (k < 0)
+                       goto error;
+               isl_assert(steps->ctx, i == k, goto error);
+               isl_int_set_si(path->div[k][0], 0);
+       }
+
+       for (i = 0; i < d; ++i) {
+               k = isl_basic_map_alloc_equality(path);
+               if (k < 0)
+                       goto error;
+               isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
+               isl_int_set_si(path->eq[k][1 + nparam + i], 1);
+               isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
+               for (j = 0; j < n; ++j)
+                       isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
+                                   steps->row[j][i]);
+       }
+
+       k = isl_basic_map_alloc_equality(path);
+       if (k < 0)
+               goto error;
+       isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
+       isl_int_set_si(path->eq[k][1 + param], -1);
+       for (j = 0; j < n; ++j)
+               isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
+
+       for (i = 0; i < n; ++i) {
+               k = isl_basic_map_alloc_inequality(path);
+               if (k < 0)
+                       goto error;
+               isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
+               isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
+               if (i == strict)
+                       isl_int_set_si(path->ineq[k][0], -1);
+       }
+
+       k = isl_basic_map_alloc_inequality(path);
+       if (k < 0)
+               goto error;
+       isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
+       isl_int_set_si(path->ineq[k][1 + param], 1);
+       isl_int_set_si(path->ineq[k][0], -1);
+
+       isl_dim_free(dim);
+
+       path = isl_basic_map_simplify(path);
+       path = isl_basic_map_finalize(path);
+       return isl_map_from_basic_map(path);
+error:
+       isl_dim_free(dim);
+       isl_basic_map_free(path);
+       return NULL;
+}
+
+/* Check whether the overapproximation of the power of "map" is exactly
+ * the power of "map".  In particular, for each path of a given length
+ * that starts of in domain or range and ends up in the range,
+ * check whether there is at least one path of the same length
+ * with a valid first segment, i.e., one in "map".
+ *
+ * "domain" and "range" are the domain and range of "map"
+ * "steps" represents the translations of "map"
+ * "path" is a path along "steps"
+ *
+ * "domain", "range", "steps" and "path" have been precomputed by the calling
+ * function.
+ */
+static int check_exactness(__isl_take isl_map *map, __isl_take isl_set *domain,
+       __isl_take isl_set *range, __isl_take isl_map *path,
+       __isl_keep isl_mat *steps, unsigned param)
+{
+       isl_map *test;
+       int ok;
+       int i;
+
+       if (!map)
+               goto error;
+
+       test = isl_map_empty(isl_map_get_dim(map));
+       for (i = 0; i < map->n; ++i) {
+               struct isl_map *path_i;
+               struct isl_set *dom_i;
+               path_i = path_along_steps(isl_map_get_dim(map), steps, param, i);
+               dom_i = isl_set_from_basic_set(
+                       isl_basic_map_domain(isl_basic_map_copy(map->p[i])));
+               path_i = isl_map_intersect_domain(path_i, dom_i);
+               test = isl_map_union(test, path_i);
+       }
+       isl_map_free(map);
+       test = isl_map_intersect_range(test, isl_set_copy(range));
+
+       domain = isl_set_union(domain, isl_set_copy(range));
+       path = isl_map_intersect_domain(path, domain);
+       path = isl_map_intersect_range(path, range);
+
+       ok = isl_map_is_subset(path, test);
+
+       isl_map_free(path);
+       isl_map_free(test);
+
+       return ok;
+error:
+       isl_map_free(map);
+       isl_set_free(domain);
+       isl_set_free(range);
+       isl_map_free(path);
+       return -1;
+}
+
+/* Compute the positive powers of "map", or an overapproximation.
+ * The power is given by parameter "param".  If the result is exact,
+ * then *exact is set to 1.
+ */
+__isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param,
+       int *exact)
+{
+       struct isl_mat *steps = NULL;
+       struct isl_set *domain = NULL;
+       struct isl_set *range = NULL;
+       struct isl_map *app = NULL;
+       struct isl_map *path = NULL;
+       int ok;
+
+       if (exact)
+               *exact = 1;
+
+       map = isl_map_remove_empty_parts(map);
+       if (!map)
+               return NULL;
+
+       if (isl_map_fast_is_empty(map))
+               return map;
+
+       isl_assert(map->ctx, param < isl_map_dim(map, isl_dim_param), goto error);
+       isl_assert(map->ctx,
+               isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
+               goto error);
+
+       domain = isl_map_domain(isl_map_copy(map));
+       range = isl_map_range(isl_map_copy(map));
+       app = isl_map_from_domain_and_range(isl_set_copy(domain),
+                                           isl_set_copy(range));
+
+       steps = extract_steps(map, &ok);
+       if (!ok)
+               goto not_exact;
+
+       path = path_along_steps(isl_map_get_dim(map), steps, param, -1);
+       app = isl_map_intersect(app, isl_map_copy(path));
+
+       if (exact &&
+           (*exact = check_exactness(isl_map_copy(map), isl_set_copy(domain),
+                                 isl_set_copy(range), isl_map_copy(path),
+                                 steps, param)) < 0)
+               goto error;
+
+       if (0) {
+not_exact:
+               if (exact)
+                       *exact = 0;
+       }
+       isl_set_free(domain);
+       isl_set_free(range);
+       isl_mat_free(steps);
+       isl_map_free(path);
+       isl_map_free(map);
+       return app;
+error:
+       isl_set_free(domain);
+       isl_set_free(range);
+       isl_mat_free(steps);
+       isl_map_free(path);
+       isl_map_free(map);
+       isl_map_free(app);
+       return NULL;
+}
+
+/* Compute the transitive closure  of "map", or an overapproximation.
+ * If the result is exact, then *exact is set to 1.
+ * Simply compute the powers of map and then project out the parameter
+ * describing the power.
+ */
+__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
+       int *exact)
+{
+       unsigned param;
+
+       if (!map)
+               goto error;
+
+       param = isl_map_dim(map, isl_dim_param);
+       map = isl_map_add(map, isl_dim_param, 1);
+       map = isl_map_power(map, param, exact);
+       map = isl_map_project_out(map, isl_dim_param, param, 1);
+
+       return map;
+error:
+       isl_map_free(map);
+       return NULL;
+}