From: Sven Verdoolaege Date: Thu, 8 Apr 2010 08:04:46 +0000 (+0200) Subject: isl_map_closure: optionally use Omega-like algorithm X-Git-Tag: isl-0.03~255 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=94a70c9c2ad985d9e55499c8759fdce7ca668670;p=platform%2Fupstream%2Fisl.git isl_map_closure: optionally use Omega-like algorithm --- diff --git a/doc/implementation.tex b/doc/implementation.tex index 5f37778..55a158e 100644 --- a/doc/implementation.tex +++ b/doc/implementation.tex @@ -181,6 +181,7 @@ $$ \end{example} \subsection{Computing an Approximation of $R^k$} +\label{s:power} There are some special cases where the computation of $R^k$ is very easy. One such case is that where $R$ does not compose with itself, @@ -723,3 +724,112 @@ As explained by \shortciteN{Beletska2009}, applying the same formula to $R$ directly, without a decomposition, would result in an overapproximation of the power. \end{example} + +\subsection{An {\tt Omega}-like implementation} + +While the main algorithm of \shortciteN{Kelly1996closure} is +designed to compute and underapproximation of the transitive closure, +the authors mention that they could also compute overapproximations. +In this section, we describe our implementation of an algorithm +that is based on their ideas. + +The main tool is Equation~(2) of \shortciteN{Kelly1996closure}. +The input relation $R$ is first overapproximated by a ``d-form'' relation +$$ +\{\, \vec i \to \vec j \mid \exists \vec \alpha : +\vec L \le \vec j - \vec i \le \vec U +\wedge +(\forall p : j_p - i_p = M_p \alpha_p) +\,\} +, +$$ +where $p$ ranges over the dimensions and $\vec L$, $\vec U$ and +$\vec M$ are constant integer vectors. The elements of $\vec U$ +may be $\infty$, meaning that there is no upper bound corresponding +to that element, and similarly for $\vec L$. +Such an overapproximation can be obtained by computing strides, +lower and upper bounds on the difference set $\Delta \, R$. +The transitive closure of such a ``d-form'' relation is +\begin{equation} +\label{eq:omega} +\{\, \vec i \to \vec j \mid \exists \vec \alpha, k : +k \ge 1 \wedge +k \, \vec L \le \vec j - \vec i \le k \, \vec U +\wedge +(\forall p : j_p - i_p = M_p \alpha_p) +\,\} +. +\end{equation} +The domain and range of this transitive closure are then +intersected with those of the input relation. +This is a special case of the algorithm in \autoref{s:power}. + +In their algorithm for computing lower bounds, the authors +use the above algorithm as a substep on the disjuncts in the relation. +At the end, they say +\begin{quote} +If an upper bound is required, it can be calculated in a manner +similar to that of a single conjunct [sic] relation. +\end{quote} +Presumably, the authors mean that a ``d-form'' approximation +of the whole input relation should be used. +However, the accuracy can be improved by also using the following +idea from the same paper. If $R$ is a union of $m$ basic maps, +$$ +R = \bigcup_i R_i +, +$$ +and if we can find an $R_i$ such that for all other $R_j$ we have +that +$$ +R_i^* \circ R_j \circ R_i^* +$$ +can be represented as a single basic map, i.e., without a union, +then we can compute $R^+$ as +$$ +R^+ = R_i^+ \cup +\left( +\bigcup_{j \ne i} +R_i^* \circ R_j \circ R_i^* +\right)^+ +, +$$ +reducing the number of disjuncts in the argument of the transitive +closure by one. +An overapproximation of $R_i^*$ can be obtained by +allowing the value zero for $k$ in \eqref{eq:omega}, +i.e., by computing +$$ +\{\, \vec i \to \vec j \mid \exists \vec \alpha, k : +k \ge 0 \wedge +k \, \vec L \le \vec j - \vec i \le k \, \vec U +\wedge +(\forall p : j_p - i_p = M_p \alpha_p) +\,\} +. +$$ +However, when we intersect domain and range of this relation +with those of the input relation, then the result only contains +the idenity mapping on the intersection of domain and range. +\shortciteN{Kelly1996closure} propose to intersect domain +and range with then {\em union} of domain and range of the input +relation instead and call the result $R_i^?$. +Now, this union of domain and range of $R_i$ may not contain +the domains and ranges of the whole of $R$. +We can therefore not always replace +$R_i^* \circ R_j \circ R_i^*$ by +$R_i^? \circ R_j \circ R_i^?$. +\shortciteN{Kelly1996closure} propose to check the following +conditions to decide whether this replacement is justified: +$R_i^? - R_i^+$ is not a union and for each $j \ne i$ +the condition +$$ +\left(R_i^? - R_i^+\right) +\circ +R_j +\circ +\left(R_i^? - R_i^+\right) += +R_j +$$ +holds. diff --git a/doc/manual.tex b/doc/manual.tex index 5e2f6d5..e9c594f 100644 --- a/doc/manual.tex +++ b/doc/manual.tex @@ -25,6 +25,7 @@ \numberwithin{example}{section} \newcommand{\exampleautorefname}{Example} +\renewcommand{\subsectionautorefname}{Section} \def\Z{\mathbb{Z}} \def\Q{\mathbb{Q}} diff --git a/include/isl_options.h b/include/isl_options.h index 455b53f..b918d03 100644 --- a/include/isl_options.h +++ b/include/isl_options.h @@ -38,6 +38,10 @@ struct isl_options { #define ISL_GBR_ALWAYS 2 unsigned gbr; unsigned gbr_only_first; + + #define ISL_CLOSURE_ISL 0 + #define ISL_CLOSURE_OMEGA 1 + unsigned closure; }; ISL_ARG_DECL(isl_options, struct isl_options, isl_options_arg) diff --git a/isl_options.c b/isl_options.c index dfebd3e..6047535 100644 --- a/isl_options.c +++ b/isl_options.c @@ -51,6 +51,12 @@ struct isl_arg_choice isl_gbr_choice[] = { {0} }; +struct isl_arg_choice isl_closure_choice[] = { + {"isl", ISL_CLOSURE_ISL}, + {"omega", ISL_CLOSURE_OMEGA}, + {0} +}; + struct isl_arg isl_options_arg[] = { ISL_ARG_CHOICE(struct isl_options, lp_solver, 0, "lp-solver", \ isl_lp_solver_choice, ISL_LP_TAB) @@ -62,6 +68,8 @@ ISL_ARG_CHOICE(struct isl_options, context, 0, "context", \ isl_pip_context_choice, ISL_CONTEXT_GBR) ISL_ARG_CHOICE(struct isl_options, gbr, 0, "gbr", \ isl_gbr_choice, ISL_GBR_ONCE) +ISL_ARG_CHOICE(struct isl_options, closure, 0, "closure", \ + isl_closure_choice, ISL_CLOSURE_ISL) ISL_ARG_BOOL(struct isl_options, gbr_only_first, 0, "gbr-only-first", 0) ISL_ARG_END }; diff --git a/isl_transitive_closure.c b/isl_transitive_closure.c index 66cb2ea..9e16184 100644 --- a/isl_transitive_closure.c +++ b/isl_transitive_closure.c @@ -11,6 +11,7 @@ #include "isl_map.h" #include "isl_map_private.h" #include "isl_seq.h" +#include /* Given a map that represents a path with the length of the path * encoded as the difference between the last output coordindate @@ -1137,6 +1138,411 @@ __isl_give isl_map *isl_map_power(__isl_take isl_map *map, unsigned param, return map_power(map, param, exact, 0); } +/* Check whether equality i of bset is a pure stride constraint + * on a single dimensions, i.e., of the form + * + * v = k e + * + * with k a constant and e an existentially quantified variable. + */ +static int is_eq_stride(__isl_keep isl_basic_set *bset, int i) +{ + int k; + unsigned nparam; + unsigned d; + unsigned n_div; + int pos1; + int pos2; + + if (!bset) + return -1; + + if (!isl_int_is_zero(bset->eq[i][0])) + return 0; + + nparam = isl_basic_set_dim(bset, isl_dim_param); + d = isl_basic_set_dim(bset, isl_dim_set); + n_div = isl_basic_set_dim(bset, isl_dim_div); + + if (isl_seq_first_non_zero(bset->eq[i] + 1, nparam) != -1) + return 0; + pos1 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam, d); + if (pos1 == -1) + return 0; + if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + pos1 + 1, + d - pos1 - 1) != -1) + return 0; + + pos2 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d, n_div); + if (pos2 == -1) + return 0; + if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d + pos2 + 1, + n_div - pos2 - 1) != -1) + return 0; + if (!isl_int_is_one(bset->eq[i][1 + nparam + pos1]) && + !isl_int_is_negone(bset->eq[i][1 + nparam + pos1])) + return 0; + + return 1; +} + +/* Given a map, compute the smallest superset of this map that is of the form + * + * { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p } + * + * (where p ranges over the (non-parametric) dimensions), + * compute the transitive closure of this map, i.e., + * + * { i -> j : exists k > 0: + * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } + * + * and intersect domain and range of this transitive closure with + * domain and range of the original map. + * + * If with_id is set, then try to include as much of the identity mapping + * as possible, by computing + * + * { i -> j : exists k >= 0: + * k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p } + * + * instead (i.e., allow k = 0) and by intersecting domain and range + * with the union of the domain and the range of the original map. + * + * In practice, we compute the difference set + * + * delta = { j - i | i -> j in map }, + * + * look for stride constraint on the individual dimensions and compute + * (constant) lower and upper bounds for each individual dimension, + * adding a constraint for each bound not equal to infinity. + */ +static __isl_give isl_map *box_closure(__isl_take isl_map *map, int with_id) +{ + int i; + int k; + unsigned d; + unsigned nparam; + unsigned total; + isl_dim *dim; + isl_set *delta; + isl_set *domain = NULL; + isl_set *range = NULL; + isl_map *app = NULL; + isl_basic_set *aff = NULL; + isl_basic_map *bmap = NULL; + isl_vec *obj = NULL; + isl_int opt; + + isl_int_init(opt); + + delta = isl_map_deltas(isl_map_copy(map)); + + aff = isl_set_affine_hull(isl_set_copy(delta)); + if (!aff) + goto error; + dim = isl_map_get_dim(map); + d = isl_dim_size(dim, isl_dim_in); + nparam = isl_dim_size(dim, isl_dim_param); + total = isl_dim_total(dim); + bmap = isl_basic_map_alloc_dim(dim, + aff->n_div + 1, aff->n_div, 2 * d + 1); + for (i = 0; i < aff->n_div + 1; ++i) { + k = isl_basic_map_alloc_div(bmap); + if (k < 0) + goto error; + isl_int_set_si(bmap->div[k][0], 0); + } + for (i = 0; i < aff->n_eq; ++i) { + if (!is_eq_stride(aff, i)) + continue; + k = isl_basic_map_alloc_equality(bmap); + if (k < 0) + goto error; + isl_seq_clr(bmap->eq[k], 1 + nparam); + isl_seq_cpy(bmap->eq[k] + 1 + nparam + d, + aff->eq[i] + 1 + nparam, d); + isl_seq_neg(bmap->eq[k] + 1 + nparam, + aff->eq[i] + 1 + nparam, d); + isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d, + aff->eq[i] + 1 + nparam + d, aff->n_div); + isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0); + } + obj = isl_vec_alloc(map->ctx, 1 + nparam + d); + if (!obj) + goto error; + isl_seq_clr(obj->el, 1 + nparam + d); + for (i = 0; i < d; ++ i) { + enum isl_lp_result res; + + isl_int_set_si(obj->el[1 + nparam + i], 1); + + res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt, + NULL, NULL); + if (res == isl_lp_error) + goto error; + if (res == isl_lp_ok) { + k = isl_basic_map_alloc_inequality(bmap); + if (k < 0) + goto error; + isl_seq_clr(bmap->ineq[k], + 1 + nparam + 2 * d + bmap->n_div); + isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1); + isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1); + isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt); + } + + res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt, + NULL, NULL); + if (res == isl_lp_error) + goto error; + if (res == isl_lp_ok) { + k = isl_basic_map_alloc_inequality(bmap); + if (k < 0) + goto error; + isl_seq_clr(bmap->ineq[k], + 1 + nparam + 2 * d + bmap->n_div); + isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1); + isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1); + isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt); + } + + isl_int_set_si(obj->el[1 + nparam + i], 0); + } + k = isl_basic_map_alloc_inequality(bmap); + if (k < 0) + goto error; + isl_seq_clr(bmap->ineq[k], + 1 + nparam + 2 * d + bmap->n_div); + if (!with_id) + isl_int_set_si(bmap->ineq[k][0], -1); + isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1); + + domain = isl_map_domain(isl_map_copy(map)); + domain = isl_set_coalesce(domain); + range = isl_map_range(isl_map_copy(map)); + range = isl_set_coalesce(range); + if (with_id) { + domain = isl_set_union(domain, range); + domain = isl_set_coalesce(domain); + range = isl_set_copy(domain); + } + app = isl_map_from_domain_and_range(domain, range); + + isl_vec_free(obj); + isl_basic_set_free(aff); + isl_map_free(map); + bmap = isl_basic_map_finalize(bmap); + isl_set_free(delta); + isl_int_clear(opt); + + map = isl_map_from_basic_map(bmap); + map = isl_map_intersect(map, app); + + return map; +error: + isl_vec_free(obj); + isl_basic_map_free(bmap); + isl_basic_set_free(aff); + isl_map_free(map); + isl_set_free(delta); + isl_int_clear(opt); + return NULL; +} + +/* Check whether app is the transitive closure of map. + * In particular, check that app is acyclic and, if so, + * check that + * + * app \subset (map \cup (map \circ app)) + */ +static int check_exactness_omega(__isl_keep isl_map *map, + __isl_keep isl_map *app) +{ + isl_set *delta; + int i; + int is_empty, is_exact; + unsigned d; + isl_map *test; + + delta = isl_map_deltas(isl_map_copy(app)); + d = isl_set_dim(delta, isl_dim_set); + for (i = 0; i < d; ++i) + delta = isl_set_fix_si(delta, isl_dim_set, i, 0); + is_empty = isl_set_is_empty(delta); + isl_set_free(delta); + if (is_empty < 0) + return -1; + if (!is_empty) + return 0; + + test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map)); + test = isl_map_union(test, isl_map_copy(map)); + is_exact = isl_map_is_subset(app, test); + isl_map_free(test); + + return is_exact; +} + +/* Check if basic map M_i can be combined with all the other + * basic maps such that + * + * (\cup_j M_j)^+ + * + * can be computed as + * + * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+ + * + * In particular, check if we can compute a compact representation + * of + * + * M_i^* \circ M_j \circ M_i^* + * + * for each j != i. + * Let M_i^? be an extension of M_i^+ that allows paths + * of length zero, i.e., the result of box_closure(., 1). + * The criterion, as proposed by Kelly et al., is that + * id = M_i^? - M_i^+ can be represented as a basic map + * and that + * + * id \circ M_j \circ id = M_j + * + * for each j != i. + * + * If this function returns 1, then tc and qc are set to + * M_i^+ and M_i^?, respectively. + */ +static int can_be_split_off(__isl_keep isl_map *map, int i, + __isl_give isl_map **tc, __isl_give isl_map **qc) +{ + isl_map *map_i, *id; + int j = -1; + + map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i])); + *tc = box_closure(isl_map_copy(map_i), 0); + *qc = box_closure(map_i, 1); + id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc)); + + if (!id || !*qc) + goto error; + if (id->n != 1 || (*qc)->n != 1) + goto done; + + for (j = 0; j < map->n; ++j) { + isl_map *map_j, *test; + int is_ok; + + if (i == j) + continue; + map_j = isl_map_from_basic_map( + isl_basic_map_copy(map->p[j])); + test = isl_map_apply_range(isl_map_copy(id), + isl_map_copy(map_j)); + test = isl_map_apply_range(test, isl_map_copy(id)); + is_ok = isl_map_is_equal(test, map_j); + isl_map_free(map_j); + isl_map_free(test); + if (is_ok < 0) + goto error; + if (!is_ok) + break; + } + +done: + isl_map_free(id); + if (j == map->n) + return 1; + + isl_map_free(*qc); + isl_map_free(*tc); + *qc = NULL; + *tc = NULL; + + return 0; +error: + isl_map_free(id); + isl_map_free(*qc); + isl_map_free(*tc); + *qc = NULL; + *tc = NULL; + return -1; +} + +/* Compute an overapproximation of the transitive closure of "map" + * using a variation of the algorithm from + * "Transitive Closure of Infinite Graphs and its Applications" + * by Kelly et al. + * + * We first check whether we can can split of any basic map M_i and + * compute + * + * (\cup_j M_j)^+ + * + * as + * + * M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+ + * + * using a recursive call on the remaining map. + * + * If not, we simply call box_closure on the whole map. + */ +static __isl_give isl_map *compute_closure_omega(__isl_take isl_map *map) +{ + int i, j; + + if (!map) + return NULL; + if (map->n == 1) + return box_closure(map, 0); + + map = isl_map_cow(map); + if (!map) + goto error; + + for (i = 0; i < map->n; ++i) { + int ok; + isl_map *qc, *tc; + ok = can_be_split_off(map, i, &tc, &qc); + if (ok < 0) + goto error; + if (!ok) + continue; + + isl_basic_map_free(map->p[i]); + if (i != map->n - 1) + map->p[i] = map->p[map->n - 1]; + map->n--; + + map = isl_map_apply_range(isl_map_copy(qc), map); + map = isl_map_apply_range(map, qc); + + return isl_map_union(tc, compute_closure_omega(map)); + } + + return box_closure(map, 0); +error: + isl_map_free(map); + return NULL; +} + +/* Compute an overapproximation of the transitive closure of "map" + * using a variation of the algorithm from + * "Transitive Closure of Infinite Graphs and its Applications" + * by Kelly et al. and check whether the result is definitely exact. + */ +static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map, + int *exact) +{ + isl_map *app; + + app = compute_closure_omega(isl_map_copy(map)); + + if (exact) + *exact = check_exactness_omega(map, app); + + isl_map_free(map); + return app; +} + /* 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 @@ -1150,6 +1556,9 @@ __isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map, if (!map) goto error; + if (map->ctx->opt->closure == ISL_CLOSURE_OMEGA) + return transitive_closure_omega(map, exact); + param = isl_map_dim(map, isl_dim_param); map = isl_map_add(map, isl_dim_param, 1); map = map_power(map, param, exact, 1);