pn: allow specification of isl options

[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/constraint.h"
#include <isl/set.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;
        assert(isl_map_n_out(map) < array->dims.size());
        unsigned s_dim = array->dims.size() - isl_map_n_out(map);
        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 *dim = isl_map_get_space(map);
        dim = isl_space_drop_dims(dim, isl_dim_in,
                                0, isl_space_dim(dim, isl_dim_in));
        dim = isl_space_drop_dims(dim, isl_dim_out,
                                0, isl_space_dim(dim, isl_dim_out));
        dim = isl_space_add_dims(dim, isl_dim_in, s_dim);
        dim = isl_space_add_dims(dim, isl_dim_out, s_dim);
        isl_basic_map *id = isl_basic_map_identity(isl_space_copy(dim));
        isl_local_space *ls = isl_local_space_from_space(dim);
        isl_int v;
        isl_int_init(v);
        for (int i = 0; i < s_dim; ++i) {
                isl_constraint *c;
                c = isl_inequality_alloc(isl_local_space_copy(ls));
                isl_int_set_si(v, 1);
                isl_constraint_set_coefficient(c, isl_dim_out, 0, v);
                id = isl_basic_map_add_constraint(id, c);
                c = isl_inequality_alloc(isl_local_space_copy(ls));
                isl_int_set_si(v, -1);
                isl_constraint_set_coefficient(c, isl_dim_out, 0, v);
                isl_int_set_si(v, array->dims[isl_map_n_out(map)+i]-1);
                isl_constraint_set_constant(c, v);
                id = isl_basic_map_add_constraint(id, c);
        }
        isl_local_space_free(ls);
        isl_int_clear(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_n_in(map) != isl_set_dim(dom, isl_dim_set))
                assert(0);
        map = isl_map_intersect_domain(map, dom);
        if (isl_map_n_out(map) != 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;

        /* 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);

        ~add_dep_info();
};

/* 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);
}

/* 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;
}

/* Add parameters to "set" identifying the last iteration of the access
 * identified by "na".
 *
 * In particular, we add a parameter
 *
 *      __last_<stmt>_<access_nr>_valid = 1
 *
 * and parameters
 *
 *      __last_<stmt>_<access_nr>_<i> = it_<i>
 *
 * with i ranging over the iterators 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 parameter, 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 = valid_bit_id(ctx, na);
        pos = isl_set_find_dim_by_id(set, isl_dim_param, id);
        if (pos < 0) {
                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);
        } else
                isl_id_free(id);
        set = isl_set_fix_si(set, isl_dim_param, pos, 1);

        for (int i = 0; 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);

                pos = isl_set_find_dim_by_name(set, isl_dim_param, name);
                if (pos < 0) {
                        id = isl_id_alloc(ctx, name, NULL);
                        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);
                }

                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 int 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 int 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 0;
}

/* 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.
 */
void add_dep_info::set_read_na(na_pair *read_na)
{
        used.clear();
        clear_mem_dep();
        read_na_pair = read_na;
}

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

/* Is parameter "i" of "space" such that it is called __last_*
 * and has a non-NULL user field?
 * 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 *prefix = "__last_";
        size_t prefix_len = strlen(prefix);
        const char *name;
        isl_id *id;
        na_pair *write_na;

        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;
        id = isl_space_get_dim_id(space, isl_dim_param, i);
        write_na = (na_pair *) isl_id_get_user(id);
        isl_id_free(id);
        if (!write_na)
                return false;

        return true;
}

/* 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;
}

/* 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.
 *
 * Call the given maps A and B.
 *
 * If the range spaces of the two maps are different, then the filter
 * is obviously never the same.
 *
 * 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 need to impose that
 *
 *      \forall j : i -> j \in A => j -> k \in B
 *
 * In other words
 *
 *      \not \exists j : i -> j \in A and j -> k \not\in B
 *
 * We therefore need to remove the elements i -> k that belong to the relation
 *
 *      A . (\bar{B})^-1
 *
 * where \bar{B} is the complement of B.
 *
 * Note that if A is single-valued, then we don't need to perform
 * this computation.  Each iteration accesses only a single element
 * and we already know that at least one element is accessed in common,
 * so this single element must be a subset of the elements accessed
 * through B.
 *
 * 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,
 * let A be of the form { S[..] -> F[..] }.  We first construct the
 * mapping { S[..] -> [F[..] -> F[..]] } and remove from the range
 * of this mapping identical pairs of F indices.  What is left
 * is a mapping from iterations to pairs of different F indices.
 * Any element in the domain of this mapping accesses more than one
 * element and is therefore removed.
 */
static __isl_give isl_map *compute_common_filter(__isl_keep isl_map *map1,
        bool subset1, __isl_keep isl_map *map2, bool subset2)
{
        isl_space *space1, *space2;
        isl_map *reverse, *common, *bad;
        int equal;

        space1 = isl_space_range(isl_map_get_space(map1));
        space2 = isl_space_range(isl_map_get_space(map2));
        equal = isl_space_is_equal(space1, space2);
        isl_space_free(space1);
        isl_space_free(space2);
        if (!equal) {
                space1 = isl_space_domain(isl_map_get_space(map1));
                space2 = isl_space_domain(isl_map_get_space(map2));
                space1 = isl_space_map_from_domain_and_range(space1, space2);
                return isl_map_empty(space1);
        }

        reverse = isl_map_reverse(isl_map_copy(map2));
        common = isl_map_apply_range(isl_map_copy(map1), isl_map_copy(reverse));
        if (subset1 && !subset2 && !isl_map_plain_is_single_valued(map1)) {
                bad = isl_map_apply_range(isl_map_copy(map1),
                                    isl_map_complement(isl_map_copy(reverse)));
                common = isl_map_subtract(common, bad);
        }
        if (subset2 && !subset1 && !isl_map_plain_is_single_valued(map2)) {
                bad = isl_map_apply_range(
                                    isl_map_complement(isl_map_copy(map1)),
                                    isl_map_copy(reverse));
                common = isl_map_subtract(common, bad);
        }
        if (subset1 && subset2 && !isl_map_plain_is_single_valued(map1)) {
                bad = isl_map_range_product(isl_map_copy(map1),
                                            isl_map_copy(map1));
                space1 = isl_space_range(isl_map_get_space(bad));
                space1 = isl_space_unwrap(space1);
                bad = isl_map_subtract_range(bad,
                        isl_map_wrap(isl_map_universe(space1)));
                common = isl_map_subtract_domain(common, isl_map_domain(bad));
        }
        if (subset1 && subset2 && !isl_map_plain_is_single_valued(map2)) {
                bad = isl_map_range_product(isl_map_copy(map2),
                                            isl_map_copy(map2));
                space1 = isl_space_range(isl_map_get_space(bad));
                space1 = isl_space_unwrap(space1);
                bad = isl_map_subtract_range(bad,
                        isl_map_wrap(isl_map_universe(space1)));
                common = isl_map_subtract_range(common, isl_map_domain(bad));
        }
        isl_map_free(reverse);

        return common;
}

/* Assuming "coa" is a (read) access, construct a map from the domain
 * of the access relation to the access relation of the corresponding
 * write.  If we are unable to determine the corresponding write, then
 * return the empty map.
 */
static __isl_give isl_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;
}

/* 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 valid bit (__last_<na>_valid).
 *
 * In particlar, if valid is set, we return the set
 *
 *      { S[i] : i = __last_<na> and __last_<na>_valid >= 1 }
 *
 * If valid is not set, we return the set
 *
 *      { S[i] : (i >> __last_<na> and __last_<na>_valid >= 1) or
 *               __last_<na>_valid <= 0 }
 *
 * That is, the iterations after the last if there is a last iteration
 * or just any iteration if there is no 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 = 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, 0);
        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.
 *
 * 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)
{
        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(space1);
        same_value = isl_map_universe(space2);

        for (int i = 0; i < node1->filters.size(); ++i) {
                isl_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]);
                if (isl_map_is_empty(map_i)) {
                        isl_map_free(map_i);
                        continue;
                }

                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_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);
                        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_map_free(map_j);
                }

                isl_map_free(map_i);
        }

        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",
 * 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, then the filter values should
 * lie outside of the possible values for any potential source iteration.
 * Note however, that it 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))
                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 = isl_map_intersect_domain(mem1, source);

        mem2 = info->get_mem_dep(write2_na);
        source = source_iterations(write2_na, id2, valid2, info);
        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".
 */
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",
 * 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 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))
                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 = 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".
 */
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;
}

/* Simplify the constraints on the parameters introduced
 * in add_parametrization().
 * We first remove all those parameters that turn out not be needed.
 * If there are any of those parameters left, then we 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 = remove_redundant_controls(dep, has_controls);
        if (*has_controls)
                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 int 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 0;
        }

        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 0;
}

/* 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_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_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_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_map *index)
{
        for (int i = 0; i < filter->index.size(); ++i) {
                int equal;

                equal = isl_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".
 */
static struct da_filter *da_filter_add(struct da_filter *filter,
        __isl_take isl_map *index)
{
        isl_id *id;
        isl_space *space;

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

        filter->value = isl_map_add_dims(filter->value, isl_dim_out, 1);
        space = isl_space_range(isl_map_get_space(index));
        if (isl_space_is_wrapping(space))
                space = isl_space_range(isl_space_unwrap(space));
        id = isl_space_get_tuple_id(space, isl_dim_set);
        isl_space_free(space);
        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_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"?
 */
static int is_total_filter(__isl_keep isl_map *map, __isl_keep isl_set *domain)
{
        isl_set *map_domain;
        int total;

        map_domain = isl_map_domain(isl_map_copy(map));
        total = isl_set_is_subset(domain, map_domain);
        isl_set_free(map_domain);

        return total;
}

/* Does "node" have any total filters?
 */
static int any_total_filter(pdg::node *node, __isl_keep isl_set *domain)
{
        for (int i = 0; i < node->filters.size(); ++i) {
                isl_map *map = extract_access_map(node->filters[i]);
                int total;

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

                if (total)
                        return total;
        }

        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.  If we are unable to determine this
 * write access, then we continue with the next source filter.
 * Otherwise 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;
        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));
        space = isl_space_map_from_domain_and_range(
                        isl_space_range(isl_map_get_space(sink_value)), 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_map *map_i;
                map_i = extract_access_map(source_na->node->filters[i]);
                if (isl_map_is_empty(map_i)) {
                        isl_map_free(map_i);
                        continue;
                }
                map_i = isl_map_apply_range(isl_map_copy(source_map), map_i);
                for (int j = 0; j < filter->index.size(); ++j) {
                        if (isl_map_is_subset(map_i, filter->index[j]))
                                filter_map = isl_map_equate(filter_map,
                                        isl_dim_in, j, isl_dim_out, i);
                }
                isl_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;
}

/* 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) \cap {i -> \cap_{k \in map1(i)} map2(k)}
 *
 * If map1 is single valued, then we can do a simple join.
 */
static __isl_give isl_map *join_non_empty_subset(__isl_take isl_map *map1,
        __isl_take isl_map *map2)
{
        isl_map *res;

        if (isl_map_plain_is_single_valued(map1))
                return isl_map_apply_range(map1, map2);

        res = isl_map_apply_range(isl_map_copy(map1), isl_map_copy(map2));
        map2 = isl_map_complement(map2);
        map1 = isl_map_apply_range(map1, map2);
        res = isl_map_subtract(res, map1);

        return res;
}

/* 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.
 * 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.
 *
 * Project out those dimensions that correspond to filters with
 * an access function that is not total, i.e., one that does not
 * access any element for some of the non-executed iterations.
 * If there are no total 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 [K -> F] -> T is a non-empty subset 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_set *sink, na_pair *source_na, add_dep_info *info)
{
        isl_space *space;
        isl_set *source_domain;
        isl_map *mem;
        isl_map *may_run, *not_run;
        isl_map *source_value;
        std::vector<isl_map *> index;
        isl_map *index_product;
        isl_map *map;
        isl_map *filter_map;
        isl_set *unaffected_sink;
        isl_map *unaffected_value;
        isl_set *domain;

        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 = remove_all_controls(may_run);
        mem = info->get_mem_dep(source_na);
        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);

        domain = isl_map_domain(isl_map_copy(not_run));

        if (isl_map_is_empty(not_run) ||
            !any_total_filter(source_na->node, domain)) {
                isl_set_free(domain);
                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_map *map = extract_access_map(source_na->node->filters[i]);
                int sv;

                sv = is_total_filter(map, domain);
                if (sv) {
                        index.push_back(map);
                } else {
                        isl_map_free(map);
                        source_value = isl_map_project_out(source_value,
                                                isl_dim_out, index.size(), 1);
                }
        }

        index_product = isl_map_copy(index[0]);
        for (int i = 1; i < index.size(); ++i)
                index_product = isl_map_range_product(index_product,
                                                        isl_map_copy(index[1]));
        index_product = isl_map_reverse(index_product);
        index_product = isl_map_domain_product(isl_map_copy(not_run),
                                                index_product);

        source_value = join_non_empty_subset(index_product, source_value);
        space = isl_space_domain(isl_map_get_space(source_value));
        map = isl_map_domain_map(isl_map_universe(isl_space_unwrap(space)));
        source_value = isl_map_apply_domain(source_value, map);

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

        space = isl_space_map_from_domain_and_range(
                        isl_space_range(isl_map_get_space(source_value)),
                        isl_space_range(isl_map_get_space(filter->value)));
        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_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_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_map_free(index[i]);
        isl_set_free(domain);
        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.
 * 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_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_valid_bit(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, sink, na, info);
        }
        isl_space_free(space);

        return filter;
}

#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.
 *
 * If there are no filters on the source, then obviously the source
 * is always executed.
 *
 * 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.
 * If, moreover, the relation is empty then the filters cannot have any value,
 * meaning that the sources cannot have executed.
 *
 * 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 sources are necessarily executed, 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;

        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_empty(sink_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, sink, info);

        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_empty(sink_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 (__last_<na>_valid = 0) or that there is a last
 * iteration (__last_<na>_valid = 1) 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 pos;

        id = valid_bit_id(isl_set_get_ctx(set), na);
        pos = isl_set_find_dim_by_id(set, isl_dim_param, id);

        if (pos < 0) {
                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);
        } else
                isl_id_free(id);

        valid = isl_set_copy(set);
        invalid = set;

        valid = isl_set_fix_si(valid, isl_dim_param, pos, 1);

        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_fix_si(invalid, isl_dim_param, pos, 0);

        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 and compute the
 * last of these definitely executed source iterations.
 *
 * For those sink iterations that have a corresponding last 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);
        definite_source_map = isl_map_lexmax(definite_source_map);

        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_le_first(space, depth);
        map_after = isl_map_intersect_range(map_after, dom);

        definite_source_map = isl_map_apply_range(definite_source_map,
                                                map_after);

        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 = valid_bit_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_map_get_space(source_map);
        n_param = isl_space_dim(space, isl_dim_param);
        for (int i = 0; i < n_param; ++i) {
                isl_id *other_id;

                if (!is_valid_bit(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 (*_valid = 0)
 * or it is executed (*_valid = 1) 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.
 */
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 = 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);
}

/* 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.
 */
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));

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

        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.
 * 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);
        }

        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;
                }
        }
}

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));

    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]);
    }
    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_fast_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 find_unique_source(pdg::PDG* pdg, pdg::access *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->array)
                                continue;
                        if (access->type != pdg::access::write)
                                continue;
                        add_unique_source(a, access);
                        return;
                }
        }

        assert(0);
}

extern "C" {
        static int 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 0.
 * Otherwise, return -1.
 */
static int 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 -1;
        }

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

        return 0;
}

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) {
                find_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);
        }
        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);
}

}