mem.cc: directly include required headers

[ppn.git] / da.cc

#include <set>
#include <vector>
#include <iostream>

#include <isa/yaml.h>
#include <isa/pdg.h>
#include "da.h"

extern "C" {
#include "isl_util.h"
}
#include <isl/id.h>
#include <isl/space.h>
#include <isl/aff.h>
#include <isl/set.h>
#include <isl/map.h>
#include <isl/union_set.h>
#include <isl/union_map.h>
#include <isl/flow.h>

using pdg::PDG;
using namespace std;

namespace da {

/* A pair of a node and an access from that node.
 * "map" is the converted access relation.
 * "projected_map" represents the converted access relation without
 * any embedded access relation or access filters
 */
struct na_pair {
    pdg::node   *node;
    pdg::access *access;
    isl_map     *map;
    isl_map     *projected_map;
    void project_out_access_filters(void);
    na_pair(pdg::node *n, pdg::access *a) : node(n), access(a), map(NULL),
        projected_map(NULL) {}
    ~na_pair() {
        isl_map_free(map);
        isl_map_free(projected_map);
    }
};

/* If the access has any embedded filters, then project them out
 * from "projected_map", initializing "projected_map" from "map"
 * if there is no "projected_map" yet.
 */
void na_pair::project_out_access_filters(void)
{
        isl_space *space;
        isl_map *proj;

        if (access->nested.size() == 0)
                return;

        if (!projected_map)
                projected_map = isl_map_copy(map);

        space = isl_space_domain(isl_map_get_space(projected_map));
        space = isl_space_unwrap(space);
        proj = isl_map_domain_map(isl_map_universe(space));

        projected_map = isl_map_apply_domain(projected_map, proj);
}

static int precedes_level_nodes(na_pair *first, na_pair *second);
static int precedes_level_accesses(na_pair *first, na_pair *second);

/* Given a map from a domain to an orthogonal projection of an array
 * (say, the rows of an array), mapping i to m(i), this function
 * extends the range of the mapping to the original array and extends
 * the domain of the mapping correspondingly such that (i,j) maps
 * to (m(i),j), with (m(i),j) identifying an element of the array.
 * The bounds on j are taken from the size of the array.
 *
 * The mapping from i to (i,j) is stored in the "extension" field
 * of the access.
 *
 * The dependences computed using these extended access mappings,
 * will map a (possibly) extended source domain to a (possibly)
 * extended sink domain.  One or both of these domains need to
 * be transformed back to the original domains using the inverse
 * of the corresponding extensions.
 */
static isl_map *extend_access(isl_map *map, na_pair *na)
{
        pdg::array *array = na->access->array;
        unsigned n_out = isl_map_dim(map, isl_dim_out);
        assert(n_out < array->dims.size());
        unsigned s_dim = array->dims.size() - n_out;
        isl_id *array_id = NULL;
        if (isl_map_has_tuple_id(map, isl_dim_out))
                array_id = isl_map_get_tuple_id(map, isl_dim_out);
        isl_space *space = isl_map_get_space(map);
        space = isl_space_drop_dims(space, isl_dim_in,
                                0, isl_space_dim(space, isl_dim_in));
        space = isl_space_drop_dims(space, isl_dim_out,
                                0, isl_space_dim(space, isl_dim_out));
        space = isl_space_add_dims(space, isl_dim_in, s_dim);
        space = isl_space_add_dims(space, isl_dim_out, s_dim);
        isl_basic_map *id = isl_basic_map_identity(space);
        for (int i = 0; i < s_dim; ++i) {
                int v;

                id = isl_basic_map_lower_bound_si(id, isl_dim_out, i, 0);
                v = array->dims[n_out + i] - 1;
                id = isl_basic_map_upper_bound_si(id, isl_dim_out, i, v);
        }
        map = isl_map_product(map, isl_map_from_basic_map(id));
        map = isl_map_flatten_range(map);
        if (array_id)
                map = isl_map_set_tuple_id(map, isl_dim_out, array_id);
        if (!na->access->extension) {
                isl_map *ext = isl_map_copy(map);
                ext = isl_set_unwrap(isl_map_domain(ext));
                ext = isl_map_reverse(isl_map_domain_map(ext));
                na->access->extension = new pdg::IslMap(ext);
        }
        if (!na->access->extended_map)
                na->access->extended_map = new pdg::IslMap(isl_map_copy(map));
        return map;
}

/* If access "access" contains any nested accesses, then the domain
 * of the access relation contains extra dimensions corresponding to
 * the values of the nested accesses.
 * Add these extra dimensions, with ranges given by the value_bounds
 * of the corresponding array to domain "dom".
 * If a nested access array does not have value_bounds, then we assume
 * an infinite interval.
 */
static isl_set *append_nested_value_domains(isl_set *dom, pdg::access *access)
{
        isl_ctx *ctx;

        ctx = isl_set_get_ctx(dom);
        for (int i = 0; i < access->nested.size(); ++i) {
                pdg::call_or_access *coa = access->nested[i];
                assert(coa->type == pdg::call_or_access::t_access);
                pdg::access *nested = coa->access;
                isl_set *bounds;
                if (nested->array->value_bounds)
                        bounds = nested->array->value_bounds->get_isl_set(ctx);
                else {
                        isl_space *dim = isl_set_get_space(dom);
                        dim = isl_space_drop_dims(dim, isl_dim_set,
                                            0, isl_space_dim(dim, isl_dim_set));
                        dim = isl_space_add_dims(dim, isl_dim_set, 1);
                        bounds = isl_set_universe(dim);
                }
                dom = isl_set_product(dom, bounds);
        }
        return dom;
}

/* Combine constraints of the "pure" mapping with the constraints
 * on the domain.  If the range of the mapping is of a dimension
 * that is lower than the dimension of the accessed array,
 * we extend the dimension of both domain and range of the mapping
 * with the missing dimension.  The size of domain and range
 * in these dimensions is set to the extent of the array in the
 * corresponding missing dimension.   Each point in the original
 * domain is therefore expanded to a hyperrectangle and each point
 * in this hyperrectangle is mapped onto a single point in the array.
 *
 * If node->source is a wrapped map, then the iteration domain
 * is the domain of this map.
 */
static isl_map *convert_access(na_pair *na)
{
        isl_map *map = na->access->map->get_isl_map();
        isl_set *dom = na->node->source->get_isl_set();

        if (isl_set_is_wrapping(dom)) {
                dom = isl_map_domain(isl_set_unwrap(dom));
                dom = isl_set_coalesce(dom);
        }

        dom = append_nested_value_domains(dom, na->access);
        if (isl_map_dim(map, isl_dim_in) != isl_set_dim(dom, isl_dim_set))
                assert(0);
        map = isl_map_intersect_domain(map, dom);
        if (isl_map_dim(map, isl_dim_out) != na->access->array->dims.size())
                map = extend_access(map, na);
        return map;
}

typedef std::map<na_pair *, isl_map *> na_pair2map;

struct add_dep_info {
        PDG                             *pdg;
        pdg::array                      *a;
        type                            t;
        enum pdg::dependence::type      dtype;
        /* The sink. */
        na_pair                         *read_na_pair;
        /* The comparison routine that was used during
         * the dependence analysis.
         */
        isl_access_level_before         precedes_level;
        /* Cache of memory based dependence relations.
         * The key of the map refers to the write.
         */
        na_pair2map                     mem_dep;
        /* How many loops are shared by the current sink and source?
         * Initialized to -1.
         */
        int                             n_shared;
        /* For each potential source, what's (up to now) the minimal
         * and maximal number of shared loops?
         * If not set, then we don't know yet.
         */
        std::map<na_pair *, int>        min_n_shared;
        std::map<na_pair *, int>        max_n_shared;

        /* Potential sources that are actually used. */
        std::set<na_pair *>             used;

        __isl_give isl_map *get_mem_dep(na_pair *write_na);
        void clear_mem_dep();
        void set_read_na(na_pair *read_na);
        void update_min_n_shared(na_pair *source_na);

        ~add_dep_info();
};

/* Update min_n_shared of "source_na" to the current number of shared loops.
 * The new value is always smaller than or equal to the old value (if any).
 * If max_n_shared hasn't been set yet, then set it as well.
 */
void add_dep_info::update_min_n_shared(na_pair *source_na)
{
        if (max_n_shared.find(source_na) == max_n_shared.end())
                max_n_shared[source_na] = n_shared;
        min_n_shared[source_na] = n_shared;
}

/* Create an id
 *
 *      __last_<stmt>_<access_nr>_valid
 *
 * corresponding to "na", with "na" attached as user pointer.
 */
static __isl_give isl_id *valid_bit_id(isl_ctx *ctx, na_pair *na)
{
        char name[60];

        snprintf(name, sizeof(name), "__last_%s_%d_valid",
                                na->node->name->s.c_str(), na->access->nr);
        return isl_id_alloc(ctx, name, na);
}

/* Create an id
 *
 *      __last_<stmt>_<access_nr>_shared
 *
 * corresponding to "na", with "na" attached as user pointer.
 */
static __isl_give isl_id *create_shared_id(isl_ctx *ctx, na_pair *na)
{
        char name[60];

        snprintf(name, sizeof(name), "__last_%s_%d_shared",
                                na->node->name->s.c_str(), na->access->nr);
        return isl_id_alloc(ctx, name, na);
}

/* Project out all the dimensions of the given type from "map" except "pos".
 */
static __isl_give isl_map *project_on(__isl_take isl_map *map,
        enum isl_dim_type type, unsigned pos)
{
        unsigned n = isl_map_dim(map, type);

        map = isl_map_project_out(map, type, pos + 1, n - (pos + 1));
        map = isl_map_project_out(map, type, 0, pos);

        return map;
}

/* Does output dimension "pos" have a fixed value in terms of the
 * input dimensions (and parameters)?
 */
static int has_fixed_value(__isl_keep isl_map *map, int pos)
{
        int sv;

        map = isl_map_copy(map);
        map = project_on(map, isl_dim_out, pos);
        sv = isl_map_is_single_valued(map);
        isl_map_free(map);

        return sv;
}

/* Return the position of the parameter with the given "id" in "set",
 * adding it if it wasn't there already.
 */
static int find_or_add_param(__isl_keep isl_set **set, __isl_take isl_id *id)
{
        int pos;

        pos = isl_set_find_dim_by_id(*set, isl_dim_param, id);
        if (pos >= 0) {
                isl_id_free(id);
                return pos;
        }

        pos = isl_set_dim(*set, isl_dim_param);
        *set = isl_set_add_dims(*set, isl_dim_param, 1);
        *set = isl_set_set_dim_id(*set, isl_dim_param, pos, id);

        return pos;
}

/* Add parameters to "set" identifying the last iteration of the access
 * identified by "na".
 *
 * In particular, we add a parameter
 *
 *      __last_<stmt>_<access_nr>_shared >= info->n_shared
 *
 * and parameters
 *
 *      __last_<stmt>_<access_nr>_<i> = it_<i>
 *
 * with i ranging over the iterators, starting at info->n_shared,
 * that are affected by the filters,
 * except those that have a fixed value according to the memory based
 * dependence.
 * "na" is attached to the first two parameters, so that it can be recovered
 * in refine_controls().  If the set already references some of these
 * parameters, then we don't add the parameter again, but instead
 * simply add the corresponding constraint.
 */
static __isl_give isl_set *add_parametrization(__isl_take isl_set *set,
        na_pair *na, add_dep_info *info)
{
        int pos;
        isl_ctx *ctx;
        isl_id *id;
        char name[60];
        int depth;
        isl_map *mem;

        depth = na->node->get_filter_depth();
        mem = info->get_mem_dep(na);
        mem = isl_map_reverse(mem);

        ctx = isl_set_get_ctx(set);
        id = create_shared_id(ctx, na);
        pos = find_or_add_param(&set, id);
        set = isl_set_lower_bound_si(set, isl_dim_param, pos, info->n_shared);

        for (int i = info->n_shared; i < depth; ++i) {
                if (has_fixed_value(mem, i))
                        continue;

                snprintf(name, sizeof(name), "__last_%s_%d_%d",
                         na->node->name->s.c_str(), na->access->nr, i);

                id = isl_id_alloc(ctx, name, NULL);
                pos = find_or_add_param(&set, id);

                set = isl_set_equate(set, isl_dim_param, pos, isl_dim_set, i);
        }

        isl_map_free(mem);
        return set;
}

/* Is the i-th parameter of "map" a control, i.e., a parameter
 * introduced by add_parametrization()?
 * In particular, is the parameter of the form __last_*?
 */
static bool is_control(__isl_keep isl_map *map, int i)
{
        const char *name;
        const char *prefix = "__last_";
        size_t prefix_len = strlen(prefix);

        if (!isl_map_has_dim_id(map, isl_dim_param, i))
                return false;
        name = isl_map_get_dim_name(map, isl_dim_param, i);
        return strncmp(name, prefix, prefix_len) == 0;
}

/* Is the i-th parameter of "set" a control, i.e., a parameter
 * introduced by add_parametrization()?
 * In particular, is the parameter of the form __last_*?
 */
static bool is_control(__isl_keep isl_set *set, int i)
{
        const char *name;
        const char *prefix = "__last_";
        size_t prefix_len = strlen(prefix);

        if (!isl_set_has_dim_id(set, isl_dim_param, i))
                return false;
        name = isl_set_get_dim_name(set, isl_dim_param, i);
        return strncmp(name, prefix, prefix_len) == 0;
}

/* Remove all controls that are redundant, i.e., that do not appear
 * in any of the constraints.
 * Set *has_controls to true if there are any controls that are not redundant.
 */
static __isl_give isl_map *remove_redundant_controls(__isl_take isl_map *dep,
        bool *has_controls)
{
        int i;
        int n_param;

        *has_controls = false;

        n_param = isl_map_dim(dep, isl_dim_param);
        for (i = n_param - 1; i >= 0; --i) {
                if (!is_control(dep, i))
                        continue;
                if (isl_map_involves_dims(dep, isl_dim_param, i, 1))
                        *has_controls = true;
                else
                        dep = isl_map_project_out(dep, isl_dim_param, i, 1);
        }

        return dep;
}

/* Rename controls of "dep" from
 *
 *      __last_<source_stmt>_<source_acc_nr>_*
 *
 * to
 *
 *      __last_<source_stmt>_<source_acc_nr>_<sink_stmt>_<sink_acc_nr>_*
 *
 * "na" represents the sink.
 */
static __isl_give isl_map *rename_controls(__isl_take isl_map *dep, na_pair *na)
{
        char buf[100];
        int n_param;
        const char *name;
        const char *underscore;

        n_param = isl_map_dim(dep, isl_dim_param);
        for (int i = 0; i < n_param; ++i) {
                int len;
                if (!is_control(dep, i))
                        continue;
                name = isl_map_get_dim_name(dep, isl_dim_param, i);
                underscore = strrchr(name, '_');
                assert(underscore);
                len = underscore + 1 - name;
                memcpy(buf, name, len);
                snprintf(buf + len, sizeof(buf) - len, "%s_%d_%s",
                        na->node->name->s.c_str(), na->access->nr, name + len);
                dep = isl_map_set_dim_name(dep, isl_dim_param, i, buf);
        }

        return dep;
}

extern "C" {
        static isl_stat extract_dep(__isl_take isl_map *dep, int must,
                void *dep_user, void *user);
}

/* Extract the single dependence relation from the result of
 * dataflow analyis and assign it to *user.
 */
static isl_stat extract_dep(__isl_take isl_map *dep, int must, void *dep_user,
        void *user)
{
        isl_map **dep_p = (isl_map **) user;
        assert(!*dep_p);
        *dep_p = dep;
        return isl_stat_ok;
}

/* Return the memory based dependence relation from write_na
 * to read_na_pair.  If the "projected_map"
 * fields are not NULL, then use the "projected_map"
 * instead of the "map" of write_na and this->read_na_pair.
 */
__isl_give isl_map *add_dep_info::get_mem_dep(na_pair *write_na)
{
        isl_access_info *acc;
        isl_flow *deps;
        isl_map *read_map, *write_map;
        isl_map *dep = NULL;

        if (mem_dep.find(write_na) != mem_dep.end())
                return isl_map_copy(mem_dep[write_na]);

        if (read_na_pair->projected_map)
                read_map = isl_map_copy(read_na_pair->projected_map);
        else
                read_map = isl_map_copy(read_na_pair->map);
        acc = isl_access_info_alloc(read_map, read_na_pair, precedes_level, 1);
        if (write_na->projected_map)
                write_map = isl_map_copy(write_na->projected_map);
        else
                write_map = isl_map_copy(write_na->map);
        acc = isl_access_info_add_source(acc, write_map, 0, write_na);
        deps = isl_access_info_compute_flow(acc);
        isl_flow_foreach(deps, &extract_dep, &dep);
        isl_flow_free(deps);

        mem_dep[write_na] = isl_map_copy(dep);

        return dep;
}

/* Clear the cache of memory based dependence relations.
 */
void add_dep_info::clear_mem_dep()
{
        na_pair2map::iterator it;

        for (it = mem_dep.begin(); it != mem_dep.end(); ++it)
                isl_map_free(it->second);
        mem_dep.clear();
}

/* Set read_na_pair to read_na.
 *
 * If the cache of memory based dependence relations contains any
 * elements then they refer to a different read, so we need to clear
 * the cache.
 *
 * We also clear the set of used potential sources and reset
 * the data that keeps track of the number of shared loops between
 * the sink (read_na_pair) and the sources.
 */
void add_dep_info::set_read_na(na_pair *read_na)
{
        used.clear();
        clear_mem_dep();
        n_shared = -1;
        min_n_shared.clear();
        max_n_shared.clear();
        read_na_pair = read_na;
}

add_dep_info::~add_dep_info()
{
        clear_mem_dep();
}

/* Is the name of parameter "i" of "space" of the form __last_*_suffix?
 */
static bool is_last_with_suffix(__isl_keep isl_space *space, int i,
        const char *suffix, size_t suffix_len)
{
        const char *prefix = "__last_";
        size_t prefix_len = strlen(prefix);
        const char *name;
        size_t len;

        if (!isl_space_has_dim_id(space, isl_dim_param, i))
                return false;
        name = isl_space_get_dim_name(space, isl_dim_param, i);
        if (strncmp(name, prefix, prefix_len))
                return false;
        len = strlen(name);
        return len > suffix_len && !strcmp(name + len - suffix_len, suffix);
}

/* Is the name of parameter "i" of "space" of the form __last_*_valid?
 * In practice, those are the parameters __last_*_valid, created
 * in add_parametrization().
 */
static bool is_valid_bit(__isl_keep isl_space *space, int i)
{
        const char *suffix = "_valid";

        return is_last_with_suffix(space, i, suffix, strlen(suffix));
}

/* Is the name of parameter "i" of "space" of the form __last_*_shared?
 * In practice, those are the parameters __last_*_shared, created
 * in add_parametrization().
 */
static bool is_shared(__isl_keep isl_space *space, int i)
{
        const char *suffix = "_shared";
        size_t suffix_len = strlen(suffix);

        return is_last_with_suffix(space, i, suffix, strlen(suffix));
}

/* Assuming "coa" is a (read) access, return the array being
 * accessed.
 */
static pdg::array *get_filter_array(pdg::call_or_access *coa)
{
        assert(coa->type == pdg::call_or_access::t_access);
        return coa->access->array;
}

/* Compute a map between domain elements (i) of "map1" and range elements
 * of "map2" (j) such that all the images of i in "map1" map to j through
 * "map2" and such that there is at least one such image element.
 *
 * In other words, the result contains those pairs of elements such that
 * map1(i) \cap map2^-1(j) is non-empty and map1(i) \subseteq map2^-1(j).
 *
 * Equivalently, compute
 *
 *      (map1 . map2) \setminus
 *              (map1 . ((\range map1 \to \range map2) \setminus map2))
 *
 * If map1 is single valued, then we can do a simple join.
 */
static __isl_give isl_union_map *join_non_empty_subset(
        __isl_take isl_union_map *umap1, __isl_take isl_union_map *umap2)
{
        isl_union_set *dom, *ran;
        isl_union_map *univ;
        isl_union_map *res;

        if (isl_union_map_is_single_valued(umap1))
                return isl_union_map_apply_range(umap1, umap2);

        res = isl_union_map_apply_range(isl_union_map_copy(umap1),
                                        isl_union_map_copy(umap2));
        dom = isl_union_map_range(isl_union_map_copy(umap1));
        ran = isl_union_map_range(isl_union_map_copy(umap2));
        univ = isl_union_map_from_domain_and_range(dom, ran);
        umap2 = isl_union_map_subtract(univ, umap2);
        umap1 = isl_union_map_apply_range(umap1, umap2);
        res = isl_union_map_subtract(res, umap1);

        return res;
}

/* Compute a map between domain elements (i) of "map1" and range elements
 * of "map2" (j) such that the images of i in "map1" include all those
 * elements that map to j through "map2" and such that there is
 * at least one such image element.
 *
 * In other words, the result contains those pairs of elements such that
 * map1(i) \cap map2^-1(j) is non-empty and map1(i) \supseteq map2^-1(j).
 *
 * Equivalently, compute
 *
 *      (map1 . map2) \setminus
 *              (((\domain map1 \to \domain map2) \setminus map1) . map2)
 *
 * If map1 is single valued, then we can do a simple join.
 */
static __isl_give isl_union_map *join_non_empty_superset(
        __isl_take isl_union_map *umap1, __isl_take isl_union_map *umap2)
{
        isl_union_set *dom, *ran;
        isl_union_map *univ;
        isl_union_map *res;

        if (isl_union_map_is_single_valued(umap2))
                return isl_union_map_apply_range(umap1, umap2);

        res = isl_union_map_apply_range(isl_union_map_copy(umap1),
                                        isl_union_map_copy(umap2));
        dom = isl_union_map_domain(isl_union_map_copy(umap1));
        ran = isl_union_map_domain(isl_union_map_copy(umap2));
        univ = isl_union_map_from_domain_and_range(dom, ran);
        umap1 = isl_union_map_subtract(univ, umap1);
        umap1 = isl_union_map_apply_range(umap1, umap2);
        res = isl_union_map_subtract(res, umap1);

        return res;
}

/* Return those elements in the domain of "umap" where "umap" is multi-valued.
 *
 * In particular, construct a mapping between domain elements of "umap"
 * and pairs of corresponding image elements.
 * Remove pairs of identical image elements from the range of this mapping.
 * The result is a mapping between domain elements and pairs of different
 * corresponding image elements.  The domain of this mapping contains those
 * domain elements of "umap" with at least two images.
 */
static __isl_give isl_union_set *multi_valued(__isl_keep isl_union_map *umap)
{
        isl_union_map *multi, *id;

        multi = isl_union_map_range_product(isl_union_map_copy(umap),
                                            isl_union_map_copy(umap));
        id = isl_union_map_universe(isl_union_map_copy(multi));
        id = isl_union_set_unwrap(isl_union_map_range(id));
        id = isl_union_set_identity(isl_union_map_domain(id));
        multi = isl_union_map_subtract_range(multi, isl_union_map_wrap(id));
        return isl_union_map_domain(multi);
}

/* Given two filter access relations, return a mapping between the domain
 * elements of these access relations such that they access "the same filter".
 * In particular, any pair of elements in the returned relation
 * accesses at least one element in common, but if subset1 is set,
 * then the set of elements accessed by the first is a subset of the
 * set of elements accessed by the second.  Similarly, if subset2 is set,
 * then the set of elements accessed by the second is a subset of the
 * set of elements accessed by the first.  If both are set, then we further
 * impose that both should access exactly one element.
 * "space" is the space in which the result should live.
 * Although "map1" and "map2" are allowed to have ranges in multiple spaces,
 * their domains should live in a single space.  "space" is the space
 * of the relation between those two domains.
 *
 * Call the given maps A and B.
 *
 * A relation between domains elements of A and B that access at least
 * one element in common can be obtained as
 *
 *      A . B^-1
 *
 * To ensure that all elements accessed through A form a subset of
 * the elements accessed through B, we compute join_non_empty_subset(A, B^-1).
 *
 * Ensuring that all elements accessed through B form a subset of
 * the elements accessed through A is handled in a similar way.
 *
 * To remove those iterations that access more that one element,
 * we compute those parts of the domains where A and B are multi-valued
 * and subtract them from domain and range of the result.
 */
static __isl_give isl_map *compute_common_filter(__isl_keep isl_union_map *map1,
        bool subset1, __isl_keep isl_union_map *map2, bool subset2,
        __isl_keep isl_space *space)
{
        isl_union_map *reverse, *common;
        isl_union_set *bad;
        isl_map *res;

        reverse = isl_union_map_reverse(isl_union_map_copy(map2));

        if (subset1 && !subset2) {
                common = join_non_empty_subset(isl_union_map_copy(map1),
                                                reverse);
        } else if (!subset1 && subset2) {
                common = join_non_empty_superset(isl_union_map_copy(map1),
                                                reverse);
        } else {
                common = isl_union_map_apply_range(isl_union_map_copy(map1),
                                                        reverse);
                if (subset1 && subset2) {
                        bad = multi_valued(map1);
                        common = isl_union_map_subtract_domain(common, bad);
                        bad = multi_valued(map2);
                        common = isl_union_map_subtract_range(common, bad);
                }
        }

        res = isl_union_map_extract_map(common, isl_space_copy(space));
        isl_union_map_free(common);
        return res;
}

/* Assuming "coa" is a (read) access, construct a union map from the domain
 * of the access relation to the access relations of the corresponding
 * writes.  If we are unable to determine the corresponding writes, then
 * return a map to the read access relation.
 */
static __isl_give isl_union_map *extract_access_map(pdg::call_or_access *coa)
{
        assert(coa->type == pdg::call_or_access::t_access);
        return coa->access->extract_access_map();
}

/* Return a set that contains all possible filter values,
 * where the possible values for a given filter is either as specified
 * by the value_bounds property of the corresponding array or the universe.
 */
static __isl_give isl_set *compute_filter_bounds(pdg::node *node)
{
        isl_ctx *ctx;
        isl_set *bounds;

        ctx = isl_set_get_ctx(node->source->set);

        bounds = isl_set_universe(isl_space_set_alloc(ctx, 0, 0));
        for (int i = 0; i < node->filters.size(); ++i) {
                isl_set *bnd;
                pdg::array *array;
                pdg::call_or_access *coa = node->filters[i];
                assert(coa->type == pdg::call_or_access::t_access);
                array = coa->access->array;
                if (array->value_bounds)
                        bnd = array->value_bounds->get_isl_set();
                else
                        bnd = isl_set_universe(isl_space_set_alloc(ctx, 0, 1));
                bounds = isl_set_flat_product(bounds, bnd);
        }

        return bounds;
}

/* Return either the filter values themselves or their complement,
 * taken with respect to the bounds on the filter values.
 */
static __isl_give isl_map *compute_filter_values(pdg::node *node,
        bool complement)
{
        isl_set *bounds;
        isl_map *value;
        isl_map *res;

        value = isl_set_unwrap(node->source->get_isl_set());
        if (!complement)
                return value;

        bounds = compute_filter_bounds(node);
        res = isl_map_from_domain_and_range(
                                isl_map_domain(isl_map_copy(value)), bounds);
        res = isl_map_subtract(res, value);

        return res;
}

/* Equate the first "n" input and output dimensions of "map"
 * and return the result.
 */
static __isl_give isl_map *share(__isl_take isl_map *map, int n)
{
        for (int i = 0; i < n; ++i)
                map = isl_map_equate(map, isl_dim_in, i, isl_dim_out, i);

        return map;
}

/* Return the set of source iterations of "na" either at the last
 * iteration (if valid is set) or after the last iteration
 * (if valid is not set).  "id" represents the control variable
 * corresponding to the number of shared loops (__last_<na>_shared).
 *
 * In particlar, if valid is set, we return the set
 *
 *      { S[i] : i = __last_<na> and
 *               __last_<na>_shared >= info->n_shared }
 *
 * If valid is not set, we return the set
 *
 *      { S[i] : (i >> __last_<na> and __last_<na>_shared >= info->n_shared) or
 *               __last_<na>_shared < info->n_shared }
 *
 * That is, the iterations after the last if there is a last iteration
 * with at least info->n_shared shared loops
 * or just any iteration if there is no such last iteration.
 *
 * The lexicographic order i >> __last_<na> is imposed on the loop iterators
 * that are affected by any filters.
 */
static __isl_give isl_set *source_iterations(na_pair *na,
        __isl_keep isl_id *id, bool valid, add_dep_info *info)
{
        isl_space *space;
        isl_set *set, *invalid;
        isl_map *map_after;
        int pos;
        int depth;

        space = isl_set_get_space(na->node->source->set);
        space = isl_space_domain(isl_space_unwrap(space));

        set = isl_set_universe(space);
        set = add_parametrization(set, na, info);

        if (valid)
                return set;

        depth = na->node->get_filter_depth();

        space = isl_space_map_from_set(isl_set_get_space(set));
        map_after = isl_map_lex_lt_first(space, depth);
        map_after = share(map_after, info->n_shared);
        map_after = isl_map_intersect_domain(map_after, set);
        set = isl_map_range(map_after);

        invalid = isl_set_universe(isl_set_get_space(set));
        pos = isl_set_find_dim_by_id(invalid, isl_dim_param, id);
        assert(pos >= 0);
        invalid = isl_set_upper_bound_si(invalid,
                                        isl_dim_param, pos, info->n_shared - 1);
        set = isl_set_union(set, invalid);

        return set;
}

/* Look for a matching between the filters of node1 and those of node2.
 * That is look for pairs of filters of the two nodes that are "the same".
 * Return true if any such matching can be found.  The correspondence between
 * the filters is returned in *same_value_p, while the pairs of iterations
 * where the filters are the same is returned in *same_filter_p.
 * The first "n_shared" dimensions of these iterations are guaranteed
 * to be equal to each other.
 *
 * Two filter accesses are considered "the same" if they access at least
 * one element in common.  Moreover, if valid1 is false then the set
 * of elements accessed by an element from node1 should be a subset
 * of the set of elements accessed by the corresponding element from node2.
 * Similarly for valid2.
 *
 * We perform a greedy search, checking if two filters could possibly
 * match given the matchings we have performed before and updating
 * the matching if it is indeed possible.
 *
 * Note that this function only computes one of the possibly many matchings.
 */
static bool compute_matching(pdg::node *node1, bool valid1,
        pdg::node *node2, bool valid2, __isl_give isl_map **same_filter_p,
        __isl_give isl_map **same_value_p, int n_shared)
{
        bool any = false;
        isl_space *space1, *space2;
        isl_map *same_filter;
        isl_map *same_value;

        space1 = isl_space_unwrap(isl_set_get_space(node1->source->set));
        space2 = isl_space_unwrap(isl_set_get_space(node2->source->set));
        space1 = isl_space_product(space1, space2);
        space2 = isl_space_unwrap(isl_space_range(isl_space_copy(space1)));
        space1 = isl_space_unwrap(isl_space_domain(space1));

        same_filter = isl_map_universe(isl_space_copy(space1));
        same_filter = share(same_filter, n_shared);
        same_value = isl_map_universe(space2);

        for (int i = 0; i < node1->filters.size(); ++i) {
                isl_union_map *map_i;
                pdg::call_or_access *filter_i = node1->filters[i];
                pdg::array *array_i = get_filter_array(filter_i);

                map_i = extract_access_map(node1->filters[i]);

                for (int j = 0; j < node2->filters.size(); ++j) {
                        pdg::call_or_access *filter_j;
                        filter_j = node2->filters[j];
                        pdg::array *array_j = get_filter_array(filter_j);
                        isl_union_map *map_j;
                        isl_map *same_filter_ij;
                        
                        if (array_i != array_j)
                                continue;

                        map_j = extract_access_map(node2->filters[j]);
                        same_filter_ij = compute_common_filter(map_i, !valid1,
                                                                map_j, !valid2,
                                                                space1);
                        same_filter_ij = isl_map_intersect(same_filter_ij,
                                                isl_map_copy(same_filter));
                        if (isl_map_is_empty(same_filter_ij))
                                isl_map_free(same_filter_ij);
                        else {
                                any = true;
                                isl_map_free(same_filter);
                                same_filter = same_filter_ij;
                                same_value = isl_map_equate(same_value,
                                                            isl_dim_in, i,
                                                            isl_dim_out, j);
                        }
                        isl_union_map_free(map_j);
                }

                isl_union_map_free(map_i);
        }
        isl_space_free(space1);

        if (any) {
                *same_value_p = same_value;
                *same_filter_p = same_filter;
        } else {
                isl_map_free(same_value);
                isl_map_free(same_filter);
                *same_value_p = NULL;
                *same_filter_p = NULL;
        }

        return any;
}

/* Given a set of sink iterations "sink", mappings "map1" and "map2"
 * from two potential sources to this sink,
 * the possible filter values "value1" and "value2" at those
 * potential sources, a relation "same_filter" between the two
 * potential sources expressing when some filters of the two
 * potential sources are the same and the correponding matching
 * "same_value" between the filter values,
 * remove those elements from the sink that have
 * corresponding pairs of potential source iterations that should
 * have the same filter values but do not.
 *
 * Let us call the sink S, the potential sources A and B and the
 * corresponding filters F and G.
 *
 * We start from the mappings A[..] -> S[..] and B[..] -> S[..],
 * combine them into
 *
 *      S[..] -> [A[..] -> B[..]]
 *
 * and intersect the range with the condition "same_filter" on A and B,
 * resulting in a mapping from sink iterations to pairs of potential
 * source iterations that should have the same filter values
 * (as specified by "same_value").
 *
 * We subtract from the range of this mapping those pairs of
 * potential source iterations that actually have the same filter values.
 * The result is a mapping from sink iterations to pairs of potential
 * source iterations that should have the same filter values but do not.
 *
 * The mapping between potential source iterations that have the
 * same filter values is obtained by combining the mappings
 * A[..] -> F[..] and B[..] -> G[..] into
 *
 *      [A[..] -> B[..]] -> [F[..] -> G[..]]
 *
 * intersecting the range with "same_value" and then computing the domain.
 */
static __isl_give isl_set *remove_conflict(__isl_take isl_set *sink,
        __isl_take isl_map *map1, __isl_take isl_map *value1,
        __isl_take isl_map *map2, __isl_take isl_map *value2,
        __isl_take isl_map *same_filter, __isl_take isl_map *same_value)
{
        isl_map *value, *conflict;
        isl_set *conflict_set;

        conflict = isl_map_domain_product(map1, map2);
        conflict = isl_map_reverse(conflict);

        conflict = isl_map_intersect_range(conflict, isl_map_wrap(same_filter));

        value = isl_map_product(value1, value2);
        value = isl_map_intersect_range(value, isl_map_wrap(same_value));
        conflict = isl_map_subtract_range(conflict, isl_map_domain(value));

        conflict_set = isl_map_domain(conflict);

        sink = isl_set_subtract(sink, conflict_set);

        return sink;
}

/* Remove inconsistencies from the set of sink iterations "sink"
 * based on two potential sources identified by "id1" and "id2"
 * (representing the number of shared loops),
 * in particular, on either the last iteration where the filters hold
 * (if valid? is set) or on later iterations (if valid? is not set).
 *
 * Let us first consider the case where both "valid1" and "valid2" are set.
 * If the last iterations of the corresponding sources access the same
 * filters, then these filters should have the same value.
 * If a filter access accesses more than one element, then these elements
 * should all have the same value.  It is therefore sufficient for the
 * two last iterations to access at least one element in common for there
 * to be a requirement that the corresponding values should be the same.
 * We therefore obtain the filter values, the mappings from the sink
 * to the last iterations, a matching between the
 * the filters of the corresponding sources and remove conflicts from "dep".
 *
 * If one or both of the valid bits are not set, then we need to make
 * some changes.  First the inconsistencies now do not arise from
 * the filter values at the last iteration, but from the filter values
 * lying _outside_ of the possible values for all iterations _after_
 * the "last" (i.e., the last iteration satisfying the filter constraints).
 * In case there is no last iteration with at least info->n_shared shared loops,
 * then the filter values should lie outside of the possible values
 * for any potential source iteration with info->n_shared shared loops.
 * Note however, that if the filter access relation accesses several
 * elements, then it is sufficient for one of those to have a value
 * outside of the possible values.  We can therefore only consider
 * any inconsistencies for those cases where the set of accessed elements
 * forms a subset of the set of accessed elements through the other potential
 * source.  If valid1 is not set, but valid2 is set, then we consider
 * those pairs of potential source iterations where the first accesses
 * a subset of the second and we impose that at least one of those
 * accessed elements has a valid outside the possible values.
 * Since those accessed elements form a subset of the elements accessed
 * by the other potential source, there is at least one element that
 * has a value outside of the posssible values on the first potential source
 * and a value belonging to the posssible values on the second potential source.
 * We can therefore impose that this value should exist.
 *
 * If both valid1 and valid2 are not set, then we can only
 * impose a constraint on those pairs of iterations that access the same
 * single element.  We then know that the value of this single element
 * accessed by both potential sources should lie outside of the possible
 * values on both sides.
 */
static __isl_give isl_set *remove_inconsistencies(__isl_take isl_set *sink,
        add_dep_info *info, __isl_keep isl_id *id1, bool valid1,
        __isl_keep isl_id *id2, bool valid2)
{
        na_pair *write1_na, *write2_na;
        isl_set *source;
        isl_map *value1, *value2;
        isl_map *same_filter;
        isl_map *same_value;
        isl_map *mem1, *mem2;

        write1_na = (na_pair *) isl_id_get_user(id1);
        write2_na = (na_pair *) isl_id_get_user(id2);

        if (!compute_matching(write1_na->node, valid1, write2_na->node, valid2,
                                &same_filter, &same_value, info->n_shared))
                return sink;

        value1 = compute_filter_values(write1_na->node, !valid1);
        value2 = compute_filter_values(write2_na->node, !valid2);

        mem1 = info->get_mem_dep(write1_na);
        source = source_iterations(write1_na, id1, valid1, info);
        mem1 = share(mem1, info->n_shared);
        mem1 = isl_map_intersect_domain(mem1, source);

        mem2 = info->get_mem_dep(write2_na);
        source = source_iterations(write2_na, id2, valid2, info);
        mem2 = share(mem2, info->n_shared);
        mem2 = isl_map_intersect_domain(mem2, source);

        sink = remove_conflict(sink, mem1, value1, mem2, value2,
                                same_filter, same_value);

        return sink;
}

/* Remove inconsistencies from the set of sink iterations "sink"
 * based on two potential sources identified by "id1" and "id2",
 * (representing the number of shared loops).
 */
static __isl_give isl_set *remove_inconsistencies(__isl_take isl_set *sink,
        add_dep_info *info, __isl_keep isl_id *id1, __isl_keep isl_id *id2)
{
        sink = remove_inconsistencies(sink, info, id1, false, id2, false);
        sink = remove_inconsistencies(sink, info, id1, false, id2, true);
        sink = remove_inconsistencies(sink, info, id1, true, id2, false);
        sink = remove_inconsistencies(sink, info, id1, true, id2, true);
        return sink;
}

/* Remove inconsistencies from the set of sink iterations "sink"
 * based on the potential source identified by "id"
 * (representing the number of shared loops),
 * in particular, on either the last iteration where the filters hold
 * (if valid is set) or on later iterations (if valid is not set).
 *
 * This function is very similar to the remove_inconsistencies
 * function above that considers two potential sources instead
 * of the sink and one potential source.  The main differences
 * are that for the sink, the filters always hold and that the mapping
 * from sink iterations to sink iterations is computed in a different
 * (and fairly trivial) way.
 */
static __isl_give isl_set *remove_inconsistencies(__isl_take isl_set *sink,
        add_dep_info *info, __isl_keep isl_id *id, bool valid)
{
        na_pair *read_na, *write_na;
        isl_set *after;
        isl_map *value1, *value2;
        isl_map *same_filter;
        isl_map *same_value;
        isl_map *id_map, *mem;

        read_na = info->read_na_pair;
        write_na = (na_pair *) isl_id_get_user(id);

        if (!compute_matching(read_na->node, true, write_na->node, valid,
                                &same_filter, &same_value, info->n_shared))
                return sink;

        value1 = compute_filter_values(read_na->node, false);
        value2 = compute_filter_values(write_na->node, !valid);

        id_map = isl_set_identity(isl_set_copy(sink));

        mem = info->get_mem_dep(write_na);
        after = source_iterations(write_na, id, valid, info);
        mem = share(mem, info->n_shared);
        mem = isl_map_intersect_domain(mem, after);

        sink = remove_conflict(sink, id_map, value1, mem, value2,
                                same_filter, same_value);

        return sink;
}

/* Remove inconsistencies from the set of sink iterations "sink"
 * based on the potential source identified by "id"
 * (representing the number of shared loops).
 */
static __isl_give isl_set *remove_inconsistencies(__isl_take isl_set *sink,
        add_dep_info *info, __isl_keep isl_id *id)
{
        sink = remove_inconsistencies(sink, info, id, false);
        sink = remove_inconsistencies(sink, info, id, true);
        return sink;
}

/* Remove all parameters that were introduced by add_parametrization().
 */
static __isl_give isl_map *remove_all_controls(__isl_take isl_map *dep)
{
        int i;
        int n_param;
        const char *name;

        n_param = isl_map_dim(dep, isl_dim_param);
        for (i = n_param - 1; i >= 0; --i) {
                if (!is_control(dep, i))
                        continue;
                dep = isl_map_project_out(dep, isl_dim_param, i, 1);
        }
        dep = isl_map_coalesce(dep);

        return dep;
}

/* Remove all parameters that were introduced by add_parametrization().
 */
static __isl_give isl_set *remove_all_controls(__isl_take isl_set *set)
{
        int i;
        int n_param;
        const char *name;

        n_param = isl_set_dim(set, isl_dim_param);
        for (i = n_param - 1; i >= 0; --i) {
                if (!is_control(set, i))
                        continue;
                set = isl_set_project_out(set, isl_dim_param, i, 1);
        }
        set = isl_set_coalesce(set);

        return set;
}

/* Create the map
 *
 *      { a -> b : (shared >= max and
 *                      the first max iterators are equal) or
 *                (shared = max - 1 and
 *                      the first max - 1 iterators are equal and
 *                      dimension max - 1 of a is smaller than that of b) or
 *                ...
 *                (shared = min and
 *                      the first min iterators are equal and
 *                      dimension min of a is smaller than that of b) }
 */
static __isl_give isl_map *compute_shared_map(__isl_take isl_space *space,
        __isl_keep isl_id *shared_id, int min, int max)
{
        isl_map *shared_map;
        int shared_pos;

        shared_pos = isl_space_find_dim_by_id(space, isl_dim_param, shared_id);
        shared_map = isl_map_universe(isl_space_copy(space));
        shared_map = isl_map_lower_bound_si(shared_map, isl_dim_param,
                                            shared_pos, max);
        shared_map = share(shared_map, max);

        for (int i = min; i < max; ++i) {
                isl_map *shared_map_i;

                shared_map_i = isl_map_universe(isl_space_copy(space));
                shared_map_i = isl_map_fix_si(shared_map_i, isl_dim_param,
                                                shared_pos, i);
                shared_map_i = share(shared_map_i, i);
                shared_map_i = isl_map_order_lt(shared_map_i, isl_dim_in, i,
                                                isl_dim_out, i);

                shared_map = isl_map_union(shared_map, shared_map_i);
        }

        isl_space_free(space);

        return shared_map;
}

/* Different parts of the final dependence relations may have been
 * created at different depths and may therefore have a different
 * number of dimensions of the last iterator.  The __last_<na>_shared
 * value determines how many of the dimensions are implicitly equal
 * to those of the sink iteration.  This function creates a set that
 * makes these equalities explicit, so that we can later remove
 * the __last_<na>_shared parameter.  It also marks those parts
 * that have a number of shared iterators that is smaller than the minimum
 * as not having any last iteration.
 *
 * In particular, we create a set in the space of the sink of the form
 *
 *      { s : __last_<na>_valid = 0 or
 *              (__last_<na>_valid = 1 and
 *               __last_<na>_<i> is a potential source iteration and
 *               ((__last_<na>_shared >= max_shared and
 *                      the first max_shared iterators are equal) or
 *                (__last_<na>_shared = max_shared - 1 and
 *                      the first max_shared - 1 iterators are equal and
 *                      iterator max_shared - 1 of the source is smaller) or
 *                ...
 *                (__last_<na>_shared = min_shared and
 *                      the first min_shared iterators are equal and
 *                      iterator min_shared of the source is smaller))) }
 *
 * That is, for those parts with __last_<na>_shared smaller than
 * max_n_shared[source_na], intersection with the set will introduce
 * __last_<na>_<i> parameters (assuming they don't have a known fixed value)
 * up until __last_<na>_shared and equate them to the corresponding iterators
 * of the sink.
 */
static __isl_give isl_set *shared_refinement(__isl_keep isl_id *shared_id,
        add_dep_info *info)
{
        isl_ctx *ctx;
        isl_space *space;
        isl_map *mem;
        na_pair *source_na;
        isl_map *shared_map;
        int valid_pos, shared_pos;
        isl_set *valid;
        isl_set *invalid;
        isl_set *domain;
        isl_id *valid_id;

        ctx = isl_id_get_ctx(shared_id);

        source_na = (na_pair *) isl_id_get_user(shared_id);
        valid_id = valid_bit_id(ctx, source_na);

        mem = info->get_mem_dep(source_na);
        space = isl_map_get_space(mem);
        isl_map_free(mem);

        shared_pos = isl_space_dim(space, isl_dim_param);
        space = isl_space_add_dims(space, isl_dim_param, 1);
        space = isl_space_set_dim_id(space, isl_dim_param,
                                        shared_pos, isl_id_copy(shared_id));

        shared_map = compute_shared_map(isl_space_copy(space), shared_id,
                info->min_n_shared[source_na], info->max_n_shared[source_na]);

        domain = isl_set_universe(isl_space_domain(space));
        valid_pos = find_or_add_param(&domain, isl_id_copy(valid_id));
        domain = isl_set_fix_si(domain, isl_dim_param, valid_pos, 1);
        domain = add_parametrization(domain, source_na, info);
        shared_map = isl_map_intersect_domain(shared_map, domain);

        valid = isl_map_range(shared_map);
        invalid = isl_set_universe(isl_set_get_space(valid));
        shared_pos = isl_set_find_dim_by_id(invalid, isl_dim_param, shared_id);
        invalid = isl_set_upper_bound_si(invalid, isl_dim_param, shared_pos,
                                        info->n_shared - 1);

        valid_pos = find_or_add_param(&invalid, valid_id);
        invalid = isl_set_fix_si(invalid, isl_dim_param, valid_pos, 0);

        return isl_set_union(valid, invalid);
}

/* Different parts of the final dependence relations may have been
 * created at different depths and may therefore have a different
 * number of dimensions of the last iterator.  The __last_*_shared
 * value determines how many of the dimensions are implicitly equal
 * to those of the sink iteration.
 *
 * For each of the __last_*_shared parameters, explicitly add
 * the implicitly equal __last_*_i iterators by intersecting
 * the sink with the set computed by shared_refinement.
 * Finally, remove the __last_*_shared parameters.
 */
static __isl_give isl_map *refine_shared(__isl_take isl_map *dep,
        add_dep_info *info)
{
        isl_space *space;
        isl_set *refinement;
        int n_param;

        space = isl_map_get_space(dep);
        n_param = isl_space_dim(space, isl_dim_param);

        refinement = isl_set_universe(isl_space_range(isl_space_copy(space)));

        for (int i = 0; i < n_param; ++i) {
                isl_set *ref_i;
                isl_id *id;

                if (!is_shared(space, i))
                        continue;
                if (!isl_map_involves_dims(dep, isl_dim_param, i, 1))
                        continue;

                id = isl_space_get_dim_id(space, isl_dim_param, i);
                ref_i = shared_refinement(id, info);
                isl_id_free(id);

                if (isl_set_is_wrapping(refinement)) {
                        isl_map *map = isl_set_unwrap(refinement);
                        map = isl_map_intersect_domain(map, ref_i);
                        refinement = isl_map_wrap(map);
                } else
                        refinement = isl_set_intersect(refinement, ref_i);
        }
        dep = isl_map_intersect_range(dep, refinement);

        isl_space_free(space);

        space = isl_map_get_space(dep);
        n_param = isl_space_dim(space, isl_dim_param);
        for (int i = n_param - 1; i >= 0; --i) {
                if (!is_shared(space, i))
                        continue;
                dep = isl_map_project_out(dep, isl_dim_param, i, 1);
        }
        isl_space_free(space);

        dep = isl_map_coalesce(dep);

        return dep;
}

/* Compute the gist of "dep" with respect to the fact that
 *
 *      0 <= __last_*_valid <= 1
 */
static __isl_give isl_map *gist_valid(__isl_take isl_map *dep)
{
        isl_space *space;
        isl_set *valid;
        int n_param;

        space = isl_space_params(isl_map_get_space(dep));
        n_param = isl_space_dim(space, isl_dim_param);

        valid = isl_set_universe(isl_space_copy(space));
        for (int i = 0; i < n_param; ++i) {
                if (!is_valid_bit(space, i))
                        continue;

                valid = isl_set_lower_bound_si(valid, isl_dim_param, i, 0);
                valid = isl_set_upper_bound_si(valid, isl_dim_param, i, 1);
        }

        isl_space_free(space);

        dep = isl_map_gist_params(dep, valid);

        return dep;
}

/* Simplify the constraints on the parameters introduced
 * in add_parametrization().
 * We first add __last_*_i iterators that are only implicitly referred
 * to through the __last_*_shared parameters.
 * Then, we remove all those parameters that turn out not be needed.
 * If there are any of those parameters left, then we compute the gist
 * with respect to the valid bit being either 0 or 1 and rename
 * the parameters to also include a reference to the sink.
 * The resulting relation is assigned to controlled_relation,
 * while the relation field is assigned the result of projecting out
 * all those parameters.
 */
static __isl_give isl_map *simplify_controls(__isl_take isl_map *dep,
        add_dep_info *info, bool *has_controls)
{
        bool hs;

        if (!has_controls)
                has_controls = &hs;
        dep = refine_shared(dep, info);
        dep = remove_redundant_controls(dep, has_controls);
        if (*has_controls) {
                dep = gist_valid(dep);
                dep = rename_controls(dep, info->read_na_pair);
        }

        return dep;
}

/* Extract a single dependence from the result of dataflow analysis.
 *
 * We first simplify the constraints on the parameters introduced
 * in add_parametrization().
 *
 * If the dependence relation turns out to be empty, we simply return.
 * Otherwise, we create a corresponding pdg::dependence and keep track
 * of the fact that the potential source is actually used
 * so that we can remove any reference to potential sources that are
 * never used from the dependence relations.
 */
static isl_stat add_dep(__isl_take isl_map *dep, int must, void *dep_user,
        void *user)
{
        bool has_controls;
        na_pair *write_na = (na_pair *) dep_user;
        add_dep_info *info = (struct add_dep_info *)user;
        isl_ctx *ctx = info->pdg->get_isl_ctx();
        pdg::dependence *d;

        dep = isl_map_coalesce(dep);
        dep = simplify_controls(dep, info, &has_controls);
        if (isl_map_is_empty(dep)) {
                isl_map_free(dep);
                return isl_stat_ok;
        }

        info->used.insert(write_na);

        d = new pdg::dependence;
        d->array = info->a;
        d->type = info->dtype;
        d->from = write_na->node;
        d->to = info->read_na_pair->node;
        d->from_access = write_na->access;
        d->to_access = info->read_na_pair->access;

        if (d->from_access->extension || d->to_access->extension)
                d->extended_relation =
                        new pdg::IslMap(remove_all_controls(isl_map_copy(dep)));
        if (d->from_access->extension)
                dep = isl_map_apply_domain(dep,
                            isl_map_reverse(
                                d->from_access->extension->get_isl_map(ctx)));
        if (d->to_access->extension)
                dep = isl_map_apply_range(dep,
                            isl_map_reverse(
                                d->to_access->extension->get_isl_map(ctx)));
        if (has_controls)
                d->controlled_relation = new pdg::IslMap(isl_map_copy(dep));
        d->relation = new pdg::IslMap(remove_all_controls(isl_map_copy(dep)));
        info->pdg->dependences.push_back(d);
        isl_map_free(dep);
        return isl_stat_ok;
}

/* This structure represents a set of filter index expressions
 * along with bounds on the correponding filter values.
 * The number of output dimensions in "value" is the same as
 * the number of elements in the "index" vector.
 */
struct da_filter {
        std::vector<isl_union_map *> index;
        isl_map *value;
};

/* Construct a da_filter object representing the filters in
 * na->node and na->access.
 */
static struct da_filter *extract_filter(na_pair *na)
{
        da_filter *filter = new da_filter;
        isl_set *domain;
        pdg::node *node = na->node;
        pdg::access *access = na->access;

        domain = node->source->get_isl_set();
        if (isl_set_is_wrapping(domain))
                filter->value = isl_set_unwrap(domain);
        else
                filter->value = isl_map_from_domain(domain);

        for (int i = 0; i < node->filters.size(); ++i)
                filter->index.push_back(extract_access_map(node->filters[i]));

        if (access->nested.size() == 0)
                return filter;

        domain = isl_map_domain(access->map->get_isl_map());
        filter->value = isl_map_flat_range_product(filter->value,
                                                    isl_set_unwrap(domain));

        for (int i = 0; i < access->nested.size(); ++i)
                filter->index.push_back(extract_access_map(access->nested[i]));

        return filter;
}

static struct da_filter *da_filter_free(struct da_filter *filter)
{
        if (!filter)
                return NULL;
        for (int i = 0; i < filter->index.size(); ++i)
                isl_union_map_free(filter->index[i]);
        isl_map_free(filter->value);
        delete filter;
        return NULL;
}

static void da_filter_dump(struct da_filter *filter)
{
        isl_printer *p;

        if (!filter)
                return;

        p = isl_printer_to_file(isl_map_get_ctx(filter->value), stderr);
        p = isl_printer_start_line(p);
        p = isl_printer_print_str(p, "value(");
        for (int i = 0; i < filter->index.size(); ++i) {
                if (i)
                        p = isl_printer_print_str(p, ", ");
                p = isl_printer_print_union_map(p, filter->index[i]);
        }
        p = isl_printer_print_str(p, ") in ");
        p = isl_printer_print_map(p, filter->value);
        p = isl_printer_end_line(p);

        isl_printer_free(p);
}

/* Look for a filter index expression in "filter" that is identical
 * to "index".  Return the index of this index expression if it is
 * found and the number of elements in filter->index otherwise.
 */
static int da_filter_find_exact_match(struct da_filter *filter,
        __isl_keep isl_union_map *index)
{
        for (int i = 0; i < filter->index.size(); ++i) {
                int equal;

                equal = isl_union_map_is_equal(filter->index[i], index);
                if (equal < 0)
                        return -1;
                if (equal)
                        return i;
        }

        return filter->index.size();
}

/* Add the index expression "index" to "filter" with an unconstrained
 * filter value.  To ease debugging we set the name of the new filter
 * value dimension to that of the array being accessed by "index".
 * Although the range of "index" is allowed to live in more than
 * one space, we assume that they are all wrapped maps to the same
 * array space.
 */
static struct da_filter *da_filter_add(struct da_filter *filter,
        __isl_take isl_union_map *index)
{
        isl_id *id;
        isl_union_map *univ;
        isl_set *set;

        if (!filter || !index)
                goto error;

        filter->value = isl_map_add_dims(filter->value, isl_dim_out, 1);
        univ = isl_union_map_universe(isl_union_map_copy(index));
        univ = isl_union_set_unwrap(isl_union_map_range(univ));
        set = isl_set_from_union_set(isl_union_map_range(univ));
        id = isl_set_get_tuple_id(set);
        isl_set_free(set);
        filter->value = isl_map_set_dim_id(filter->value, isl_dim_out,
                                            filter->index.size(), id);
        if (!filter->value)
                goto error;
        filter->index.push_back(index);

        return filter;
error:
        isl_union_map_free(index);
        da_filter_free(filter);
        return NULL;
}

/* Intersect the set of possible filter values in "filter" with "value".
 */
static struct da_filter *da_filter_restrict(struct da_filter *filter,
        __isl_take isl_map *value)
{
        if (!filter || !value)
                goto error;

        filter->value = isl_map_intersect(filter->value, value);
        if (!filter->value)
                return da_filter_free(filter);

        return filter;
error:
        isl_map_free(value);
        da_filter_free(filter);
        return NULL;
}

/* Does "map" represent a total filter on "domain", i.e., one that is defined
 * on every element of "domain"?
 *
 * Although the range of "map" may live in different spaces, we assume
 * that the domain of "map" lives in a single space.
 */
static int is_total_filter(__isl_keep isl_union_map *map,
        __isl_keep isl_set *domain)
{
        isl_union_set *map_domain;
        isl_set *set;
        int total;

        map_domain = isl_union_map_domain(isl_union_map_copy(map));
        set = isl_set_from_union_set(map_domain);
        total = isl_set_is_subset(domain, set);
        isl_set_free(set);

        return total;
}

/* Does "node" have only total filters on "domain"?
 */
static int all_total_filters(pdg::node *node, __isl_take isl_set *domain)
{
        int total = 1;

        for (int i = 0; i < node->filters.size(); ++i) {
                isl_union_map *map = extract_access_map(node->filters[i]);

                total = is_total_filter(map, domain);
                isl_union_map_free(map);

                if (!total)
                        break;
        }

        isl_set_free(domain);
        return total;
}

/* Does "node" have only total filters on the domain of "map"?
 */
static bool filters_total_on_domain(pdg::node *node, __isl_keep isl_map *map)
{
        return all_total_filters(node, isl_map_domain(isl_map_copy(map)));
}

/* Does "node" have only total filters on the range of "map"?
 */
static bool filters_total_on_range(pdg::node *node, __isl_keep isl_map *map)
{
        return all_total_filters(node, isl_map_range(isl_map_copy(map)));
}

/* Are the filters of "source_node" total on the range of "source_map"
 * and those of "sink_node" (if any) total on the domain of "soruce_map"?
 */
static bool total_filters(pdg::node *source_node, pdg::node *sink_node,
        __isl_keep isl_map *source_map)
{
        int total;

        total = filters_total_on_range(source_node, source_map);
        if (!total)
                return false;

        if (!isl_set_is_wrapping(sink_node->source->set))
                return true;

        total = filters_total_on_domain(sink_node, source_map);
        if (!total)
                return false;

        return true;
}

/* Does "node" have any single valued filters?
 */
static int any_single_valued_filter(pdg::node *node)
{
        for (int i = 0; i < node->filters.size(); ++i) {
                isl_union_map *map = extract_access_map(node->filters[i]);
                int sv;

                sv = isl_union_map_is_single_valued(map);
                isl_union_map_free(map);

                if (sv)
                        return sv;
        }

        return 0;
}

/* Construct a mapping from the source iterations to all source filter values
 * that allow some corresponding sink iteration(s) (according to "source_map")
 * to be executed.  In other words, if any sink iteration is executed,
 * then we know that the filters of the corresponding source iterations
 * satisfy the returned relation.
 *
 * The "filter" input represents what is known about the filter
 * values at the sink.
 *
 * The "source" domain of the source is a wrapped map,
 * mapping iteration vectors to filter values.
 * We first construct a relation between the sink filter values (available
 * in "filter") and the source filter values.  The purpose of this relation
 * is to find as much information about the source filter values
 * as possible.  We can start with the value bounds on the arrays
 * accessed in the filters as they always hold.
 * Next, we loop over the source filters and check whether there
 * is any sink filter that covers the source filter.
 * In particular, for each source filter, we construct a map
 * from the source iteration domain to a wrapped access relation,
 * representing the write access relation that corresponds to
 * the filter read access.  Note that if we were unable to determine this
 * write access, then the mapping returned by extract_access_map
 * maps to the original read access, which will not match with
 * any filter access relations of the sink.
 * We combine the constructed map with the proto-dependence
 * (source_map) to obtain a mapping from sink iterations to
 * access relations such that there is some source iteration that
 * may be the source of the given sink iteration (based on source_map)
 * and that the filter value at this source iteration was written by
 * that access.  If the result is a subset of the mapping from the
 * sink iterations to the corresponding write access relation for some filter,
 * then we know that any constraint on this filter value also applies
 * to the source filter value.  We therefore introduce an equality
 * in our mapping from sink filter values to source filter values.
 *
 * When we apply the mapping from sink filter values to source filter values
 * to the mapping from source iterations to sink filter values, we
 * obtain a mapping from sink iterations to source filter values
 * that represents what we know about the source filter values.
 * That is, for each sink iteration in the domain of this map, if this
 * sink iteration is executed, then the actual source filter values
 * are an element of the image of the sink iteration.
 * In other words, the sink iteration is only executed if the source
 * filter values are an element of the image.
 * Mapping this relation back to the source through the proto-dependence
 * source_map, we obtain a relation from source iterations to all source
 * filter values for which any sink iteration is executed.
 * In particular, for values of the filters outside this relation,
 * no corresponding (according to source_map) sink is executed.
 */
static __isl_give isl_map *known_filter_values(struct da_filter *filter,
        na_pair *source_na, __isl_keep isl_map *source_map)
{
        isl_space *space, *space2;
        isl_map *source_value, *sink_value;
        isl_map *filter_map;
        isl_set *bounds;

        sink_value = isl_map_copy(filter->value);

        space = isl_set_get_space(source_na->node->source->set);
        space = isl_space_range(isl_space_unwrap(space));
        space2 = isl_space_range(isl_map_get_space(sink_value));
        space = isl_space_align_params(space, isl_space_copy(space2));
        space2 = isl_space_align_params(space2, isl_space_copy(space));
        space = isl_space_map_from_domain_and_range(space2, space);
        filter_map = isl_map_universe(space);

        bounds = compute_filter_bounds(source_na->node);
        filter_map = isl_map_intersect_range(filter_map, bounds);

        for (int i = 0; i < source_na->node->filters.size(); ++i) {
                isl_union_map *map_i;
                map_i = extract_access_map(source_na->node->filters[i]);
                map_i = isl_union_map_apply_range(
                    isl_union_map_from_map(isl_map_copy(source_map)), map_i);
                for (int j = 0; j < filter->index.size(); ++j) {
                        if (isl_union_map_is_subset(map_i, filter->index[j]))
                                filter_map = isl_map_equate(filter_map,
                                        isl_dim_in, j, isl_dim_out, i);
                }
                isl_union_map_free(map_i);
        }

        source_value = isl_map_apply_range(sink_value, filter_map);
        source_value = isl_map_apply_domain(source_value,
                                                isl_map_copy(source_map));
        source_value = isl_map_coalesce(source_value);

        return source_value;
}

/* Add a new output dimension to "map" with constraints that are the
 * same as those on output dimension "pos".
 *
 * Given map { [i] -> [j] }, we first and an extra dimension,
 *
 *      { [i] -> [j,*] }
 *
 * extract out j_pos,
 *
 *      { [[i] -> [j_0,...,j_{pos-1},*,j_{pos+1},...,*]] -> [j_pos] }
 *
 * create a copy
 *
 *      { [[i] -> [j_0,...,j_{pos-1},*,j_{pos+1},...,*]] -> [j_pos,j_pos'] }
 *
 * and then move the dimensions back
 *
 *      { [i] -> [j,j_pos'] }
 */
static __isl_give isl_map *copy_dim(__isl_take isl_map *map, int pos)
{
        int pos_new;

        pos_new = isl_map_dim(map, isl_dim_out);
        pos += isl_map_dim(map, isl_dim_in);
        pos_new += isl_map_dim(map, isl_dim_in);
        map = isl_map_add_dims(map, isl_dim_out, 1);
        map = isl_map_from_domain(isl_map_wrap(map));
        map = isl_map_add_dims(map, isl_dim_out, 1);
        map = isl_map_equate(map, isl_dim_in, pos, isl_dim_out, 0);
        map = isl_map_eliminate(map, isl_dim_in, pos, 1);
        map = isl_map_range_product(map, isl_map_copy(map));
        map = isl_map_equate(map, isl_dim_in, pos, isl_dim_out, 0);
        map = isl_map_equate(map, isl_dim_in, pos_new, isl_dim_out, 1);
        map = isl_set_unwrap(isl_map_domain(map));

        return map;
}

/* Given constraints on the filter values "filter" at the sink iterations
 * "sink", derive additional constraints from the filter values of source
 * node "source_na".  In particular, consider the iterations of "source_na"
 * that have _not_ been executed based on the constraints of the corresponding
 * last iteration parameters in "sink" and what this implies about the
 * filter values at those iterations.
 *
 * Essentially, we consider all pairs of sink iterations and filter
 * elements, together with the corresponding non-executed source iterations
 * and the possible values of those filters.  We universally quantify
 * the non-executed source iterations so that we obtain the intersection
 * of the constraints on the filter values over all those source iterations
 * and then existentially quantify the filter elements to obtain constraints
 * that are valid for all filter elements.
 *
 * In more details, the computation is performed as follows.
 *
 * Compute a mapping from potential last iterations of the other source
 * to sink iterations, taking into account the contraints on
 * the last executed iterations encoded in the parameters of "sink",
 * but projecting out all parameters encoding last iterations from the result.
 * Include all earlier iterations of the other source, resulting in
 * a mapping with a domain that includes all potential iterations of the
 * other source.
 *
 * Subtract these iterations from all possible iterations of the other
 * source for a given sink iteration, resulting in a mapping from
 * potential source iterations that are definitely not executed
 * to the corresponding sink iteration.
 *
 * If this map is empty, then this means we can't find any iterations
 * of the other source that are certainly not executed and then we
 * can't derive any further information.
 * Similarly, if the filters of source_na->node are not total
 * on the set of non-executed iterations, then we cannot draw any conclusions.
 * (Note that we already tested that the filters on sink->node are total
 * on the domain of "source_map".  By intersecting the set of corresponding
 * sink iterations with this domain, we ensure that this property also holds
 * on those sink iterations.)
 * Otherwise, keep track of those sink iterations without any corresponding
 * non-executed other source iterations.  We will lose these sink iterations
 * in subsequent computations, so we need to add them back in at the end.
 *
 * Compute bounds on the filter values at the non executed iterations
 * based on what we know about filters at the sink and the fact that
 * the iterations are not executed, meaning that the filter values
 * do not satisfy the constraints that allow the iteration to be executed.
 * The result is a mapping T -> V.
 *
 * Note that we only know that there is some accessed filter element
 * that does not satisfy the constraints that allow the iteration to be
 * executed.  We therefore project out those dimensions that correspond
 * to filters with an access relation that is not single-valued, i.e.,
 * one that may access more than one element for some iterations.
 * If there are no single-valued filters, then we can skip the rest of
 * the computation.
 *
 * Construct a mapping [K -> F] -> T, with K the sink iterations,
 * T the corresponding non-executed iterations of the other source and
 * F the filters accessed at those iterations.
 *
 * Combine the above two mappings into a mapping [K -> F] -> V
 * such that the set of possible filter values (V) is the intersection
 * over all iterations of the other source that access the filter,
 * and such that there is at least one such iteration.
 * In other words, ensure that the range of [K -> F] -> T
 * is a non-empty subset of the range of V -> T.
 * We require a non-empty subset to ensure that the domain of
 * [K -> F] -> V is equal to the result of composing K -> T with T -> F.
 * Projecting out F from [K -> F] -> V, we obtain a map K -> V that is
 * the union of all possible values of the filters K -> F, i.e.,
 * the constraints that the values of K -> F satisfy.
 * If we didn't impose non-emptiness above, then the domain of [K -> F] -> V
 * would also include pairs that we are not interested in, related to
 * arbitrary values.  Projecting out F would then also lead to arbitrary
 * values.
 *
 * Compose the mapping K -> T with the index expressions, pulling them
 * back to the sink.
 * For each of these pulled back index expressions, we check if it
 * is equal to one of the sink filter index expressions.  If not, we
 * add it to the sink filter index expressions.
 * In both cases, we keep track of the fact that this sink filter
 * should have a value that satisfies the constraints in K -> V.
 * We further check if there is any sink filter index expression
 * that is a (strict) subset of the pulled back index expression.
 * The value of any such sink filter should also satisfy those
 * constraints, so we duplicate the filter value in K -> V.
 *
 * Finally, we intersect the possible filter values with the constraints
 * obtained above on the affected sink iterations and a universe range
 * on the unaffected sink iterations.
 */
static struct da_filter *include_other_source(struct da_filter *filter,
        __isl_keep isl_map *source_map, __isl_keep isl_set *sink,
        na_pair *source_na, add_dep_info *info)
{
        isl_space *space, *space2;
        isl_set *source_domain;
        isl_map *mem;
        isl_map *may_run, *not_run;
        isl_map *source_value;
        isl_union_map *usource, *umap;
        std::vector<isl_union_map *> index;
        isl_union_map *index_product;
        isl_map *map;
        isl_map *filter_map;
        isl_set *unaffected_sink;
        isl_map *unaffected_value;

        if (source_na->node->filters.size() == 0)
                return filter;

        space = isl_set_get_space(source_na->node->source->set);
        space = isl_space_domain(isl_space_unwrap(space));
        source_domain = isl_set_universe(space);
        source_domain = add_parametrization(source_domain, source_na, info);
        may_run = isl_map_from_domain_and_range(source_domain,
                                                isl_set_copy(sink));
        may_run = share(may_run, info->n_shared);

        may_run = remove_all_controls(may_run);
        mem = info->get_mem_dep(source_na);
        mem = share(mem, info->n_shared);
        mem = isl_map_intersect_range(mem,
                isl_map_domain(remove_all_controls(isl_map_copy(source_map))));
        space = isl_space_domain(isl_map_get_space(may_run));
        may_run = isl_map_apply_domain(may_run, isl_map_lex_ge(space));
        not_run = isl_map_subtract(mem, may_run);

        if (isl_map_is_empty(not_run) ||
            !any_single_valued_filter(source_na->node) ||
            !filters_total_on_domain(source_na->node, not_run)) {
                isl_map_free(not_run);
                return filter;
        }

        not_run = isl_map_reverse(not_run);
        source_value = known_filter_values(filter, source_na, not_run);
        source_value = isl_map_subtract(source_value,
                        isl_set_unwrap(source_na->node->source->get_isl_set()));
        unaffected_sink = isl_set_copy(sink);
        unaffected_sink = remove_all_controls(unaffected_sink);
        unaffected_sink = isl_set_subtract(unaffected_sink,
                                        isl_map_domain(isl_map_copy(not_run)));

        for (int i = 0; i < source_na->node->filters.size(); ++i) {
                isl_union_map *map;
                int sv;

                map = extract_access_map(source_na->node->filters[i]);
                sv = isl_union_map_is_single_valued(map);
                if (sv) {
                        index.push_back(map);
                } else {
                        isl_union_map_free(map);
                        source_value = isl_map_project_out(source_value,
                                                isl_dim_out, index.size(), 1);
                }
        }

        index_product = isl_union_map_copy(index[0]);
        for (int i = 1; i < index.size(); ++i)
                index_product = isl_union_map_range_product(index_product,
                                                isl_union_map_copy(index[1]));
        index_product = isl_union_map_reverse(index_product);
        index_product = isl_union_map_domain_product(
                isl_union_map_from_map(isl_map_copy(not_run)), index_product);

        usource = isl_union_map_from_map(source_value);
        usource = join_non_empty_subset(index_product, usource);
        umap = isl_union_map_universe(isl_union_map_copy(usource));
        umap = isl_union_set_unwrap(isl_union_map_domain(umap));
        umap = isl_union_map_domain_map(umap);
        usource = isl_union_map_apply_domain(usource, umap);
        source_value = isl_map_from_union_map(usource);

        for (int i = 0; i < index.size(); ++i)
                index[i] = isl_union_map_apply_range(
                    isl_union_map_from_map(isl_map_copy(not_run)), index[i]);

        space = isl_space_range(isl_map_get_space(source_value));
        space2 = isl_space_range(isl_map_get_space(filter->value));
        space = isl_space_align_params(space, isl_space_copy(space2));
        space2 = isl_space_align_params(space2, isl_space_copy(space));
        space = isl_space_map_from_domain_and_range(space, space2);
        filter_map = isl_map_universe(space);

        for (int i = 0; i < index.size(); ++i) {
                int exact = da_filter_find_exact_match(filter, index[i]);
                assert(exact >= 0);
                if (exact == filter->index.size()) {
                        filter = da_filter_add(filter,
                                                isl_union_map_copy(index[i]));
                        filter_map = isl_map_add_dims(filter_map,
                                                        isl_dim_out, 1);
                }
                filter_map = isl_map_equate(filter_map, isl_dim_in, i,
                                                        isl_dim_out, exact);
                for (int j = 0; j < filter->index.size(); ++j) {
                        int pos;

                        if (j == exact)
                                continue;
                        if (!isl_union_map_is_subset(filter->index[j],
                                                        index[i]))
                                continue;
                        pos = isl_map_dim(source_value, isl_dim_out);
                        source_value = copy_dim(source_value, i);
                        filter_map = isl_map_add_dims(filter_map,
                                                        isl_dim_in, 1);
                        filter_map = isl_map_equate(filter_map, isl_dim_in, pos,
                                                                isl_dim_out, j);
                }
        }

        source_value = isl_map_apply_range(source_value, filter_map);
        source_value = isl_map_coalesce(source_value);
        unaffected_value = isl_map_from_domain(unaffected_sink);
        unaffected_value = isl_map_add_dims(unaffected_value, isl_dim_out,
                                                filter->index.size());
        source_value = isl_map_union(source_value, unaffected_value);
        filter = da_filter_restrict(filter, source_value);

        for (int i = 0; i < index.size(); ++i)
                isl_union_map_free(index[i]);
        isl_map_free(not_run);

        return filter;
}

/* Given constraints on the filter values "filter" at the sink iterations
 * "sink", derive additional constraints from the filter values of those
 * source nodes for which "sink" contains a reference to its last iteration,
 * for use in determining whether parametrization is needed on "source_map".
 * In particular, we try and derive extra information from the fact that
 * some iterations of those source nodes have _not_ been executed.
 */
static struct da_filter *include_other_sources(struct da_filter *filter,
        __isl_keep isl_map *source_map, __isl_keep isl_set *sink,
        add_dep_info *info)
{
        isl_space *space;
        int n_param;

        space = isl_set_get_space(sink);
        n_param = isl_space_dim(space, isl_dim_param);
        for (int i = 0; i < n_param; ++i) {
                isl_id *id;
                na_pair *na;

                if (!is_shared(space, i))
                        continue;

                id = isl_space_get_dim_id(space, isl_dim_param, i);
                na = (na_pair *) isl_id_get_user(id);
                isl_id_free(id);

                filter = include_other_source(filter, source_map, sink,
                                                na, info);
        }
        isl_space_free(space);

        return filter;
}

#define RESTRICT_ERROR  -1
#define RESTRICT_NO     0
#define RESTRICT_EMPTY  1
#define RESTRICT_INPUT  2
#define RESTRICT_OUTPUT 3

/* Given a map from sinks to potential sources (source_map)
 * and the set of sink iterations (sink),
 * check if any parametrization is needed on the sources.
 * That is, check whether the possible filter values at the sink
 * imply that the filter values at the source are always valid.
 * If so, the source is executed whenever the sink is executed
 * and no parametrization is required.
 *
 * Return RESTRICT_NO if no parametrization is required.
 * Return RESTRICT_INPUT if parametrization is required on the input
 * of the computation of the last iteration.
 * Return RESTRICT_OUTPUT if parametrization is required on the output
 * of the computation of the last iteration.  This means that we know
 * that the source will be executed, but we want to introduce a parameter
 * to represent the last iteration anyway, because the knowledge depends
 * on the parameters representing last iterations of other nodes.
 * Return RESTRICT_EMPTY if the potential sources cannot possibly
 * be executed, assuming that the sink is executed.
 *
 * If there are no filters on the source, then obviously the source
 * is always executed.
 *
 * If the filters of the sink and the source are not all total
 * on domain and range of "source_map", then we cannot draw any conclusion.
 * In principle, we could split up "source_map" according to whether
 * the filters would be total on the domain and range.
 *
 * We first construct a mapping from source iterations to source filter
 * values that allow some corresponding sink iteration(s) (according to
 * "source_map") to be executed.
 * If this relation is a subset of the actual mapping from iteration
 * vectors to filter values at the source, then we know that a corresponding
 * sink is only executed when the source is executed and no parametrization
 * is required.  However, we postpone the decision until we have considered
 * the other potential sources below.
 * If, on the other hand, the constructed relation is disjoint
 * from the source filter relation, then the sources cannot have
 * executed if the sink is executed.  If so, we return
 * RESTRICT_EMPTY immediately.
 *
 * Otherwise, we check if we can find out more information by considering
 * information derived from knowledge about the last iterations of other
 * nodes.  If, by considering this extract information, we can find
 * that the potential source is never executed (given that the sink
 * is executed), then we return RESTRICT_EMPTY.
 * Otherwise, if we had already determined that the relation based
 * on only the sink is a subset of the filter values, then we return
 * RESTRICT_NO.  If we can only draw this conclusion when taking into
 * account the other potential sources, then we return RESTRICT_OUTPUT.
 * Otherwise, we return RESTRICT_INPUT.
 */
static int need_parametrization(__isl_keep isl_map *source_map,
        __isl_keep isl_set *sink, na_pair *source_na, add_dep_info *info)
{
        bool filtered_source;
        isl_space *space;
        isl_map *source_value, *sink_value;
        da_filter *filter;
        na_pair *sink_na = info->read_na_pair;
        int res;

        if (isl_map_plain_is_empty(source_map) ||
            isl_set_plain_is_empty(sink))
                return RESTRICT_NO;

        filtered_source = isl_set_is_wrapping(source_na->node->source->set);

        if (!filtered_source)
                return RESTRICT_NO;

        if (!total_filters(source_na->node, sink_na->node, source_map))
                return RESTRICT_INPUT;

        filter = extract_filter(sink_na);
        res = -1;
        if (filter->index.size() != 0) {
                sink_value = known_filter_values(filter, source_na, source_map);
                source_value = isl_set_unwrap(
                                        source_na->node->source->get_isl_set());

                if (isl_map_is_disjoint(sink_value, source_value))
                        res = RESTRICT_EMPTY;
                else if (isl_map_is_subset(sink_value, source_value))
                        res = RESTRICT_NO;
                isl_map_free(source_value);
                isl_map_free(sink_value);
        }

        if (res == RESTRICT_EMPTY) {
                da_filter_free(filter);
                return res;
        }

        filter = include_other_sources(filter, source_map, sink, info);

        if (!filter)
                return RESTRICT_ERROR;
        if (filter->index.size() == 0) {
                da_filter_free(filter);
                return RESTRICT_INPUT;
        }

        sink_value = known_filter_values(filter, source_na, source_map);
        source_value = isl_set_unwrap(source_na->node->source->get_isl_set());

        if (isl_map_is_disjoint(sink_value, source_value))
                res = RESTRICT_EMPTY;
        else if (res == RESTRICT_NO)
                ;
        else if (isl_map_is_subset(sink_value, source_value))
                res = RESTRICT_OUTPUT;
        else
                res = RESTRICT_INPUT;

        isl_map_free(source_value);
        isl_map_free(sink_value);
        da_filter_free(filter);

        return res;
}

extern "C" {
        static __isl_give isl_restriction *do_restrict(
                __isl_keep isl_map *source_map, __isl_keep isl_set *sink,
                void *source_user, void *user);
}

/* Add parameters corresponding to the last iteration of "na" to "set"
 * (assuming they don't already appear in "set")
 * and add constraints to them to express that there either is no
 * last iteration with info->n_shared shared loops
 * (__last_<na>_shared < info->n_shared) or that there is a last
 * iteration with at least info->n_shared shared loops
 * (__last_<na>_shared >= info->n_shared) and that the iteration is a possible
 * source of the current sink (based on the memory dependence between
 * the source and the sink).
 */
static __isl_give isl_set *set_parameter_bounds(__isl_take isl_set *set,
        na_pair *na, add_dep_info *info)
{
        isl_map *mem;
        isl_id *id;
        isl_set *valid;
        isl_set *invalid;
        isl_set *domain;
        int shared_pos;

        id = create_shared_id(isl_set_get_ctx(set), na);
        shared_pos = find_or_add_param(&set, id);

        valid = isl_set_copy(set);
        invalid = set;

        valid = isl_set_lower_bound_si(valid, isl_dim_param, shared_pos,
                                        info->n_shared);

        mem = info->get_mem_dep(na);
        domain = isl_set_universe(isl_space_domain(isl_map_get_space(mem)));
        domain = add_parametrization(domain, na, info);
        mem = isl_map_intersect_domain(mem, domain);
        valid = isl_set_intersect(valid, isl_map_range(mem));

        invalid = isl_set_upper_bound_si(invalid, isl_dim_param, shared_pos,
                                        info->n_shared - 1);

        return isl_set_union(valid, invalid);
}

/* Check if there are any iterations of "source_na" in "source_map"
 * that are definitely executed, based solely on the possible filter values.
 * If so, add constraints to "sink" to indicate that the last execution
 * cannot be earlier than those definitely executed iterations.
 *
 * We first compute the set of source iterations that are definitely
 * executed because there are no filter values that would prohibit
 * their execution.  If there are no such source iterations then we are done.
 *
 * Then we construct a map from sink iterations to associated (through
 * "source_map") definitely executed source iterations.
 *
 * For those sink iterations that have a corresponding definitely
 * executed source iteration, we add constraints that express that
 * this last definitely executed source iteration is lexicographically
 * smaller than or equal to the last executed source iteration
 * (and that there definitely is a last executed source iteration).
 */
static __isl_give isl_set *mark_definite_source(__isl_take isl_set *sink,
        add_dep_info *info, na_pair *source_na, __isl_keep isl_map *source_map)
{
        isl_space *space;
        isl_set *dom;
        isl_map *invalid;
        isl_map *definite_source_map;
        isl_set *with_source;
        isl_map *map_after;
        int depth;

        dom = source_na->node->source->get_isl_set();
        dom = isl_map_domain(isl_set_unwrap(dom));
        invalid = compute_filter_values(source_na->node, true);
        dom = isl_set_subtract(dom, isl_map_domain(invalid));

        if (isl_set_is_empty(dom)) {
                isl_set_free(dom);
                return sink;
        }

        space = isl_set_get_space(dom);

        definite_source_map = isl_map_copy(source_map);
        definite_source_map = isl_map_intersect_range(source_map, dom);

        dom = isl_map_domain(isl_map_copy(definite_source_map));

        with_source = isl_set_copy(sink);
        with_source = isl_set_intersect(with_source, isl_set_copy(dom));
        sink = isl_set_subtract(sink, dom);

        dom = isl_set_universe(isl_space_copy(space));
        dom = add_parametrization(dom, source_na, info);

        depth = source_na->node->get_filter_depth();

        space = isl_space_map_from_set(space);
        map_after = isl_map_lex_ge_first(space, depth);
        dom = isl_set_apply(dom, map_after);

        definite_source_map = isl_map_intersect_range(definite_source_map, dom);

        dom = isl_map_domain(definite_source_map);
        with_source = isl_set_intersect(with_source, dom);

        sink = isl_set_union(sink, with_source);

        return sink;
}

/* Remove inconsistencies from the set of sink iterations "sink"
 * based on the current potential source "source_na" and other
 * potential sources referenced by "sink".
 *
 * We first identify those iterations of "source_na" that are
 * definitely executed based solely on the possible filter values.
 *
 * If the sink has filters, then we remove inconsistencies based
 * on the sink and the current potential source.
 *
 * Finally, we go through the references to other potential sources
 * in "sink" and remove inconsistencies based on this other potential
 * source and the current potential source.
 */
static __isl_give isl_set *remove_inconsistencies(__isl_take isl_set *sink,
        add_dep_info *info, na_pair *source_na, __isl_keep isl_map *source_map)
{
        isl_space *space;
        isl_id *source_id;
        int n_param;

        source_id = create_shared_id(isl_map_get_ctx(source_map), source_na);

        sink = mark_definite_source(sink, info, source_na, source_map);

        if (isl_set_is_wrapping(info->read_na_pair->node->source->set))
                sink = remove_inconsistencies(sink, info, source_id);

        space = isl_set_get_space(sink);
        n_param = isl_space_dim(space, isl_dim_param);
        for (int i = 0; i < n_param; ++i) {
                isl_id *other_id;

                if (!is_shared(space, i))
                        continue;
                other_id = isl_space_get_dim_id(space, isl_dim_param, i);

                if (other_id != source_id) {
                        na_pair *other_na;

                        other_na = (na_pair *) isl_id_get_user(other_id);
                        sink = set_parameter_bounds(sink, other_na, info);
                        sink = remove_inconsistencies(sink, info, other_id,
                                                        source_id);
                }

                isl_id_free(other_id);
        }
        isl_space_free(space);

        isl_id_free(source_id);
        return sink;
}

/* Compute a restriction for the given sink.
 * That is, add constraints to the parameters expressing
 * that the source is either not executed with info->n_shared shared
 * iterators (*_shared < info->n_shared)
 * or it is executed (*_shared >= info->n_shared) and then the last iteration
 * satisfies the corresponding memory based dependence.
 *
 * Only do this for that part of "sink" that has any corresponding
 * sources in "source_map".  The remaining part of "sink" is not affected.
 *
 * Note that the "sink" set may have undergone a refinement based
 * on the _shared parameters and we want this refinement to also
 * be present in the sink restriction.  We therefore need
 * to intersect the affected part with "sink".
 */
static __isl_give isl_set *compute_sink_restriction(
        __isl_keep isl_map *source_map, __isl_keep isl_set *sink,
        na_pair *source_na, add_dep_info *info)
{
        isl_set *with_source, *without_source;

        with_source = isl_map_domain(isl_map_copy(source_map));
        without_source = isl_set_subtract(isl_set_copy(sink),
                                            isl_set_copy(with_source));
        with_source = isl_set_intersect(with_source, isl_set_copy(sink));

        with_source = set_parameter_bounds(with_source, source_na, info);
        with_source = remove_inconsistencies(with_source, info, source_na,
                                                source_map);

        return isl_set_union(with_source, without_source);
}

/* How many loops do "node1" and "node2" share?
 */
static int max_shared(pdg::node *node1, pdg::node *node2)
{
        int shared = 0;
        int size = node1->prefix.size();

        if (node2->prefix.size() < size)
                size = node2->prefix.size();

        for (int i = 0; i < size && node1->prefix[i] == node2->prefix[i]; ++i)
                if (node1->prefix[i] == -1)
                        ++shared;

        return shared;
}

/* Determine the number of shared loop iterators between sink
 * and source domains in "sink2source".  That is, find out how
 * many of the initial input and output dimensions are equal
 * to each other.  Return the result.
 *
 * The first time this function is called for a given sink access
 * (info->n_shared is set to -1 in add_dep_info::set_read_na),
 * we check for equal dimensions up to the shared nesting depth.
 * Later call check dimensions up to the result of the previous call.
 */
static int extract_shared_levels(__isl_keep isl_map *sink2source,
        na_pair *source_na, add_dep_info *info)
{
        int i;
        int max;
        isl_space *space;

        max = info->n_shared;
        if (max < 0)
                max = max_shared(source_na->node, info->read_na_pair->node);

        if (max == 0)
                return info->n_shared = 0;

        space = isl_map_get_space(sink2source);
        for (i = 0; i < max; ++i) {
                isl_map *test;
                int subset;

                test = isl_map_universe(isl_space_copy(space));
                test = isl_map_equate(test, isl_dim_in, i, isl_dim_out, i);
                subset = isl_map_is_subset(sink2source, test);
                isl_map_free(test);

                if (!subset)
                        break;
        }
        isl_space_free(space);

        return i;
}

/* The last iteration referred to by the sink may have been added
 * at a different nesting level.  This means that __last_<na>_shared
 * is greater than or equal to a value greater than info->n_shared
 * and that therefore the iterators between info->n_shared and
 * __last_<na>_shared are not represented as they are implicitly
 * considered to be equal to the corresponding sink iterator.
 * For consistency, we need to explicitly add those iterators
 * and set them to be equal to the corresponding sink iterator.
 *
 * In particular, we create a set in the space of the sink of the form
 *
 *      { s : __last_<na>_shared < info->n_shared or
 *              (__last_<na>_<i> is a potential source iteration for s and
 *               (__last_<na>_shared >= max_shared or
 *                (__last_<na>_shared = max_shared - 1 and
 *                      the first max_shared - 1 iterators are equal and
 *                      iterator max_shared - 1 of the source is smaller) or
 *                ...
 *                (__last_<na>_shared = info->n_shared and
 *                      the first n_shared iterators are equal and
 *                      iterator n_shared of the source is smaller))) }
 *
 * That is, for those parts with __last_<na>_shared smaller than
 * max_n_shared[source_na], intersection with the set will introduce
 * __last_<na>_<i> parameters (assuming they don't have a known fixed value)
 * up until __last_<na>_shared and equate them to the corresponding iterators
 * of the sink.
 */
static __isl_give isl_set *internal_shared_refinement(
        __isl_keep isl_id *shared_id, add_dep_info *info)
{
        isl_map *mem;
        na_pair *source_na;
        isl_set *domain;
        isl_map *shared_map;
        int shared_pos;
        isl_set *valid;
        isl_set *invalid;

        source_na = (na_pair *) isl_id_get_user(shared_id);

        mem = info->get_mem_dep(source_na);
        domain = isl_set_universe(isl_space_domain(isl_map_get_space(mem)));
        domain = add_parametrization(domain, source_na, info);
        mem = isl_map_intersect_domain(mem, domain);

        shared_map = compute_shared_map(isl_map_get_space(mem), shared_id,
                info->n_shared, info->max_n_shared[source_na]);

        mem = isl_map_intersect(mem, shared_map);

        valid = isl_map_range(mem);
        invalid = isl_set_universe(isl_set_get_space(valid));
        shared_pos = isl_set_find_dim_by_id(invalid, isl_dim_param, shared_id);
        invalid = isl_set_upper_bound_si(invalid, isl_dim_param, shared_pos,
                                        info->n_shared - 1);

        return isl_set_union(valid, invalid);
}

/* The last iteration of some source referred to by the sink may have been
 * added at a different nesting level.  This means that __last_*_shared
 * is greater than or equal to a value greater than info->n_shared
 * and that therefore the iterators between info->n_shared and
 * __last_*_shared are not represented as they are implicitly
 * considered to be equal to the corresponding sink iterator.
 * For consistency, we need to explicitly add those iterators
 * and set them to be equal to the corresponding sink iterator.
 *
 * For each of the __last_*_shared parameters, explicitly add
 * the implicitly equal __last_*_i iterators by intersecting
 * the sink with the set computed by internal_shared_refinement.
 */
static __isl_give isl_set *refine_shared_internal(__isl_take isl_set *sink,
        add_dep_info *info)
{
        isl_space *space;
        isl_set *refinement;
        int n_param;

        space = isl_set_get_space(sink);
        refinement = isl_set_universe(isl_space_copy(space));

        n_param = isl_space_dim(space, isl_dim_param);
        for (int i = 0; i < n_param; ++i) {
                isl_set *ref_i;
                isl_id *id;

                if (!is_shared(space, i))
                        continue;

                id = isl_space_get_dim_id(space, isl_dim_param, i);
                ref_i = internal_shared_refinement(id, info);
                isl_id_free(id);
                refinement = isl_set_intersect(refinement, ref_i);
        }
        isl_space_free(space);

        sink = isl_set_intersect(sink, refinement);

        return sink;
}

/* Given a map from sinks to potential sources (sink2source),
 * check if any parametrization is needed.
 * Depending on the result, return either a universe restriction,
 * an empty restriction (if the sources cannot have executed),
 * a restriction that parametrizes the source and the sink
 * of the input of the computation of the last source
 * or a restriction that parametrizes the source of the output.
 *
 * sink_map maps the domain of sink2source to the sink iteration domain.
 * source_map maps the range of sink2source to the source iteration domain.
 *
 * Before we check if we need any parametrization, we update the number of
 * shared loop levels and add possibly missing __last_*_i iterators
 * (in refine_shared_internal).  If parametrization turns out to be required,
 * we also update the minimal number of shared loop levels for the given
 * source.
 */
static __isl_give isl_restriction *compute_restriction_core(
        __isl_keep isl_map *sink2source,
        __isl_take isl_map *sink_map, __isl_take isl_map *source_map,
        __isl_keep isl_set *sink, na_pair *source_na, add_dep_info *info)
{
        isl_space *space;
        isl_set *source_restr;
        isl_set *sink_restr;
        int need;

        sink2source = isl_map_copy(sink2source);
        sink = isl_set_copy(sink);

        sink2source = isl_map_apply_range(sink2source,
                                isl_map_copy(source_map));
        sink2source = isl_map_apply_domain(sink2source,
                                isl_map_copy(sink_map));
        sink = isl_set_apply(sink, isl_map_copy(sink_map));

        info->n_shared = extract_shared_levels(sink2source, source_na, info);
        sink = refine_shared_internal(sink, info);

        need = need_parametrization(sink2source, sink, source_na, info);
        if (need == RESTRICT_ERROR ||
            need == RESTRICT_NO || need == RESTRICT_EMPTY) {
                isl_map_free(source_map);
                isl_map_free(sink_map);
                isl_set_free(sink);
                if (need == RESTRICT_ERROR) {
                        isl_map_free(sink2source);
                        return NULL;
                } else if (need == RESTRICT_NO)
                        return isl_restriction_none(sink2source);
                else
                        return isl_restriction_empty(sink2source);
        }

        info->update_min_n_shared(source_na);

        space = isl_map_get_space(source_map);
        source_restr = isl_set_universe(isl_space_range(space));

        source_restr = add_parametrization(source_restr, source_na, info);
        source_restr = isl_set_apply(source_restr, isl_map_reverse(source_map));

        if (need == RESTRICT_OUTPUT) {
                isl_map_free(sink_map);
                isl_set_free(sink);
                isl_map_free(sink2source);
                return isl_restriction_output(source_restr);
        }

        sink_restr = compute_sink_restriction(sink2source, sink,
                                                source_na, info);
        sink_restr = isl_set_apply(sink_restr, isl_map_reverse(sink_map));

        isl_set_free(sink);
        isl_map_free(sink2source);

        return isl_restriction_input(source_restr, sink_restr);
}

/* Compute a restriction for the given map from sinks to potential sources
 * (sink2source).
 *
 * First check if the sink access has any filters.  If so, compose the original
 * sink_map with a mapping that projects out these access filters.
 * Handle the source access similarly.
 * Then call compute_restriction_core to perform the main computation.
 */
static __isl_give isl_restriction *compute_restriction(
        __isl_keep isl_map *sink2source,
        __isl_take isl_map *sink_map, __isl_take isl_map *source_map,
        __isl_keep isl_set *sink, na_pair *source_na, add_dep_info *info)
{
        na_pair *sink_na = info->read_na_pair;

        if (sink_na->access->nested.size() > 0) {
                isl_space *space;
                isl_map *map;

                space = isl_space_range(isl_map_get_space(sink_map));
                space = isl_space_unwrap(space);
                map = isl_map_domain_map(isl_map_universe(space));

                sink_map = isl_map_apply_range(sink_map, map);
        }

        if (source_na->access->nested.size() > 0) {
                isl_space *space;
                isl_map *map;

                space = isl_space_range(isl_map_get_space(source_map));
                space = isl_space_unwrap(space);
                map = isl_map_domain_map(isl_map_universe(space));

                source_map = isl_map_apply_range(source_map, map);
        }

        return compute_restriction_core(sink2source, sink_map, source_map,
                                        sink, source_na, info);
}

/* Compute a restriction for the given map from sinks to potential sources
 * (sink2source).  We simply call compute_restriction to compute the
 * restriction.  Since, unlike the case of do_restrict_domain_map bewloe,
 * we didn't encode the entire access relation in the domains of the input
 * to isl_access_info_compute_flow, we pass identity mappings on the source
 * and sink to compute_restriction.
 */
static __isl_give isl_restriction *do_restrict(__isl_keep isl_map *sink2source,
        __isl_keep isl_set *sink, void *source_user, void *user)
{
        na_pair *source_na = (na_pair *) source_user;
        add_dep_info *info = (struct add_dep_info *) user;
        isl_space *space;
        isl_map *source_map;
        isl_map *sink_map;

        space = isl_space_domain(isl_map_get_space(sink2source));
        sink_map = isl_map_identity(isl_space_map_from_set(space));
        space = isl_space_range(isl_map_get_space(sink2source));
        source_map = isl_map_identity(isl_space_map_from_set(space));

        return compute_restriction(sink2source, sink_map, source_map,
                                        sink, source_na, info);
}

/* Does the iteration domain of any of the writes involve any filters?
 */
static bool any_filters(vector<na_pair> &writers)
{
        for (int i = 0; i < writers.size(); ++i)
                if (isl_set_is_wrapping(writers[i].node->source->set))
                        return true;
        return false;
}

/* Does any of the dependences starting at "first" have
 * controlled dependence relation?
 */
static bool any_controlled_dependences(PDG *pdg, int first)
{
        for (int i = first; i < pdg->dependences.size(); ++i)
                if (pdg->dependences[i]->controlled_relation)
                        return true;

        return false;
}

/* Remove parameters from "map" that start with "prefix".
 */
static __isl_give isl_map *remove_source(__isl_take isl_map *map,
        const char *prefix)
{
        int n = isl_map_dim(map, isl_dim_param);
        size_t len = strlen(prefix);

        for (int i = n - 1; i >= 0; --i) {
                const char *name;

                name = isl_map_get_dim_name(map, isl_dim_param, i);
                if (strncmp(name, prefix, len))
                        continue;

                map = isl_map_project_out(map, isl_dim_param, i, 1);
        }

        map = isl_map_coalesce(map);

        return map;
}

/* Remove parameters from the dependences starting at "first"
 * that refer to any of the unused potential sources, i.e.,
 * those potential sources that are in writers, but not in info->used.
 *
 * Since the current sink does not depend on those unused potential sources,
 * the corresponding dependence relations cannot depend on them and
 * any reference to them can simply be projected out.
 */
static void remove_unused_sources(PDG *pdg, int first,
        vector<na_pair> &writers, add_dep_info *info)
{
        char name[60];
        std::vector<na_pair *> unused;

        for (int i = 0; i < writers.size(); ++i) {
                if (info->used.find(&writers[i]) == info->used.end())
                        unused.push_back(&writers[i]);
        }

        if (unused.size() == 0)
                return;
        if (!any_controlled_dependences(pdg, first))
                return;

        for (int i = 0; i < unused.size(); ++i) {
                na_pair *na = unused[i];

                snprintf(name, sizeof(name), "__last_%s_%d_",
                        na->node->name->s.c_str(), na->access->nr);

                for (int j = first; j < pdg->dependences.size(); ++j) {
                        pdg::dependence *dep = pdg->dependences[j];
                        isl_map *map;

                        if (!dep->controlled_relation)
                                continue;

                        map = dep->controlled_relation->map;
                        map = remove_source(map, name);
                        dep->controlled_relation->map = map;
                }
        }
}

/* Look for the unique write access that writes to the array accessed
 * by "a" and return an na_pair consisting of the node in which the
 * access is performed and the access itself.
 */
static na_pair find_unique_source(pdg::PDG* pdg, pdg::array *a)
{
        for (int i = 0; i < pdg->nodes.size(); ++i) {
                pdg::node *node = pdg->nodes[i];
                pdg::statement *s = pdg->nodes[i]->statement;
                for (int j = 0; j < s->accesses.size(); ++j) {
                        pdg::access *access = s->accesses[j];
                        if (access->array != a)
                                continue;
                        if (access->type != pdg::access::write)
                                continue;
                        return na_pair(node, access);
                }
        }

        assert(0);
}

/* Add a dependence from (source_node, source_access) to
 * (sink_node, sink_access) to pdg->dependences, for the
 * case where the array being accessed is marked uniquely_defined.
 *
 * Since the array is marked uniquely_defined, the value based
 * dependence is equal to the memory based dependence, so we
 * simply need to compose the access relations to obtain
 * the dependence relation.
 * This dependence relation is then specialized with respect to
 * the context and the iteration domain of the sink.
 */
static void add_unique_dep(PDG *pdg, pdg::node *source_node,
        pdg::access *source_access, pdg::node *sink_node,
        pdg::access *sink_access)
{
        pdg::dependence *d;
        isl_set *dom;
        isl_map *dep;
        isl_map *read;

        d = new pdg::dependence;
        d->array = source_access->array;
        d->type = pdg::dependence::flow;
        d->from = source_node;
        d->to = sink_node;
        d->from_access = source_access;
        d->to_access = sink_access;

        dep = source_access->map->get_isl_map();
        read = sink_access->map->get_isl_map();
        dep = isl_map_apply_range(dep, isl_map_reverse(read));
        dom = sink_node->source->get_isl_set();
        if (isl_set_is_wrapping(dom))
                dom = isl_map_domain(isl_set_unwrap(dom));
        dep = isl_map_intersect_range(dep, dom);
        dep = isl_map_intersect_params(dep, pdg->get_context_isl_set());
        d->relation = new pdg::IslMap(dep);

        pdg->dependences.push_back(d);
}

/* Find the flow dependences associated to the array "a", which is marked
 * uniquely_defined, and add them to pdg->dependences.
 *
 * First determine the unique source and then iterate through all the reads,
 * adding dependences from the unique source to each of the reads.
 */
static void find_unique_deps(PDG *pdg, pdg::array *a)
{
        na_pair na = find_unique_source(pdg, a);

        for (int i = 0; i < pdg->nodes.size(); ++i) {
                pdg::node *node = pdg->nodes[i];
                pdg::statement *s = pdg->nodes[i]->statement;
                for (int j = 0; j < s->accesses.size(); ++j) {
                        pdg::access *access = s->accesses[j];
                        if (access->array != a)
                                continue;
                        if (access->type != pdg::access::read)
                                continue;
                        add_unique_dep(pdg, na.node, na.access, node, access);
                }
        }
}

/* Find the dependence of type "t" associated to array "a" and add them
 * to pdg->dependences.
 *
 * If we are looking for flow dependences for an array that is marked
 * uniquely_defined, then we do not need to compute anything, but instead
 * can simply read off the dependences in find_unique_deps.
 */
void find_deps(PDG* pdg, pdg::array *a, type t)
{
    isl_ctx *ctx = pdg->get_isl_ctx();
    int nparam = pdg->params.size();
    int selfinput = 0;
    int reuse = 0;
    int firstuse = 0;
    add_dep_info info = { pdg, a, t };
    isl_set *context;
    bool need_parametrization;

    switch (t) {
    case flow:
        info.dtype = pdg::dependence::flow;
        firstuse = 1;
        break;
    case anti:
        info.dtype = pdg::dependence::anti;
        break;
    case data_reuse:
        info.dtype = pdg::dependence::reuse;
        reuse = 1;
        firstuse = 1;
        break;
    case reuse_pair:
        info.dtype = pdg::dependence::reuse_pair;
        break;
    case output:
        info.dtype = pdg::dependence::output;
        break;
    }

    a->analysis_performed.push_back(new pdg::dependence_type(info.dtype));

    if (t == flow && a->uniquely_defined) {
        find_unique_deps(pdg, a);
        return;
    }

    vector<na_pair> readers;
    vector<na_pair> writers;
    for (int i = 0; i < pdg->nodes.size(); ++i) {
        pdg::node *node = pdg->nodes[i];
        pdg::statement *s = pdg->nodes[i]->statement;
        for (int j = 0; j < s->accesses.size(); ++j) {
            pdg::access *access = s->accesses[j];
            if (access->array != a)
                continue;
            switch (t) {
            case flow:
            case data_reuse:
            case anti:
                if ((access->type == pdg::access::read) ^ (t == anti))
                    readers.push_back(na_pair(node, access));
                else
                    writers.push_back(na_pair(node, access));
                break;
            case output:
                if (access->type == pdg::access::read)
                    continue;
            case reuse_pair:
                readers.push_back(na_pair(node, access));
                writers.push_back(na_pair(node, access));
                break;
            }
        }
    }

    int maxsize = (selfinput || reuse) ? writers.size() + readers.size()
                                       : writers.size();
    context = pdg->get_context_isl_set();
    info.precedes_level = (isl_access_level_before)
            ((t == reuse_pair) ? precedes_level_accesses
                                : precedes_level_nodes);
    need_parametrization = any_filters(writers);
    for (int i = 0; i < writers.size(); ++i) {
        writers[i].map = convert_access(&writers[i]);
        writers[i].project_out_access_filters();
    }
    for (int i = 0; i < readers.size(); ++i) {
        readers[i].map = convert_access(&readers[i]);
        readers[i].map = isl_map_intersect_params(readers[i].map,
                                                isl_set_copy(context));
        readers[i].project_out_access_filters();
    }
    for (int i = 0; i < readers.size(); ++i) {
        isl_access_info *acc;
        int n_dep = pdg->dependences.size();

        acc = isl_access_info_alloc(isl_map_copy(readers[i].map),
                                    &readers[i], info.precedes_level, maxsize);
        if (need_parametrization)
                acc = isl_access_info_set_restrict(acc, &do_restrict, &info);
        for (int j = 0; j < writers.size(); ++j)
                acc = isl_access_info_add_source(acc,
                        isl_map_copy(writers[j].map), 1, &writers[j]);
        if (selfinput && writers.size()) {
            pdg::node *readnode = readers[i].node;
            for (int j = 0; j < readers.size(); ++j) {
                if (readers[j].node == readnode)
                    acc = isl_access_info_add_source(acc,
                        isl_map_copy(readers[j].map), 1, &readers[j]);
            }
        }
        if (reuse) {
            for (int j = 0; j < readers.size(); ++j)
                acc = isl_access_info_add_source(acc,
                    isl_map_copy(readers[j].map), 1, &readers[j]);
        }
        info.set_read_na(&readers[i]);
        isl_flow *deps = isl_access_info_compute_flow(acc);
        isl_flow_foreach(deps, add_dep, &info);
        isl_map *no_source;
        no_source = isl_flow_get_no_source(deps, 1);
        no_source = isl_map_from_range(isl_map_domain(no_source));
        no_source = simplify_controls(no_source, &info, NULL);
        if (!isl_map_plain_is_empty(no_source) && firstuse) {
            pdg::dependence *d = new pdg::dependence;
            d->array = a;
            d->type = pdg::dependence::uninitialized;
            d->to = readers[i].node;
            d->to_access = readers[i].access;
            if (d->to_access->extension) {
                d->extended_relation = new pdg::IslMap(isl_map_copy(no_source));
                no_source = isl_map_apply_range(no_source,
                            isl_map_reverse(
                                d->to_access->extension->get_isl_map(ctx)));
            }
            d->relation = new pdg::IslMap(isl_map_copy(no_source));
            pdg->dependences.push_back(d);
        }
        isl_map_free(no_source);
        isl_flow_free(deps);
        remove_unused_sources(pdg, n_dep, writers, &info);
    }
    isl_set_free(context);
}

/* Add a dependence from "write" to "a" to a->sources.
 */
static void add_unique_source(pdg::access *a, pdg::access *write)
{
        isl_map *dep;
        isl_map *read;

        dep = isl_map_range_map(write->map->get_isl_map());
        read = isl_map_range_map(a->map->get_isl_map());
        dep = isl_map_apply_range(dep, isl_map_reverse(read));

        a->sources.push_back(new pdg::IslMap(dep));
}

/* Look for the unique write access that writes to the array accessed
 * by "a" and then add a dependence from that write to a->sources.
 */
static void add_unique_source(pdg::PDG* pdg, pdg::access *a)
{
        na_pair na = find_unique_source(pdg, a->array);
        add_unique_source(a, na.access);
}

extern "C" {
        static isl_stat add_source(__isl_take isl_map *dep, int must,
                void *dep_user, void *user);
}

/* Add "dep" to a->sources, provided it is exact, and return isl_stat_ok.
 * Otherwise, return isl_stat_error.
 */
static isl_stat add_source(__isl_take isl_map *dep, int must, void *dep_user,
        void *user)
{
        bool has_controls;
        pdg::access *a = (pdg::access *) user;

        dep = remove_redundant_controls(dep, &has_controls);
        if (has_controls) {
                isl_map_free(dep);
                return isl_stat_error;
        }

        a->sources.push_back(new pdg::IslMap(dep));

        return isl_stat_ok;
}

extern "C" {
        static __isl_give isl_restriction *do_restrict_domain_map(
                __isl_keep isl_map *source_map, __isl_keep isl_set *sink,
                void *source_user, void *user);
}

/* Compute a restriction for the given map from sinks to potential sources
 * (sink2source).  We simply call compute_restriction to compute the
 * restriction.  This function is used from within find_sources,
 * which encodes the entire access relation into the domains of
 * the access relations passed to isl_access_info_compute_flow.
 * That is, the access relations passed to isl_access_info_compute_flow
 * are the result of applying isl_map_range_map to the original access
 * relations.  We therefore pass mappings that undo this encoding
 * to compute_restriction.
 */
static __isl_give isl_restriction *do_restrict_domain_map(
        __isl_keep isl_map *source_map, __isl_keep isl_set *sink,
        void *source_user, void *user)
{
        na_pair *source_na = (na_pair *) source_user;
        add_dep_info *info = (struct add_dep_info *) user;
        isl_space *space;
        isl_map *source_domain_map, *sink_domain_map;

        space = isl_space_range(isl_map_get_space(source_map));
        space = isl_space_unwrap(space);
        source_domain_map = isl_map_domain_map(isl_map_universe(space));
        space = isl_space_domain(isl_map_get_space(source_map));
        space = isl_space_unwrap(space);
        sink_domain_map = isl_map_domain_map(isl_map_universe(space));

        return compute_restriction(source_map, sink_domain_map,
                                source_domain_map, sink, source_na, info);
}

/* Find the sources of (read) access "a" in node "node".
 * If they are complete (no uninitialized accesses) and exact,
 * then put them in a->sources.  Otherwise, discard them.
 *
 * If the array is marked uniquely_defined, then we simply look
 * for the defining write in find_unique_source.
 *
 * Otherwise, we look for all writes that write to the same array,
 * perform dependence analysis and then check whether
 * the result is complete and exact.
 *
 * The sources record not only the node iteration, but also
 * the index of the array element.  We therefore apply
 * isl_map_range_map to the access relations, to obtain
 * a relation from the access (iteration -> element)
 * to the array element and feed that to the dependence analysis engine.
 */
void find_sources(pdg::PDG* pdg, pdg::node *node, pdg::access *a)
{
        isl_set *context;
        isl_map *no_source;
        isl_access_info *acc;
        isl_flow *deps;
        vector<na_pair> writers;
        na_pair reader(node, a);
        add_dep_info info = { pdg, a->array, flow, pdg::dependence::flow };

        if (a->array->uniquely_defined) {
                add_unique_source(pdg, a);
                return;
        }

        info.precedes_level = (isl_access_level_before) precedes_level_nodes;

        for (int i = 0; i < pdg->nodes.size(); ++i) {
                pdg::node *node = pdg->nodes[i];
                pdg::statement *s = pdg->nodes[i]->statement;
                for (int j = 0; j < s->accesses.size(); ++j) {
                        pdg::access *access = s->accesses[j];
                        if (access->array != a->array)
                                continue;
                        if (access->type != pdg::access::write)
                                continue;
                        writers.push_back(na_pair(node, access));
                }
        }

        context = pdg->get_context_isl_set();
        for (int i = 0; i < writers.size(); ++i) {
                writers[i].projected_map = convert_access(&writers[i]);
                writers[i].map = isl_map_copy(writers[i].projected_map);
                writers[i].map = isl_map_range_map(writers[i].map);
                writers[i].project_out_access_filters();
        }
        reader.projected_map = convert_access(&reader);
        reader.projected_map = isl_map_intersect_params(reader.projected_map,
                                                        context);
        reader.map = isl_map_range_map(isl_map_copy(reader.projected_map));
        reader.project_out_access_filters();

        acc = isl_access_info_alloc(isl_map_copy(reader.map),
                                &reader, info.precedes_level, writers.size());
        for (int j = 0; j < writers.size(); ++j)
                acc = isl_access_info_add_source(acc,
                        isl_map_copy(writers[j].map), 1, &writers[j]);

        if (any_filters(writers))
                acc = isl_access_info_set_restrict(acc,
                                                &do_restrict_domain_map, &info);

        info.set_read_na(&reader);
        deps = isl_access_info_compute_flow(acc);
        no_source = isl_flow_get_no_source(deps, 1);
        if (isl_map_plain_is_empty(no_source)) {
                if (isl_flow_foreach(deps, add_source, a) < 0) {
                        for (int i = 0; i < a->sources.size(); ++i)
                                delete a->sources[i];
                        a->sources.resize(0);
                }
        }
        isl_map_free(no_source);
        isl_flow_free(deps);
}

/* Find the source (if possible) of the filter "coa" in node "node".
 * We assume that the filter is an access rather than a function call.
 */
static void find_sources(pdg::PDG *pdg, pdg::node *node,
        pdg::call_or_access *coa)
{
        pdg::access *access;

        assert(coa->type == pdg::call_or_access::t_access);
        access = coa->access;

        find_sources(pdg, node, access);
}

/* Compute the sources (if possible) for all the filters in all the
 * nodes and accesses in "pdg".
 */
void compute_filter_sources(pdg::PDG *pdg)
{
        for (int i = 0; i < pdg->nodes.size(); ++i) {
                pdg::node *node = pdg->nodes[i];
                pdg::statement *s = pdg->nodes[i]->statement;
                int n_filter = node->filters.size();

                for (int j = 0; j < n_filter; ++j)
                        find_sources(pdg, node, node->filters[j]);

                for (int j = 0; j < s->accesses.size(); ++j) {
                        pdg::access *access = s->accesses[j];

                        for (int k = 0; k < access->nested.size(); ++k)
                                find_sources(pdg, node, access->nested[k]);
                }
        }
}

static int precedes_level_nodes(na_pair *first, na_pair *second)
{
    int d = 0;
    int cmp = 0;
    for (int i = 0; i < first->node->prefix.size(); ++i) {
        if (i >= second->node->prefix.size())
            break;
        cmp = first->node->prefix[i] - second->node->prefix[i];
        if (cmp)
            break;
        if (first->node->prefix[i] == -1)
            ++d;
    }
    return 2*d + (cmp<0);
}

static int precedes_level_accesses(na_pair *first, na_pair *second)
{
    int cmp = 0;
    int d = 0;
    for (int i = 0; i < first->node->prefix.size(); ++i) {
        if (i >= second->node->prefix.size())
            break;
        cmp = first->node->prefix[i] - second->node->prefix[i];
        if (cmp)
            break;
        if (first->node->prefix[i] == -1)
            ++d;
    }
    /* same node; now compare accesses */
    if (!cmp)
        cmp = first->access->nr - second->access->nr;
    return 2*d + (cmp<0);
}

}