update isl for change in lexicographic optimization

[isa.git] / dependence_graph.cc

#include <isl/val.h>
#include <isl/space.h>
#include <isl/aff.h>
#include <isl/set.h>
#include <isl/map.h>
#include <isl/printer.h>
#include <isa/pdg.h>
#include "dependence_graph.h"

void edge::dump(void)
{
        isl_ctx *ctx;
        isl_printer *prn;
        const char *s;
        switch (type) {
        case normal:            s = ""; break;
        case expansion:         s = "expansion "; break;
        }
        fprintf(stderr, "\tpos: %d, source: %d,%s %s[%p]\n",
                        pos, source->location, source->operation, s, source);
        ctx = isl_map_get_ctx(relation);
        prn = isl_printer_to_file(ctx, stderr);
        prn = isl_printer_set_indent(prn, 8);
        prn = isl_printer_print_map(prn, relation);
        prn = isl_printer_end_line(prn);
        isl_printer_free(prn);
}

void computation::dump(void)
{
        isl_ctx *ctx;
        isl_printer *prn;
        const char *s;
        if (input)
                s = "input ";
        else if (is_expanded)
                s = "is_expanded ";
        else
                s = "";
        fprintf(stderr, "op: %d,%s/%d %s[%p]\n",
                        location, operation, arity, s, this);
        ctx = isl_set_get_ctx(domain);
        prn = isl_printer_to_file(ctx, stderr);
        prn = isl_printer_print_set(prn, domain);
        prn = isl_printer_end_line(prn);
        isl_printer_free(prn);

        for (int i = 0; i < edges.size(); ++i) {
                fprintf(stderr, "edge %d:\n", i);
                edges[i]->dump();
        }
}

void dependence_graph::dump(void)
{
        for (int i = 0; i < vertices.size(); ++i)
                vertices[i]->dump();
        out->dump();
}

static void add_input_edges(isl_ctx *ctx, int pos, pdg::dependence *dep,
                        computation *comp, computation *source)
{
        pdg::access *a = dep->to_access;
        assert(source);
        isl_map *map = a->map->get_isl_map(ctx);
        isl_map *relation = dep->relation->get_isl_map(ctx);
        isl_set *domain = isl_map_range(relation);
        map = isl_map_intersect_domain(map, domain);

        edge *e = new edge;
        e->source = source;
        e->pos = pos;
        e->relation = map;
        comp->edges.push_back(e);
}

/* For each basic map in relation, add an edge to comp with
 * given source and pos.
 */
static void add_split_edge(computation *comp, computation *source,
        int pos, isl_map *relation, std::vector<computation **> *missing,
        enum edge::type type = edge::normal)
{
        edge *e = new edge;
        e->source = source;
        if (!source) {
                assert(missing);
                missing->push_back(&e->source);
        }
        e->pos = pos;
        e->relation = relation;
        e->type = type;
        comp->edges.push_back(e);
}

/* For each (flow) dependence with write access a, and an
 * edge to computation comp.
 * Note that in the pdg, the dependences have as "input" the source (read)
 * and as "output" the sink (write), while in the dependence_graph,
 * we want edges that point from sink to source.
 *
 * If the (read) access has an extension, i.e., if the access
 * reads a row or an even larger chunk of an array, then we introduces
 * an extra "expansion" edge to an extra node that corresponds to the
 * reading of the whole array.  The iteration domain of this extra node
 * has extra dimensions, as many as were missing from the access.
 * The extra edge maps each iteration of the statement in which
 * the read appears to all iterations in the extended iteration domain
 * with the same values for the first iterators.
 * The remaining iterators range over the elements of the array chunk.
 * All dependences in which the given access is involved are attached
 * to the new node, where we need to be careful to use the extended_relation
 * instead of the projected "standard" relations.
 * The operation associated with the new node is the
 * empty string since the new node can be treated as a copy operation.
 * The edge connecting the original computation and the new node
 * is marked as an "expansion" edge because it needs to be treated
 * in a special way.
 */
static void add_dep_edges(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                int pos, pdg::access *a, computation *comp)
{
        isl_ctx *ctx = pdg->get_isl_ctx();
        if (a->extension) {
                isl_map *ext = a->extension->get_isl_map(ctx);

                computation *child = new computation();
                child->domain = isl_set_apply(isl_set_copy(comp->domain),
                                                isl_map_copy(ext));
                child->operation = strdup("");
                child->arity = 1;
                child->location = comp->location;
                dg->vertices.push_back(child);

                add_split_edge(comp, child, pos, ext, NULL, edge::expansion);
                comp->is_expanded = 1;
                comp = child;
                pos = 0;
        }
        for (int i = 0; i < pdg->dependences.size(); ++i) {
                pdg::dependence *dep = pdg->dependences[i];
                if (dep->to_access != a)
                        continue;

                if (dep->type == pdg::dependence::uninitialized) {
                        add_input_edges(ctx, pos, dep, comp, a2c[dep->array]);
                        continue;
                }
                assert(dep->type == pdg::dependence::flow);

                isl_map *relation;
                if (!dep->extended_relation)
                        relation = dep->relation->get_isl_map(ctx);
                else
                        relation = dep->extended_relation->get_isl_map(ctx);
                relation = isl_map_reverse(relation);
                add_split_edge(comp, s2c[dep->from_access], pos, relation,
                                &missing[dep->from_access]);
        }
}

/* Each access without an associated array corresponds to a read
 * of an (affine combination of) iterator value(s).
 * The expression in terms of the iterators is given by the access map.
 * For each of these accesses, we add an edge to a special
 * "input computation" that is not associated to any particular
 * (input) array.  The map of the edge is the access map,
 * with the domain intersected with the domain of the current computation.
 * This ensures that corresponding reads are considered equal iff the
 * value they read is the same for both programs.
 *
 * In the special case where the affine expression is a mere integer constant
 * that appears as an argument to a #test operation, the domain of "comp"
 * may reference an argument that may not be present in "a" if it was
 * created by c2pdg from an integer expression.
 * In this case, we not only need to intersect the domain of
 * the access relation, but also introduce a reference to the argument.
 */
static void add_iterator_edges(pdg::PDG *pdg,
                            std::map<pdg::array *, computation *> &a2c,
                            int pos,
                            pdg::access *a, computation *comp)
{
        isl_ctx *ctx = pdg->get_isl_ctx();
        computation *Z = a2c[NULL];
        isl_map *map = a->map->get_isl_map(ctx);
        isl_set *domain = isl_set_copy(comp->domain);
        if (isl_set_is_wrapping(domain) && !isl_map_domain_is_wrapping(map)) {
                isl_map *domain_map;
                domain_map = isl_map_domain_map(isl_set_unwrap(domain));
                map = isl_map_apply_range(domain_map, map);
        } else
                map = isl_map_intersect_domain(map, domain);
        add_split_edge(comp, Z, pos, map, NULL);
}

static void add_access_edges(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                int pos, pdg::access *a, computation *comp)
{
        if (a->array)
                add_dep_edges(dg, pdg, s2c, a2c, missing, pos, a, comp);
        else
                add_iterator_edges(pdg, a2c, pos, a, comp);
}

static int add_edges_at_pos(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                int pos, pdg::call_or_access *coa, computation *comp);

/* Add a new copy computation to dependence graph dg with given
 * domain and location and return a pointer to the copy computation.
 */
static computation *add_copy_computation(dependence_graph *dg,
                isl_set *domain, int location)
{
        computation *copy;
        copy = new computation();
        copy->domain = domain;
        copy->operation = strdup("");
        copy->arity = 1;
        copy->location = location;
        dg->vertices.push_back(copy);
        return copy;
}

/* Create a "#NEST computation, with n + 1 arguments, where n is the
 * number of nested accesses.
 * The first n arguments correspond to the nested accesses and
 * edges for these arguments are added as usual.
 * Note that if the access "a" appears in a "#test" function,
 * then some of the nested accesses will already have been
 * added by an outer #NEST computation and should therefore be skipped.
 * The final argument corresponds to the expansion and for this position
 * a single expansion edge is added to the given copy computation.
 * The mapping on the expansion edge is the identity
 * mapping on the domain, with the range extended by one dimension for
 * each of nested accesses.  The range of values in each of these dimensions
 * is given by the value_bounds on the corresponding array.
 * The iteration domain of the given copy computation is adjusted accordingly.
 */
static computation *add_nest_computation(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                pdg::access *a, computation *copy)
{
        isl_ctx *ctx = pdg->get_isl_ctx();
        computation *nest = new computation();
        nest->domain = isl_set_copy(copy->domain);
        nest->operation = strdup("#NEST");
        nest->arity = a->nested.size() + 1;
        nest->location = copy->location;
        dg->vertices.push_back(nest);

        int nested_dim = a->map->input() - isl_set_dim(copy->domain, isl_dim_set);
        assert(nested_dim > 0 && nested_dim <= a->nested.size());
        int skip_dim = a->nested.size() - nested_dim;

        isl_space *dim;
        isl_set *bounds;
        isl_map *ext;
        int pos = 0;

        dim = isl_set_get_space(copy->domain);
        dim = isl_space_drop_dims(dim, isl_dim_set,
                                        0, isl_space_dim(dim, isl_dim_set));
        bounds = isl_set_universe(dim);

        for (int i = skip_dim; i < a->nested.size(); ++i) {
                pdg::call_or_access *coa = a->nested[i];

                if (add_edges_at_pos(dg, pdg, s2c, a2c,
                                     missing, pos, coa, nest))
                        ++pos;

                assert(coa->type == pdg::call_or_access::t_access);
                pdg::access *nested = coa->access;
                assert(nested->array->value_bounds);
                isl_set *bounds_i;
                bounds_i = nested->array->value_bounds->get_isl_set(ctx);

                bounds = isl_set_flat_product(bounds, bounds_i);
        }

        dim = isl_set_get_space(copy->domain);
        copy->domain = isl_set_product(copy->domain, isl_set_copy(bounds));
        bounds = isl_set_product(isl_set_universe(dim), bounds);
        ext = isl_map_reverse(isl_map_domain_map(isl_set_unwrap(bounds)));

        add_split_edge(nest, copy, pos, ext, NULL, edge::expansion);
        nest->is_expanded = 1;

        return nest;
}

/* Add an edge to "comp" at position "pos" according to access "a",
 * where "a" is has nested accesses.
 * We first create a copy computation where the actual access
 * will be attached.  Then we create a "#NEST" computation that
 * sits between "comp" and the copy computation.
 */
static void add_nested_access_edge(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                int pos, pdg::access *a, computation *comp)
{
        computation *copy;
        copy = add_copy_computation(dg, isl_set_copy(comp->domain),
                                    comp->location);

        computation *nest;
        nest = add_nest_computation(dg, pdg, s2c, a2c, missing, a, copy);

        edge *e = new edge;
        e->pos = pos;
        e->source = nest;
        isl_space *dim = isl_space_map_from_set(isl_set_get_space(comp->domain));
        e->relation = isl_map_identity(dim);
        comp->edges.push_back(e);

        add_access_edges(dg, pdg, s2c, a2c, missing, 0, a, copy);
}

/* "Invocations" of integer constants are replaced by accesses
 * to the "input computation", so we can treat them like and combine them
 * with iterator accesses.
 */
static void add_constant_edge(pdg::PDG *pdg,
                            std::map<pdg::array *, computation *> &a2c,
                            int pos, long cst, computation *comp)
{
        isl_ctx *ctx = pdg->get_isl_ctx();
        computation *Z = a2c[NULL];

        isl_map *map = isl_map_from_range(isl_set_copy(comp->domain));
        map = isl_map_reverse(map);
        map = isl_map_add_dims(map, isl_dim_out, 1);
        map = isl_map_fix_si(map, isl_dim_out, 0, cst);
        add_split_edge(comp, Z, pos, map, NULL);
}

static void add_edges(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                pdg::function_call *call, computation *comp);

/* Scale the single output dimension of bmap by a factor of f */
static __isl_give isl_map *scale(__isl_take isl_map *map, __isl_take isl_val *f)
{
        isl_multi_aff *scaling;
        isl_space *space;

        space = isl_map_get_space(map);
        space = isl_space_range(space);
        space = isl_space_map_from_set(space);
        scaling = isl_multi_aff_identity(space);
        scaling = isl_multi_aff_scale_val(scaling, f);

        map = isl_map_apply_range(map, isl_map_from_multi_aff(scaling));

        return map;
}

/* Compute the mapping along the given edge to the Z computation,
 * where the Z computation may be either the source of the edge
 * or it may be the source of the single edge leaving the copy
 * computation that is the source of the given edge.
 */
static isl_map *map_to_Z(edge *edge, computation *Z)
{
        if (edge->source == Z)
                return isl_map_copy(edge->relation);
        assert(edge->source->is_copy());
        assert(edge->source->edges.size() == 1);
        assert(edge->source->edges[0]->source == Z);
        return isl_map_apply_range(
                isl_map_copy(edge->relation),
                isl_map_copy(edge->source->edges[0]->relation));
}

static isl_map *affine_expression_div(computation *comp, computation *Z)
{
        isl_map *map = NULL;
        isl_val *v;
        isl_map *rel = map_to_Z(comp->edges[1], Z);

        v = isl_map_plain_get_val_if_fixed(rel, isl_dim_out, 0);
        if (isl_val_is_int(v)) {
                map = isl_map_floordiv_val(
                    isl_map_copy(comp->edges[0]->relation), v);
        } else
                isl_val_free(v);
        isl_map_free(rel);
        return map;
}

/* Construct a map for CLooG's ceild operation.
 *
 * We simply apply the identity
 *
 *      ceil(a/b) = - floor(-a/b)
 */
static isl_map *affine_expression_ceildiv(computation *comp, computation *Z)
{
        isl_map *map = NULL;
        isl_map *map2;
        isl_val *v;

        map2 = map_to_Z(comp->edges[1], Z);
        v = isl_map_plain_get_val_if_fixed(map2, isl_dim_out, 0);
        if (isl_val_is_int(v)) {
                map = isl_map_neg(map_to_Z(comp->edges[0], Z));
                map = isl_map_neg(isl_map_floordiv_val(map, v));
        } else
                isl_val_free(v);

        isl_map_free(map2);
        return map;
}

/* Construct a mapping for the maximum of two mappings.
 *
 * In particular, given two mappings in the arguments of "comp"
 *
 *      [input] -> [a]
 *      [input] -> [b]
 *
 * construct
 *
 *      [input] -> [a] : a >= b; [input] -> [b] : a < b
 *
 * We first construct
 *
 *      [input] -> [[a] -> [b]]
 *
 * intersect the range with [a] -> [b] : a >= b and [a] -> [b] : a < b
 * and then map the range to its domain and range in the respective results.
 */
static isl_map *affine_expression_max(computation *comp, computation *Z)
{
        isl_space *dim;
        isl_map *map, *ge, *lt;
        isl_map *map1, *map2;

        map1 = map_to_Z(comp->edges[0], Z);
        map2 = map_to_Z(comp->edges[1], Z);

        dim = isl_map_get_space(map1);
        dim = isl_space_range(dim);

        map = isl_map_range_product(map1, map2);

        ge = isl_map_lex_ge(isl_space_copy(dim));
        lt = isl_map_lex_lt(isl_space_copy(dim));

        ge = isl_map_intersect_range(isl_map_copy(map), isl_map_wrap(ge));
        lt = isl_map_intersect_range(map, isl_map_wrap(lt));

        dim = isl_space_map_from_set(dim);
        map = isl_map_universe(isl_space_copy(dim));
        map = isl_map_domain_map(map);
        ge = isl_map_apply_range(ge, map);

        map = isl_map_universe(dim);
        map = isl_map_range_map(map);
        lt = isl_map_apply_range(lt, map);

        map = isl_map_union(ge, lt);

        return map;
}

/* Given two maps A -> f(A) and B -> g(B), construct a map
 * A -> f(A) * g(B) or B -> f(A) * g(B), if at least one of
 * "rel0" or "rel1" maps to a constant value.
 * Otherwise, return NULL.
 */
static __isl_give isl_map *map_mul(__isl_take isl_map *rel0,
        __isl_take isl_map *rel1)
{
        isl_val *v;

        v = isl_map_plain_get_val_if_fixed(rel0, isl_dim_out, 0);
        if (isl_val_is_int(v)) {
                isl_map_free(rel0);
                return scale(rel1, v);
        }
        isl_val_free(v);

        v = isl_map_plain_get_val_if_fixed(rel1, isl_dim_out, 0);
        if (isl_val_is_int(v)) {
                isl_map_free(rel1);
                return scale(rel0, v);
        }
        isl_val_free(v);

        isl_map_free(rel0);
        isl_map_free(rel1);

        return NULL;
}

/* If comp represents an operation on affine expressions that can
 * be represented by a single affine expression, then return a mapping
 * encoding this single affine expression.  Otherwise return NULL.
 */
static isl_map *affine_expression(computation *comp, computation *Z)
{
        if (comp->arity != comp->edges.size())
                return NULL;
        for (int i = 0; i < comp->edges.size(); ++i) {
                computation *source = comp->edges[i]->source;
                if (source == Z)
                        continue;
                if (source && source->is_copy() && source->edges.size() == 1 &&
                    source->edges[0]->source == Z)
                        continue;
                return NULL;
        }
        if (!strcmp(comp->operation, "+") && comp->arity == 2)
                return isl_map_sum(map_to_Z(comp->edges[0], Z),
                                         map_to_Z(comp->edges[1], Z));
        if (!strcmp(comp->operation, "-") && comp->arity == 2)
                return isl_map_sum(
                        map_to_Z(comp->edges[0], Z),
                        isl_map_neg(map_to_Z(comp->edges[1], Z)));
        if (!strcmp(comp->operation, "-") && comp->arity == 1)
                return isl_map_neg(map_to_Z(comp->edges[0], Z));
        if (!strcmp(comp->operation, "/") && comp->arity == 2)
                return affine_expression_div(comp, Z);
        if (!strcmp(comp->operation, "floord") && comp->arity == 2)
                return affine_expression_div(comp, Z);
        if (!strcmp(comp->operation, "ceild") && comp->arity == 2)
                return affine_expression_ceildiv(comp, Z);
        if (!strcmp(comp->operation, "max") && comp->arity == 2)
                return affine_expression_max(comp, Z);
        if (!strcmp(comp->operation, "*") && comp->arity == 2) {
                isl_map *rel0 = map_to_Z(comp->edges[0], Z);
                isl_map *rel1 = map_to_Z(comp->edges[1], Z);
                return map_mul(rel0, rel1);
        }
        return NULL;
}

static void add_computation_edge(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                pdg::function_call *call, computation *comp, int pos)
{
        char *endp;
        long int cst;
        cst = strtol(call->name->s.c_str(), &endp, 0);
        if (*call->name->s.c_str() && !*endp) {
                add_constant_edge(pdg, a2c, pos, cst, comp);
                return;
        }

        computation *child = new computation();
        child->domain = isl_set_copy(comp->domain);
        child->operation = strdup(call->name->s.c_str());
        child->arity = call->arguments.size();
        child->location = comp->location;

        edge *e = new edge;
        e->pos = pos;
        comp->edges.push_back(e);

        add_edges(dg, pdg, s2c, a2c, missing, call, child);

        isl_map *expr = affine_expression(child, a2c[NULL]);
        if (expr) {
                delete child;
                e->source = a2c[NULL];
                e->relation = expr;
        } else {
                dg->vertices.push_back(child);
                e->source = child;
                isl_space *dim = isl_set_get_space(comp->domain);
                e->relation = isl_map_identity(isl_space_map_from_set(dim));
        }
}

/* Add and edge from computation "at" to the integer computation "Z".
 * The edge maps the domain of "at" to the value of its dimension "index_dim".
 */
static void add_index_edge(dependence_graph *dg,
        computation *at, computation *Z, int pos, unsigned index_dim)
{
        unsigned total_out = isl_set_dim(at->domain, isl_dim_out);
        isl_space *dim = isl_set_get_space(at->domain);
        isl_map *id = isl_map_identity(isl_space_map_from_set(dim));
        id = isl_map_remove_dims(id, isl_dim_out, index_dim + 1,
                                     total_out - index_dim - 1);
        id = isl_map_remove_dims(id, isl_dim_out, 0, index_dim);
        add_split_edge(at, Z, pos, id, NULL);
}

/* Given a write access "a" inside the arguments of a function call
 * that was translated into computation "comp", make "a" point to
 * "comp" as its computation and fill in locations collected in
 * "missing" that should refer to the computation of "a".
 *
 * If the write access has an extension, i.e., if a whole row or
 * an even larger chunk of an array is written, then we introduce
 * a new node that corresponds to selecting the individual elements.
 * Dependences that flow from the write are then turned into edges that
 * point to the new computation.
 * If the write access extends the dimension by s, i.e., if the
 * chunk written is of dimension s, then the new computation with
 * operation "#at" has 1 + s arguments, one for the chunk itself
 * and one for each index in the chunk.
 * There is exactly one edge for each arguments position, each
 * selecting the appropriate dimensions from the domain.
 * The first argument ensures that the values of the chunks
 * are the same in both programs, while the remaining arguments
 * ensure that the same element from this chunk is selected.
 */
static void set_write_access(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                pdg::access *a, computation *comp)
{
        if (a->extension) {
                isl_ctx *ctx = pdg->get_isl_ctx();
                isl_map *ext = a->extension->get_isl_map(ctx);
                unsigned dim = isl_map_dim(ext, isl_dim_in);
                unsigned e_dim = isl_map_dim(ext, isl_dim_out);

                computation *at = new computation();
                at->domain = isl_set_apply(isl_set_copy(comp->domain),
                                                isl_map_copy(ext));
                at->operation = strdup("#at");
                at->arity = 1 + (e_dim - dim);
                at->location = comp->location;
                dg->vertices.push_back(at);

                ext = isl_map_reverse(ext);
                add_split_edge(at, comp, 0, ext, NULL);
                for (int i = 0; i < e_dim - dim; ++i)
                        add_index_edge(dg, at, a2c[NULL], 1 + i, dim + i);
                comp = at;
        }
        s2c[a] = comp;
        for (int k = 0; k < missing[a].size(); ++k)
                *missing[a][k] = comp;
}

/* For a given call or access "coa" which appears as argument "pos" of
 * a function call, add one or more edges to computation comp.
 * In case the argument is an access, the edges corresponding to
 * the dependences are added in add_dep_edges.
 * If it is itself a function call, then a single edge is added
 * to a new child and add_edges is called recursively on the
 * nested function call.
 * Return 1 if any edges were added for this position and 0 if not.
 * In particular, a write access should not be considered as an
 * argument and is therefore skipped.
 */
static int add_edges_at_pos(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                int pos, pdg::call_or_access *coa, computation *comp)
{
        if (coa->type == pdg::call_or_access::t_access) {
                pdg::access *a = coa->access;
                if (a->type == pdg::access::write) {
                        set_write_access(dg, pdg, s2c, a2c, missing, a, comp);
                        comp->arity--;
                        return 0;
                }
                int nested_dim = a->map->input() -
                                        isl_set_dim(comp->domain, isl_dim_set);
                if (nested_dim > 0)
                        add_nested_access_edge(dg, pdg, s2c, a2c, missing,
                                                pos, coa->access, comp);
                else
                        add_access_edges(dg, pdg, s2c, a2c, missing,
                                         pos, coa->access, comp);
        } else {
                add_computation_edge(dg, pdg, s2c, a2c, missing,
                                        coa->call, comp, pos);
        }
        return 1;
}

/* For each argument of function call call, we add one or more edges to
 * computation comp.
 */
static void add_edges(dependence_graph *dg, pdg::PDG *pdg,
                std::map<pdg::access *, computation *> &s2c,
                std::map<pdg::array *, computation *> &a2c,
                std::map<pdg::access *, std::vector<computation **> > &missing,
                pdg::function_call *call, computation *comp)
{
        int pos = 0;
        for (int i = 0; i < call->arguments.size(); ++i) {
                pdg::call_or_access *coa = call->arguments[i];

                if (add_edges_at_pos(dg, pdg, s2c, a2c,
                                     missing, pos, coa, comp))
                        ++pos;
        }

        return;
}

/* Add an edge from the output computation (out) to a computation computing
 * (part of) an output array.  node is pdg node corresponding to this
 * computation.  It is assumed to have a single write access, corresponding
 * to the output array.  exposed is that part of the node domain that
 * is not the source of an output dependence, i.e., those for which
 * the computed array elements are not overwritten.
 * Let f(i) be the index expression of the write access and let
 * pos be the position of the output array in the sorted list of
 * output arrays, then an edge is added to out with relation
 *
 *      { [ f(i), 0, pos ] -> [i] | i \in exposed }
 *
 * where the total dimension of the domain is out_dim, and this domain
 * is added to the domain of the output computation.
 */
static void add_output_edge(isl_ctx *ctx, computation *out, pdg::node *node,
        std::map<pdg::array *, int> &array2pos,
        computation *last,
        isl_set *exposed, unsigned nparam, unsigned out_dim)
{
        pdg::statement *s = node->statement;
        int j;
        for (j = 0;  j < s->accesses.size(); ++j)
                if (s->accesses[j]->type == pdg::access::write &&
                    s->accesses[j]->array->type == pdg::array::output)
                        break;
        pdg::array *array = s->accesses[j]->array;
        int pos = array2pos[array];
        isl_map *rel;
        if (!s->accesses[j]->extended_map)
                rel = isl_map_reverse(s->accesses[j]->map->get_isl_map(ctx));
        else
                rel = isl_map_reverse(s->accesses[j]->extended_map->get_isl_map(ctx));
        rel = isl_map_add_dims(rel, isl_dim_in,
                                out_dim - isl_map_dim(rel, isl_dim_in));
        rel = isl_map_intersect_range(rel, exposed);
        rel = isl_map_fix_input_si(rel, out_dim-1, pos);
        for (int k = array->dims.size(); k < out_dim-1; ++k)
                rel = isl_map_fix_input_si(rel, k, 0);

        out->domain = isl_set_union_disjoint(out->domain,
                                        isl_map_domain(isl_map_copy(rel)));
        add_split_edge(out, last, 0, rel, NULL);
}

static void connect_top_level_outputs(computation *comp,
        std::map<pdg::access *, computation *> &s2c,
        std::map<pdg::access *, std::vector<computation **> > &missing,
        pdg::node *node)
{
        pdg::statement *s = node->statement;

        for (int j = 0; j < s->top_outputs.size(); ++j) {
                pdg::access *a = s->top_outputs[j];
                s2c[a] = comp;
                for (int k = 0; k < missing[a].size(); ++k)
                        *missing[a][k] = comp;
        }
}

/* Add some computations to model a "#test" operation generated
 * by pers for data dependent assignments.
 * The test operation has three arguments.
 * The first is an "iteration" access with nested accesses corresponding
 * to the data dependent constructs in the conditions of the assignment.
 * The access relation maps iteration in the extended domain that
 * satisfy the conditions to 1, and those that do not satisfy the
 * conditions to 0.
 * The other two arguments are the operations or accesses performed
 * when the conditions hold (second argument) or do not hold (third
 * argument).
 *
 * We first create a copy computation and put it inside a #NEST
 * computation that add the required dimensions.
 * Then we split the extended iteration domain according to
 * the data dependent conditions using two child copy computations
 * with the appropriate dependence relations.
 * The second and third arguments of the "#test" operation are
 * then attached to these two child copy computations.
 */
static void add_test_computation(dependence_graph *dg, pdg::PDG *pdg,
        std::map<pdg::access *, computation *> &s2c,
        std::map<pdg::array *, computation *> &a2c,
        std::map<pdg::access *, std::vector<computation **> > &missing,
        pdg::node *node)
{
        isl_ctx *ctx = pdg->get_isl_ctx();
        pdg::statement *s = node->statement;

        computation *copy;
        copy = add_copy_computation(dg, node->source->get_isl_set(ctx),
                                    s->line);

        assert(s->top_function->arguments.size() == 3);
        pdg::call_or_access *coa = s->top_function->arguments[0];
        assert(coa->type == pdg::call_or_access::t_access);
        pdg::access *a = coa->access;

        computation *nest;
        nest = add_nest_computation(dg, pdg, s2c, a2c, missing, a, copy);

        isl_map *map_true = a->map->get_isl_map(ctx);
        isl_map *map_false = isl_map_copy(map_true);

        map_true = isl_map_fix_si(map_true, isl_dim_out, 0, 1);
        map_false = isl_map_fix_si(map_false, isl_dim_out, 0, 0);

        isl_space *dim = isl_space_domain(isl_map_get_space(map_true));
        dim = isl_space_map_from_set(dim);
        map_true = isl_map_intersect_domain(isl_map_identity(isl_space_copy(dim)),
                                            isl_map_domain(map_true));
        map_false = isl_map_intersect_domain(isl_map_identity(dim),
                                            isl_map_domain(map_false));

        computation *copy_true;
        copy_true = add_copy_computation(dg, isl_set_copy(copy->domain),
                                         copy->location);
        computation *copy_false;
        copy_false = add_copy_computation(dg, isl_set_copy(copy->domain),
                                         copy->location);

        add_split_edge(copy, copy_true, 0, map_true, NULL);
        add_split_edge(copy, copy_false, 0, map_false, NULL);

        add_edges_at_pos(dg, pdg, s2c, a2c, missing, 0,
                                s->top_function->arguments[1], copy_true);
        add_edges_at_pos(dg, pdg, s2c, a2c, missing, 0,
                                s->top_function->arguments[2], copy_false);

        connect_top_level_outputs(nest, s2c, missing, node);
}

static void add_top_level_computation(dependence_graph *dg, pdg::PDG *pdg,
        std::map<pdg::access *, computation *> &s2c,
        std::map<pdg::array *, computation *> &a2c,
        std::map<pdg::access *, std::vector<computation **> > &missing,
        pdg::node *node, __isl_keep isl_set *context)
{
        isl_ctx *ctx = pdg->get_isl_ctx();
        pdg::statement *s = node->statement;

        if (!strcmp(s->top_function->name->s.c_str(), "#test")) {
                add_test_computation(dg, pdg, s2c, a2c, missing, node);
                return;
        }

        computation *comp = new computation();
        comp->domain = node->source->get_isl_set(ctx);
        comp->domain = isl_set_intersect_params(comp->domain,
                                                isl_set_copy(context));
        comp->location = s->line;

        char *endp;
        long int cst;
        cst = strtol(s->top_function->name->s.c_str(), &endp, 0);
        if (*s->top_function->name->s.c_str() && !*endp) {
                comp->arity = 1;
                comp->operation = strdup("");
                add_constant_edge(pdg, a2c, 0, cst, comp);
        } else {
                comp->arity = s->top_function->arguments.size();
                comp->operation = strdup(s->top_function->name->s.c_str());
                add_edges(dg, pdg, s2c, a2c, missing, s->top_function, comp);
                isl_map *expr = affine_expression(comp, a2c[NULL]);
                if (expr) {
                        delete comp;
                        comp = new computation();
                        comp->domain = node->source->get_isl_set(ctx);
                        comp->domain = isl_set_intersect_params(comp->domain,
                                                        isl_set_copy(context));
                        comp->location = s->line;
                        comp->arity = 1;
                        comp->operation = strdup("");

                        edge *e = new edge;
                        e->pos = 0;
                        e->source = a2c[NULL];
                        e->relation = expr;
                        comp->edges.push_back(e);
                }
        }

        connect_top_level_outputs(comp, s2c, missing, node);

        dg->vertices.push_back(comp);
}

static void ensure_output_deps(pdg::PDG *pdg, pdg::array *array)
{
        for (int i = 0; i < array->analysis_performed.size(); ++i)
                if (array->analysis_performed[i]->type ==
                                                    pdg::dependence::output)
                        return;
        fprintf(stderr, "no output dependence analysis performed on array %s\n",
                array->name->s.c_str());
        exit(-1);
}

struct less_name {
        bool operator()(pdg::array* x, pdg::array* y) {
                return strcmp(x->name->s.c_str(), y->name->s.c_str()) < 0;
        }
};

/* Convert a pdg to a dependence_graph containing all the information
 * needed for equivalence checking.
 *
 * We first compute the total dimension of the output computation
 * and create input computations for each input array.
 * Then we create computations for each node in the pdg and
 * a special output computation.
 * For each node, we subsequently compute all the edges emanating
 * from the corresponding computation.
 * Finally, for each node writing to an output array, we compute
 * that part of the domain that computes part of the final
 * contents of the output array and add corresponding edges
 * to the output computation.
 */
dependence_graph *pdg_to_dg(pdg::PDG *pdg, unsigned out_dim,
        __isl_take isl_set *context)
{
        isl_space *dims;
        isl_ctx *ctx = pdg->get_isl_ctx();
        unsigned nparam = pdg->params.size();
        dependence_graph *dg = new dependence_graph;
        std::map<pdg::access *, computation *> s2c;
        std::map<pdg::array *, computation *> a2c;
        std::map<pdg::access *, std::vector<computation **> > missing;

        std::vector<pdg::array *> output_arrays;
        std::map<pdg::array *, int> array2pos;
        std::map<pdg::access *, isl_set *> exposed;

        if (context)
                context = isl_set_intersect(context, pdg->get_context_isl_set());
        else
                context = pdg->get_context_isl_set();

        /* add a special "input computation"; see add_iterator_edges */
        computation *comp;
        dims = isl_space_set_alloc(ctx, nparam, 1);
        isl_dim_set_parameter_names(dims, pdg->params);
        comp = new input_array("", dims);
        a2c[NULL] = comp;
        dg->vertices.push_back(comp);

        for (int i = 0; i < pdg->arrays.size(); ++i) {
                pdg::array *array = pdg->arrays[i];
                if (array->type == pdg::array::output) {
                        output_arrays.push_back(array);
                        ensure_output_deps(pdg, array);
                }

                if (array->type != pdg::array::input)
                        continue;

                computation *comp;
                dims = isl_space_set_alloc(ctx, nparam, array->dims.size());
                isl_dim_set_parameter_names(dims, pdg->params);
                comp = new input_array(array->name->s.c_str(), dims);
                a2c[array] = comp;
                dg->vertices.push_back(comp);
                dg->input_computations.push_back(comp);
        }
        if (!out_dim) {
                fprintf(stderr, "no output arrays\n");
                exit(-1);
        }
        std::sort(output_arrays.begin(), output_arrays.end(), less_name());
        for (int i = 0; i < output_arrays.size(); ++i) {
                dg->output_arrays.push_back(strdup(output_arrays[i]->name->s.c_str()));
                dg->output_array_dims.push_back(output_arrays[i]->dims.size());
                array2pos[output_arrays[i]] = i;
        }
        dg->out = new computation();
        dims = isl_space_set_alloc(ctx, nparam, out_dim);
        isl_dim_set_parameter_names(dims, pdg->params);
        dg->out->domain = isl_set_empty(dims);
        dg->out->operation = strdup("#Out");
        dg->out->arity = 1;
        dg->out->location = -1;
        for (int i = 0; i < pdg->nodes.size(); ++i) {
                pdg::node *node = pdg->nodes[i];
                pdg::statement *s = node->statement;

                add_top_level_computation(dg, pdg, s2c, a2c, missing, node,
                                            context);

                for (int j = 0; j < s->accesses.size(); ++j) {
                        if (s->accesses[j]->type != pdg::access::write)
                                continue;
                        pdg::access *a = s->accesses[j];

                        if (s->accesses[j]->array->type != pdg::array::output)
                                continue;

                        exposed[a] = node->source->get_isl_set(ctx);
                        if (a->extension) {
                                isl_map *ext = a->extension->get_isl_map(ctx);
                                exposed[a] = isl_set_apply(exposed[a], ext);
                        }
                        exposed[a] = isl_set_intersect_params(exposed[a],
                                                        isl_set_copy(context));
                }
        }
        for (int i = 0; i < pdg->dependences.size(); ++i) {
                pdg::dependence *dep = pdg->dependences[i];
                if (dep->type != pdg::dependence::output)
                        continue;
                pdg::access *a = dep->from_access;
                if (!exposed[a] || isl_set_is_empty(exposed[a]))
                        continue;
                isl_map *rel;
                if (!dep->extended_relation)
                        rel = dep->relation->get_isl_map(ctx);
                else
                        rel = dep->extended_relation->get_isl_map(ctx);
                exposed[a] = isl_set_subtract(exposed[a], isl_map_domain(rel));
        }
        for (int i = 0; i < pdg->nodes.size(); ++i) {
                pdg::node *node = pdg->nodes[i];
                pdg::statement *s = node->statement;
                for (int j = 0;  j < s->accesses.size(); ++j) {
                        pdg::access *a = s->accesses[j];
                        if (!exposed[a] || isl_set_is_empty(exposed[a])) {
                                isl_set_free(exposed[a]);
                                continue;
                        }
                        add_output_edge(ctx, dg->out, node, array2pos,
                                        s2c[a], exposed[a], nparam, out_dim);
                }
        }
        dg->context = context;
        return dg;
}