update isl for change in lexicographic optimization

[ppn.git] / eqv.cc

#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <algorithm>
#include <map>
#include <set>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>

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

#include <isl/ctx.h>
#include <isl/space.h>
#include <isl/set.h>
#include <isl/map.h>
#include <isl/constraint.h>
#include <isl/union_set.h>
#include <isl/printer.h>
#include "eqv_options.h"
#include "dependence_graph.h"

struct ops {
        std::vector<const char *> associative;
        std::vector<const char *> commutative;
};
int ops_init(void *user)
{
        struct ops **ops = (struct ops **)user;
        *ops = new struct ops;
        return 0;
}
void ops_clear(void *user)
{
        struct ops **ops = (struct ops **)user;
        delete *ops;
}

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

static unsigned update_out_dim(pdg::PDG *pdg, unsigned out_dim)
{
        for (int i = 0; i < pdg->arrays.size(); ++i) {
                pdg::array *array = pdg->arrays[i];
                if (array->type != pdg::array::output)
                        continue;
                if (array->dims.size() + 1 > out_dim)
                        out_dim = array->dims.size() + 1;
        }

        return out_dim;
}

static bool is_associative(const std::vector<const char *> &associative,
        const char *op)
{
        std::vector<const char *>::const_iterator iter;

        for (iter = associative.begin(); iter != associative.end(); ++iter)
                if (!strcmp(*iter, op))
                        return 1;
        return 0;
}

/* Return a computation that has an associative operation
 * that takes a different computation with the same operation
 * as one of its arguments.  Also return the edge from the
 * first to the second computation in *e.
 *
 * If no such computation exists, then return NULL.
 */
computation *dependence_graph::associative_node(edge **e,
        const std::vector<const char *> &associative)
{
        for (int i = 0; i < vertices.size(); ++i) {
                computation *comp = vertices[i];
                if (!is_associative(associative, comp->operation))
                        continue;
                for (int j = 0; j < comp->edges.size(); ++j) {
                        computation *other = comp->edges[j]->source;
                        if (comp->has_same_source(other))
                                continue;
                        if (strcmp(comp->operation, other->operation))
                                continue;
                        *e = comp->edges[j];
                        return comp;
                }
        }
        return NULL;
}

/* Splice the source of edge edge into the position of edge e,
 * removing all other edges with the same positions and bumping
 * up edges with a greater position.
 */
void dependence_graph::splice(computation *comp, edge *e)
{
        computation *source = e->source;
        int pos = e->pos;
        isl_map *map = isl_map_copy(e->relation);
        std::vector<struct edge *> new_edges;

        for (int i = 0; i < comp->edges.size(); ++i) {
                if (comp->edges[i]->pos == pos) {
                        delete comp->edges[i];
                        continue;
                }
                if (comp->edges[i]->pos > pos)
                        comp->edges[i]->pos += source->arity - 1;
                new_edges.push_back(comp->edges[i]);
        }
        for (int i = 0; i < source->edges.size(); ++i) {
                edge *old_e = source->edges[i];
                edge *e = new edge;
                e->source = old_e->source;
                e->pos = old_e->pos + pos;
                e->relation = isl_map_apply_range(
                                        isl_map_copy(map),
                                        isl_map_copy(old_e->relation));
                new_edges.push_back(e);
        }
        isl_map_free(map);
        comp->edges = new_edges;
        comp->arity += source->arity - 1;
}

/* Split computation comp into a part that always takes the edge e
 * and one that never takes that edge.
 * The second is returned and the initial computation is modified
 * to match the first.
 *
 * We need to be careful not to destroy e, as it is still used
 * by the calling method.
 */
computation *dependence_graph::split_comp(computation *comp, edge *e)
{
        computation *dup = new computation;
        dup->original = comp->original ? comp->original : comp;
        dup->operation = strdup(comp->operation);
        dup->arity = comp->arity;
        dup->location = comp->location;

        isl_map *map = isl_map_copy(e->relation);
        map = isl_map_reverse(map);
        isl_set *dom = isl_set_apply(isl_set_copy(e->source->domain), map);
        dup->domain = isl_set_subtract(isl_set_copy(comp->domain),
                                        isl_set_copy(dom));
        comp->domain = isl_set_intersect(comp->domain, dom);

        std::vector<struct edge *> old_edges = comp->edges;
        comp->edges.clear();

        for (int i = 0; i < old_edges.size(); ++i) {
                if (old_edges[i] == e) {
                        comp->edges.push_back(e);
                        continue;
                }

                edge *e = old_edges[i];
                isl_map *map, *map_dup;
                map = isl_map_copy(e->relation);
                map_dup = isl_map_copy(map);

                map = isl_map_intersect_domain(map, isl_set_copy(comp->domain));
                map_dup = isl_map_intersect_domain(map_dup,
                                                    isl_set_copy(dup->domain));
                add_split_edge(comp, e->source, e->pos, map, NULL);
                add_split_edge(dup, e->source, e->pos, map_dup, NULL);
                delete e;
        }

        vertices.push_back(dup);
        return dup;
}

/* If any edge from comp points to comp_orig, then split it
 * into two edges, one still pointing to comp_orig and the
 * other pointing to comp_dup.
 * comp_orig and comp_dup are assumed to have disjoint domains
 * and the edge relations are adjusted according to these domains.
 */
void dependence_graph::split_edges(computation *comp,
                 computation *comp_orig, computation *comp_dup)
{
        std::vector<struct edge *> old_edges = comp->edges;
        comp->edges.clear();

        for (int i = 0; i < old_edges.size(); ++i) {
                edge *e = old_edges[i];

                if (e->source != comp_orig) {
                        comp->edges.push_back(e);
                        continue;
                }

                isl_map *map_orig, *map_dup;
                map_orig = isl_map_copy(e->relation);
                map_dup = isl_map_copy(map_orig);
                map_orig = isl_map_intersect_range(map_orig,
                                            isl_set_copy(comp_orig->domain));
                map_dup = isl_map_intersect_range(map_dup,
                                            isl_set_copy(comp_dup->domain));
                add_split_edge(comp, comp_orig, e->pos, map_orig, NULL);
                add_split_edge(comp, comp_dup, e->pos, map_dup, NULL);
                delete e;
        }
}

void dependence_graph::split_edges(computation *comp_orig,
                                        computation *comp_dup)
{
        split_edges(out, comp_orig, comp_dup);
        for (int i = 0; i < vertices.size(); ++i)
                split_edges(vertices[i], comp_orig, comp_dup);
}

/* Replace all nested calls of an associative operator,
 * by a call with the nested call spliced into the first call.
 *
 * For each nested call we find, we first check if there are
 * any other edges with the same position.  If not, then
 * the nested call is performed for each iteration of the computation
 * and we can simplify splice the nested call.

 * Otherwise, we first create a duplicate computation for the iterations
 * that do not take the found edge and adjust the edges of both computations
 * to their domains.  This meand that the edge corresponding to the
 * nested call will no longer appear in the duplicated computation
 * and the other edges with the same position will no longer appear
 * in the original computation.
 * Then we splice the nested call in the original computation.
 * Finally, we split all edges that pointed to the original computation
 * into two edges, one going to the original computation and one
 * going to the duplicated computation.
 */
void dependence_graph::flatten_associative_operators(
        const std::vector<const char *> &associative)
{
        edge *e;
        computation *comp;

        while ((comp = associative_node(&e, associative)) != NULL) {
                computation *comp_dup = NULL;
                int j;
                for (j = 0; j < comp->edges.size(); ++j) {
                        if (comp->edges[j] == e)
                                continue;
                        if (comp->edges[j]->pos == e->pos)
                                break;
                }
                if (j != comp->edges.size())
                        comp_dup = split_comp(comp, e);
                splice(comp, e);
                if (comp_dup)
                        split_edges(comp, comp_dup);
        }
}

struct eq_node {
private:
        isl_map                 *got;
        isl_map                 *lost;
public:
        computation             *comp[2];
        isl_map                 *want;
        isl_map                 *need;
        std::set<eq_node *>      assumed;
        unsigned                 closed : 1;
        unsigned                 narrowing : 1;
        unsigned                 widening : 1;
        unsigned                 invalidated : 1;
        unsigned                 reset : 1;
        std::string              trace;
        isl_map                 *lost_sample;

        eq_node(computation *c1, computation *c2,
                        isl_map *w) : want(w), closed(0),
                        widening(0), narrowing(0), invalidated(0), reset(0),
                        need(NULL), got(NULL), lost(NULL), lost_sample(NULL) {
                comp[0] = c1;
                comp[1] = c2;
        }
        bool is_still_valid();
        void collect_open_assumed(std::set<eq_node *> &c);
        ~eq_node() {
                isl_map_free(want);
                isl_map_free(need);
                isl_map_free(got);
                isl_map_free(lost);
                isl_map_free(lost_sample);
        }
        void compute_got() {
                assert(lost);
                got = isl_map_copy(want);
                got = isl_map_subtract(got, isl_map_copy(lost));
        }
        isl_map *get_got() {
                if (!got)
                        compute_got();
                return isl_map_copy(got);
        }
        isl_map *peek_got() {
                if (!got)
                        compute_got();
                return got;
        }
        void compute_lost() {
                assert(got);
                lost = isl_map_copy(want);
                lost = isl_map_subtract(lost, isl_map_copy(got));
        }
        isl_map *get_lost() {
                if (!lost)
                        compute_lost();
                return isl_map_copy(lost);
        }
        isl_map *peek_lost() {
                if (!lost)
                        compute_lost();
                return lost;
        }
        void set_got(isl_map *got) {
                isl_map_free(this->lost);
                isl_map_free(this->got);
                this->lost = NULL;
                this->got = got;
        }
        void set_lost(isl_map *lost) {
                isl_map_free(this->lost);
                isl_map_free(this->got);
                this->got = NULL;
                this->lost = lost;
        }
        std::string to_string();
};

std::string eq_node::to_string()
{
        std::ostringstream strm;
        strm << comp[0]->location << "," << comp[0]->operation
             << "/" << comp[0]->arity;
        strm << " <-> ";
        strm << comp[1]->location << "," << comp[1]->operation
             << "/" << comp[1]->arity;
        strm << std::endl;
        return strm.str();
}

bool eq_node::is_still_valid()
{
        if (invalidated)
                return 0;

        std::set<eq_node *>::iterator i;
        for (i = assumed.begin(); i != assumed.end(); ++i) {
                assert(*i != this);
                if ((*i)->reset ||
                    ((*i)->closed && !(*i)->is_still_valid())) {
                        invalidated = 1;
                        return 0;
                }
        }
        return 1;
}

void eq_node::collect_open_assumed(std::set<eq_node *> &c)
{
        std::set<eq_node *>::iterator i;
        for (i = assumed.begin(); i != assumed.end(); ++i) {
                if ((*i)->closed)
                        (*i)->collect_open_assumed(c);
                else
                        c.insert(*i);
        }
}

/* A comp_pair contains all the edges that have the same pair
 * of computations.
 */
struct comp_pair {
        computation             *comp[2];

        std::vector<eq_node *>   nodes;

        eq_node *tabled(eq_node *node);
        eq_node *last_ancestor(eq_node *node);
        ~comp_pair();
};

comp_pair::~comp_pair()
{
        std::vector<eq_node *>::iterator i;
        for (i = nodes.begin(); i != nodes.end(); ++i)
                delete *i;
}

eq_node *comp_pair::tabled(eq_node *node)
{
        std::vector<eq_node *>::iterator i;

        for (i = nodes.begin(); i != nodes.end(); ++i) {
                if (*i == node)
                        continue;
                if (!(*i)->closed)
                        continue;
                if (!(*i)->is_still_valid())
                        continue;
                int is_subset;
                is_subset = isl_map_is_subset(node->want, (*i)->want);
                assert(is_subset >= 0);
                if (!is_subset)
                        continue;
                return *i;
        }
        return NULL;
}

eq_node *comp_pair::last_ancestor(eq_node *node)
{
        std::vector<eq_node *>::reverse_iterator i;
        for (i = nodes.rbegin(); i != nodes.rend(); ++i) {
                if (*i == node)
                        continue;
                if ((*i)->closed)
                        continue;
                return *i;
        }
        return NULL;
}

typedef std::pair<computation *, computation *> computation_pair;
typedef std::map<computation_pair, comp_pair *> c2p_t;

struct equivalence_checker {
        c2p_t c2p;
        struct options *options;
        int widenings;
        int narrowings;

        equivalence_checker(struct options *o) :
                options(o), widenings(0), narrowings(0) {}
        comp_pair *get_comp_pair(eq_node *node);
        void handle(eq_node *node);
        void dismiss(eq_node *node);
        void handle_propagation(eq_node *node);
        void handle_copy(eq_node *node, int i);
        void handle_widening(eq_node *node, eq_node *a);
        void handle_with_ancestor(eq_node *node, eq_node *a);
        isl_map *lost_at_pos(eq_node *node, int pos1, int pos2);
        void handle_propagation_same(eq_node *node);
        void handle_propagation_comm(eq_node *node);
        void handle_narrowing(eq_node *node);

        isl_map *lost_from_propagation(eq_node *parent,
                eq_node *node, edge *e1, edge *e2);

        void init_trace(eq_node *node);
        void extend_trace_propagation(eq_node *node, eq_node *child,
                                        edge *e1, edge *e2);
        void extend_trace_widening(eq_node *node, eq_node *child);

        bool is_commutative(const char *op);

        ~equivalence_checker();
};

equivalence_checker::~equivalence_checker()
{
        c2p_t::iterator i;

        for (i = c2p.begin(); i != c2p.end(); ++i)
                delete (*i).second;
}

bool equivalence_checker::is_commutative(const char *op)
{
        std::vector<const char *>::iterator iter;

        for (iter = options->ops->commutative.begin();
             iter != options->ops->commutative.end(); ++iter)
                if (!strcmp(*iter, op))
                        return 1;
        return 0;
}

comp_pair *equivalence_checker::get_comp_pair(eq_node *node)
{
        c2p_t::iterator i;
        comp_pair *cp;

        i = c2p.find(computation_pair(node->comp[0], node->comp[1]));
        if (i == c2p.end()) {
                cp = new comp_pair;
                c2p[computation_pair(node->comp[0], node->comp[1])] = cp;
        } else
                cp = (*i).second;
        return cp;
}

void equivalence_checker::dismiss(eq_node *node)
{
        /* unclosed nodes weren't added to c2p->nodes */
        if (node && !node->closed)
                delete node;
}

void equivalence_checker::init_trace(eq_node *node)
{
        isl_ctx *ctx;
        isl_printer *prn;

        if (!options->trace_error)
                return;
        if (isl_map_is_empty(node->peek_lost()))
                return;
        node->trace = node->to_string();
        std::cerr << node->trace;
        isl_map_free(node->lost_sample);
        node->lost_sample = isl_map_from_basic_map(
                                isl_map_sample(node->get_lost()));
        ctx = isl_map_get_ctx(node->lost_sample);
        prn = isl_printer_to_file(ctx, stderr);
        prn = isl_printer_print_map(prn, node->lost_sample);
        prn = isl_printer_end_line(prn);
        isl_printer_free(prn);
}

void equivalence_checker::extend_trace_propagation(eq_node *node,
        eq_node *child, edge *e1, edge *e2)
{
        isl_ctx *ctx;
        isl_printer *prn;

        if (!options->trace_error)
                return;
        if (!child->lost_sample)
                return;
        if (node->lost_sample)
                return;
        node->trace = child->trace + node->to_string();
        std::cerr << node->trace;
        node->lost_sample = isl_map_copy(child->lost_sample);
        if (e1)
                node->lost_sample = isl_map_apply_domain(node->lost_sample,
                                        isl_map_reverse(
                                            isl_map_copy(e1->relation)));
        if (e2)
                node->lost_sample = isl_map_apply_range(node->lost_sample,
                                        isl_map_reverse(
                                            isl_map_copy(e2->relation)));
        node->lost_sample = isl_map_intersect(node->lost_sample,
                                                isl_map_copy(node->want));
        node->lost_sample = isl_map_from_basic_map(
                                        isl_map_sample(node->lost_sample));
        ctx = isl_map_get_ctx(node->lost_sample);
        prn = isl_printer_to_file(ctx, stderr);
        prn = isl_printer_print_map(prn, node->lost_sample);
        prn = isl_printer_end_line(prn);
        isl_printer_free(prn);
}

void equivalence_checker::extend_trace_widening(eq_node *node, eq_node *child)
{
        isl_ctx *ctx;
        isl_printer *prn;

        if (!options->trace_error)
                return;
        if (!child->lost_sample)
                return;
        if (isl_map_is_empty(node->peek_lost()))
                return;
        node->trace = child->trace + node->to_string();
        std::cerr << node->trace;
        node->lost_sample = isl_map_from_basic_map(
                                        isl_map_sample(node->get_lost()));
        ctx = isl_map_get_ctx(node->lost_sample);
        prn = isl_printer_to_file(ctx, stderr);
        prn = isl_printer_print_map(prn, node->lost_sample);
        prn = isl_printer_end_line(prn);
        isl_printer_free(prn);
}

/* Check for which subset (got) of the want relation equivelance
 * holds for the pair of computations in the equivalence node.
 *
 * We first handle the easy cases: empty want, input computations
 * and tabled nodes.
 *
 * If the current node is not a narrowing or a widening node
 * and we can find an ancestor with the same pair of computations,
 * then we will try to apply induction or widening in handle_with_ancestor.
 * However, if the requested relation (want) of the ancestor is a strict
 * subset of that of the current node, then we have already applied
 * widening on an intermediate node (with a different pair of computations)
 * so we shouldn't apply widening again (and we can't apply induction
 * because the relation of the ancestor is a strict subset).
 *
 * In all other cases we try to apply propagation in handle_propagation.
 */
void equivalence_checker::handle(eq_node *node)
{
        computation *comp1 = node->comp[0];
        computation *comp2 = node->comp[1];

        if (!comp1->is_copy() && !comp2->is_copy() &&
            (strcmp(comp1->operation, comp2->operation) ||
             comp1->arity != comp2->arity)) {
                node->set_lost(isl_map_copy(node->want));
                init_trace(node);
                return;
        }

        eq_node *s = NULL;
        eq_node *a = NULL;

        isl_map *want = node->want;
        node->need = isl_map_empty(isl_map_get_space(want));
        int empty = isl_map_is_empty(want);
        assert(empty >= 0);
        if (empty) {
                node->set_lost(isl_map_empty(isl_map_get_space(want)));
                return;
        }
        if (node->comp[0]->is_input() && node->comp[1]->is_input()) {
                if (strcmp(node->comp[0]->operation, node->comp[1]->operation)) {
                        node->set_lost(isl_map_copy(want));
                        return;
                }
                node->set_lost(isl_map_subtract(
                        isl_map_copy(node->want),
                        isl_map_identity(isl_map_get_space(want))));
                return;
        }

        comp_pair *cp = get_comp_pair(node);
        if ((s = cp->tabled(node)) != NULL) {
                node->set_lost(isl_map_intersect(s->get_lost(),
                                isl_map_copy(node->want)));
                s->collect_open_assumed(node->assumed);
                return;
        }

        cp->nodes.push_back(node);

        if (!node->narrowing && !node->widening &&
                    (a = cp->last_ancestor(node)) != NULL) {
                int is_subset;
                is_subset = isl_map_is_strict_subset(a->want, node->want);
                assert(is_subset >= 0);
                if (is_subset)
                        handle_propagation(node);
                else
                        handle_with_ancestor(node, a);
        } else
                handle_propagation(node);
        node->closed = 1;
}

/* Check if we can apply propagation to prove equivalence of the given node.
 * First check if we need to apply copy propagation and if not
 * check if the operations are the same and apply "regular" propagation.
 *
 * After we get back, we need to check that any induction hypotheses
 * we have used in the process of proving the node hold.
 * If not, we replace the obtained relation (got) by that
 * of a narrowing node in handle_narrowin.
 */
void equivalence_checker::handle_propagation(eq_node *node)
{
        if (node->comp[0]->is_copy())
                handle_copy(node, 0);
        else if (node->comp[1]->is_copy())
                handle_copy(node, 1);
        else
                handle_propagation_same(node);
        node->assumed.erase(node);
        isl_map *lost = node->get_lost();
        lost = isl_map_intersect(lost, isl_map_copy(node->need));
        int is_empty = isl_map_is_empty(lost);
        isl_map_free(lost);
        assert(is_empty >= 0);
        if (!is_empty)
                handle_narrowing(node);
}

/* When applying the mappings on expansion edges to both sides of
 * an equivalence relation, each element in the original equivalence
 * relation is mapped to many elements on both sides.  For example,
 * if the original equivalence relation has i R i' for i = i', then the new
 * equivalence relation may have (i,j) R (i',j') for i = i' and
 * 0 <= j,j' <= 10, expressing that all elements of the row read by
 * iteration i should be equal to all elements of the row read by i'.
 * Instead, we want to express that each individual element read by i'
 * should be equal to the corresponding element read by i'.
 * In the example, we need to introduce the equality j = j'.
 *
 * We first check that the number of dimensions added by the expansions
 * is the same in both programs.  Then we construct an identity relation
 * between this number of dimensions, lift it to the space of the
 * new equivalence relation and take the intersection.
 *
 * We have to be careful, though, that we don't loose any array elements
 * by taking the intersection.  In particular, the expansion maps
 * a given read operation to iterations of an extended domain that
 * reads all the individual array elements.  The read operation is
 * only equivalent to some other read operation if all the reads
 * of the individual array elements are equivalent.
 * We therefore need to make sure that adding the equalities does
 * not remove any of the reads in the extended domain.
 * We do this by projecting out the extra dimensions on one side
 * from both "want" and "want \cap eq".  The resulting maps should
 * be the same.  We do this check for both sides separately.
 *
 * If either of the above tests fails, then we simply return
 * the original over-ambitious want.
 */
static isl_map *expansion_want(isl_map *want, edge *e1, edge *e2)
{
        unsigned n_in, n_out;
        unsigned s_dim_1 = isl_map_dim(e1->relation, isl_dim_out) -
                                isl_map_dim(e1->relation, isl_dim_in);
        unsigned s_dim_2 = isl_map_dim(e2->relation, isl_dim_out) -
                                isl_map_dim(e2->relation, isl_dim_in);
        if (s_dim_1 != s_dim_2)
                return want;

        isl_space *dim = isl_map_get_space(e1->relation);
        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_1);
        dim = isl_space_add_dims(dim, isl_dim_out, s_dim_1);
        isl_basic_map *id = isl_basic_map_identity(dim);

        dim = isl_space_range(isl_map_get_space(e1->relation));
        isl_basic_map *s_1 = isl_basic_map_identity(isl_space_map_from_set(dim));
        s_1 = isl_basic_map_remove_dims(s_1, isl_dim_in, 0,
                                isl_basic_map_dim(s_1, isl_dim_in) - s_dim_1);
        id = isl_basic_map_apply_domain(id, s_1);

        dim = isl_space_range(isl_map_get_space(e2->relation));
        isl_basic_map *s_2 = isl_basic_map_identity(isl_space_map_from_set(dim));
        s_2 = isl_basic_map_remove_dims(s_2, isl_dim_in, 0,
                                isl_basic_map_dim(s_2, isl_dim_in) - s_dim_2);
        id = isl_basic_map_apply_range(id, s_2);

        bool unmatched = false;
        isl_map *matched_want;
        isl_map *proj_want, *proj_matched;
        matched_want = isl_map_intersect(isl_map_copy(want),
                                         isl_map_from_basic_map(id));

        n_in = isl_map_dim(want, isl_dim_in);
        proj_want = isl_map_remove_dims(isl_map_copy(want),
                    isl_dim_in, n_in - s_dim_1, s_dim_1);
        proj_matched = isl_map_remove_dims(isl_map_copy(matched_want),
                    isl_dim_in, n_in - s_dim_1, s_dim_1);
        if (!isl_map_is_equal(proj_want, proj_matched))
                unmatched = true;
        isl_map_free(proj_want);
        isl_map_free(proj_matched);

        n_out = isl_map_dim(want, isl_dim_out);
        proj_want = isl_map_remove_dims(isl_map_copy(want),
                    isl_dim_out, n_out - s_dim_2, s_dim_2);
        proj_matched = isl_map_remove_dims(isl_map_copy(matched_want),
                    isl_dim_out, n_out - s_dim_2, s_dim_2);
        if (!isl_map_is_equal(proj_want, proj_matched))
                unmatched = true;
        isl_map_free(proj_want);
        isl_map_free(proj_matched);

        if (unmatched) {
                isl_map_free(matched_want);
                return matched_want;
        }
        isl_map_free(want);
        return matched_want;
}

static isl_map *propagation_want(eq_node *node, edge *e1, edge *e2)
{
        isl_map *new_want;

        new_want = isl_map_copy(node->want);
        if (e1)
                new_want = isl_map_apply_domain(new_want,
                            isl_map_copy(e1->relation));
        if (e2)
                new_want = isl_map_apply_range(new_want,
                            isl_map_copy(e2->relation));

        if (e1 && e2 && e1->type == edge::expansion)
                new_want = expansion_want(new_want, e1, e2);

        return isl_map_detect_equalities(new_want);
}

static eq_node *propagation_node(eq_node *node, edge *e1, edge *e2)
{
        isl_map *new_want;
        computation *comp1 = e1 ? e1->source : node->comp[0];
        computation *comp2 = e2 ? e2->source : node->comp[1];

        new_want = propagation_want(node, e1, e2);
        eq_node *child = new eq_node(comp1, comp2, new_want);
        return child;
}

static isl_map *set_to_empty(isl_map *map)
{
        isl_space *space;
        space = isl_map_get_space(map);
        isl_map_free(map);
        return isl_map_empty(space);
}

/* Propagate "lost" from child back to parent.
 */
isl_map *equivalence_checker::lost_from_propagation(eq_node *parent,
        eq_node *node, edge *e1, edge *e2)
{
        isl_map *new_lost;
        new_lost = node->get_lost();

        if (e1)
                new_lost = isl_map_apply_domain(new_lost,
                                    isl_map_reverse(isl_map_copy(e1->relation)));
        if (e2)
                new_lost = isl_map_apply_range(new_lost,
                                    isl_map_reverse(isl_map_copy(e2->relation)));

        new_lost = isl_map_intersect(new_lost, isl_map_copy(parent->want));

        extend_trace_propagation(parent, node, e1, e2);

        return new_lost;
}

void equivalence_checker::handle_copy(eq_node *node, int i)
{
        std::vector<edge *>::iterator iter;
        computation *copy = node->comp[i];
        isl_map *want = node->want;
        isl_map *lost;

        lost = isl_map_empty(isl_map_get_space(want));

        for (iter = copy->edges.begin(); iter != copy->edges.end(); ++iter) {
                edge *e = *iter;
                eq_node *child;
                if (i == 0)
                        child = propagation_node(node, e, NULL);
                else
                        child = propagation_node(node, NULL, e);

                if (!child)
                        continue;

                handle(child);

                isl_map *new_lost;
                if (i == 0)
                        new_lost = lost_from_propagation(node, child, e, NULL);
                else
                        new_lost = lost_from_propagation(node, child, NULL, e);
                dismiss(child);

                lost = isl_map_union_disjoint(lost, new_lost);
        }
        node->set_lost(lost);
}

/* Compute and return the part of "want" that is lost when propagating
 * over all edges of argument position pos1 in the first program
 * and argument position pos2 in the second program.
 * Since the domains of the dependence mappings on edges with the same
 * argument position partition the domain of the computation (and are
 * therefore disjoint), we simply need to take the disjoint union
 * of all losts over all pairs of edges.
 */
isl_map *equivalence_checker::lost_at_pos(eq_node *node, int pos1, int pos2)
{
        isl_map *lost;
        lost = isl_map_empty(isl_map_get_space(node->want));

        for (int i = 0; i < node->comp[0]->edges.size(); ++i) {
                edge *e1 = node->comp[0]->edges[i];
                if (e1->pos != pos1)
                        continue;
                for (int j = 0; j < node->comp[1]->edges.size(); ++j) {
                        edge *e2 = node->comp[1]->edges[j];
                        if (e2->pos != pos2)
                                continue;

                        eq_node *child;
                        isl_map *new_lost = NULL;

                        child = propagation_node(node, e1, e2);
                        handle(child);
                        new_lost = lost_from_propagation(node, child, e1, e2);

                        lost = isl_map_union_disjoint(lost, new_lost);

                        std::set<eq_node *>::iterator k;
                        for (k = child->assumed.begin();
                             k != child->assumed.end(); ++k)
                                node->assumed.insert(*k);
                        dismiss(child);
                }
        }
        return lost;
}

/* Compute the lost that results from propagation on a pair of
 * computations with a commutative operation.
 *
 * We first compute the losts for each pair of argument positions
 * and store the result in lost[pos1][pos2].
 * Then we perform a backtracking search over all permutations
 * of the arguments. For each permutation, we compute the lost
 * relation as the union of the losts over all arguments.
 * The final lost is the intersection of all these losts over all
 * permutations.
 */
void equivalence_checker::handle_propagation_comm(eq_node *node)
{
        int trace_error = options->trace_error;
        trace_error = 0;

        unsigned r = node->comp[0]->arity;

        std::vector<std::vector<isl_map *> > lost;

        for (int i = 0; i < r; ++i) {
                std::vector<isl_map *> row;
                for (int j = 0; j < r; ++j) {
                        isl_map *pos_lost;
                        pos_lost = lost_at_pos(node, i, j);
                        row.push_back(pos_lost);
                }
                lost.push_back(row);
        }

        int level;
        std::vector<int> perm;
        std::vector<isl_map *> lost_at;
        for (level = 0; level < r; ++level) {
                perm.push_back(0);
                lost_at.push_back(NULL);
        }

        isl_map *total_lost;
        total_lost = isl_map_copy(node->want);

        level = 0;
        while (level >= 0) {
                if (perm[level] == r) {
                        perm[level] = 0;
                        --level;
                        if (level >= 0)
                                ++perm[level];
                        continue;
                }
                int l;
                for (l = 0; l < level; ++l)
                        if (perm[l] == perm[level])
                                break;
                if (l != level) {
                        ++perm[level];
                        continue;
                }

                isl_map_free(lost_at[level]);
                lost_at[level] = isl_map_copy(lost[level][perm[level]]);
                if (level != 0)
                        lost_at[level] = isl_map_union(lost_at[level],
                                            isl_map_copy(lost_at[level - 1]));

                if (level < r - 1) {
                        ++level;
                        continue;
                }

                lost_at[level] = isl_map_coalesce(lost_at[level]);
                total_lost = isl_map_intersect(total_lost,
                                            isl_map_copy(lost_at[level]));
                ++perm[level];
        }

        for (int i = 0; i < r; ++i)
                isl_map_free(lost_at[i]);

        for (int i = 0; i < r; ++i)
                for (int j = 0; j < r; ++j)
                        isl_map_free(lost[i][j]);

        node->set_lost(total_lost);

        options->trace_error = trace_error;
        init_trace(node);
}

/* Compute the lost that results from propagation on a pair of
 * computations.
 *
 * First, for functions of zero arity (i.e., constants), equivalence
 * always holds and the lost relation is empty.
 * For commutative operations, the computation is delegated
 * to handle_propagation_comm.
 * Otherwise, we simply take the union of the losts over each
 * argument position (always taking the same argument position
 * in both programs).
 */
void equivalence_checker::handle_propagation_same(eq_node *node)
{
        unsigned r = node->comp[0]->arity;
        if (r == 0) {
                node->set_lost(isl_map_empty(isl_map_get_space(node->want)));
                return;
        }
        if (is_commutative(node->comp[0]->operation)) {
                handle_propagation_comm(node);
                return;
        }

        isl_map *lost;
        lost = isl_map_empty(isl_map_get_space(node->want));
        for (int i = 0; i < r; ++i) {
                isl_map *pos_lost;
                pos_lost = lost_at_pos(node, i, i);
                lost = isl_map_union(lost, pos_lost);
        }
        lost = isl_map_coalesce(lost);
        node->set_lost(lost);
}

void equivalence_checker::handle_widening(eq_node *node, eq_node *a)
{
        isl_map *wants;
        isl_map *aff;

        wants = isl_map_union(isl_map_copy(node->want), isl_map_copy(a->want));
        aff = isl_map_from_basic_map(isl_map_affine_hull(wants));
        aff = isl_map_intersect_domain(aff,
                                isl_set_copy(node->comp[0]->domain));
        aff = isl_map_intersect_range(aff,
                                isl_set_copy(node->comp[1]->domain));

        eq_node *child = new eq_node(node->comp[0], node->comp[1], aff);
        child->widening = 1;
        widenings++;
        handle(child);
        node->set_lost(isl_map_intersect(child->get_lost(),
                                         isl_map_copy(node->want)));
        extend_trace_widening(node, child);
        std::set<eq_node *>::iterator i;
        for (i = child->assumed.begin(); i != child->assumed.end(); ++i)
                node->assumed.insert(*i);
        dismiss(child);
}

/* Perform induction, if possible, and widening if we have to.
 */
void equivalence_checker::handle_with_ancestor(eq_node *node, eq_node *a)
{
        if (a->narrowing || isl_map_is_subset(node->want, a->want)) {
                isl_map *need;
                need = isl_map_intersect(isl_map_copy(node->want),
                                         isl_map_copy(a->want));
                node->set_lost(isl_map_subtract(isl_map_copy(node->want),
                                                isl_map_copy(a->want)));
                node->assumed.insert(a);
                a->need = isl_map_union(a->need, need);
        } else {
                handle_widening(node, a);
        }
}

struct narrowing_data {
        eq_node *node;
        equivalence_checker *ec;
        isl_map *new_got;
};

static isl_stat basic_handle_narrowing(__isl_take isl_basic_map *bmap,
        void *user)
{
        narrowing_data *data = (narrowing_data *)user;

        eq_node *child;
        child = new eq_node(data->node->comp[0], data->node->comp[1],
                            isl_map_from_basic_map(bmap));
        child->narrowing = 1;
        data->ec->narrowings++;
        data->ec->handle(child);
        data->new_got = isl_map_union(data->new_got, child->get_got());

        std::set<eq_node *>::iterator i;
        for (i = child->assumed.begin(); i != child->assumed.end(); ++i)
                data->node->assumed.insert(*i);
        data->ec->dismiss(child);

        return isl_stat_ok;
}

/* Construct and handle narrowing nodes for the given node.
 *
 * If the node itself was already a narrowing node, then we
 * simply return the empty relation.
 * Otherwise, we consider each basic relation in the obtained
 * relation, construct a new node with that basic relation as
 * requested relation and take the union of all obtained relations.
 */
void equivalence_checker::handle_narrowing(eq_node *node)
{
        node->reset = 1;
        node->assumed.clear();

        isl_map *want = node->want;
        isl_map *new_got = isl_map_empty(isl_map_get_space(want));
        if (!options->narrowing || node->narrowing) {
                node->set_got(new_got);
                return;
        }

        narrowing_data data = { node, this, new_got };
        isl_map_foreach_basic_map(node->peek_got(), &basic_handle_narrowing,
                                  &data);
        node->set_got(data.new_got);
}

static __isl_give isl_union_set *update_result(__isl_take isl_union_set *res,
        const char *array_name, isl_map *map, int first, int n)
{
        if (isl_map_is_empty(map)) {
                isl_map_free(map);
                return res;
        }

        isl_set *range = isl_map_range(map);
        range = isl_set_remove_dims(range, isl_dim_set, first, n);
        range = isl_set_coalesce(range);
        range = isl_set_set_tuple_name(range, array_name);
        res = isl_union_set_add_set(res, range);

        return res;
}

struct check_equivalence_data {
        dependence_graph *dg1;
        dependence_graph *dg2;
        equivalence_checker *ec;
        isl_map *got;
        isl_map *lost;
};

static isl_stat basic_check(__isl_take isl_basic_map *bmap, void *user)
{
        check_equivalence_data *data = (check_equivalence_data *)user;

        eq_node *root = new eq_node(data->dg1->out, data->dg2->out,
                                    isl_map_from_basic_map(bmap));
        data->ec->handle(root);
        data->got = isl_map_union_disjoint(data->got, root->get_got());
        data->lost = isl_map_union_disjoint(data->lost, root->get_lost());
        data->ec->dismiss(root);

        return isl_stat_ok;
}

static int check_equivalence_array(isl_ctx *ctx, equivalence_checker *ec,
        dependence_graph *dg1, dependence_graph *dg2, int array1, int array2,
        isl_union_set **proved, isl_union_set **not_proved)
{
        unsigned n_in;
        const char *array_name = dg1->output_arrays[array1];
        isl_set *out1 = isl_set_copy(dg1->out->domain);
        isl_set *out2 = isl_set_copy(dg2->out->domain);
        unsigned dim1 = isl_set_n_dim(out1);
        unsigned dim2 = isl_set_n_dim(out2);
        unsigned array_dim = dg1->output_array_dims[array1];
        out1 = isl_set_fix_dim_si(out1, dim1-1, array1);
        out2 = isl_set_fix_dim_si(out2, dim2-1, array2);
        out1 = isl_set_remove_dims(out1, isl_dim_set, dim1 - 1, 1);
        out2 = isl_set_remove_dims(out2, isl_dim_set, dim2 - 1, 1);
        int equal = isl_set_is_equal(out1, out2);
        isl_set_free(out1);
        isl_set_free(out2);
        assert(equal >= 0);
        if (!equal) {
                fprintf(stderr, "different output domains for array %s\n",
                                array_name);
                return -1;
        }

        out1 = isl_set_copy(dg1->out->domain);
        out2 = isl_set_copy(dg2->out->domain);
        out1 = isl_set_fix_dim_si(out1, dim1-1, array1);
        out2 = isl_set_fix_dim_si(out2, dim2-1, array2);
        out1 = isl_set_coalesce(out1);
        out2 = isl_set_coalesce(out2);
        isl_space *dim = isl_space_map_from_set(isl_set_get_space(out2));
        isl_map *id = isl_map_identity(dim);
        n_in = isl_map_dim(id, isl_dim_in);
        id = isl_map_remove_dims(id, isl_dim_in, n_in - 1, 1);
        id = isl_map_apply_domain(id, isl_map_copy(id));
        id = isl_map_intersect_domain(id, out1);
        id = isl_map_intersect_range(id, out2);

        isl_map *got = isl_map_empty(isl_map_get_space(id));
        isl_map *lost = isl_map_copy(got);
        check_equivalence_data data = { dg1, dg2, ec, got, lost };
        isl_map_foreach_basic_map(id, &basic_check, &data);
        isl_map_free(id);
        *proved = update_result(*proved,
                        array_name, data.got, array_dim, dim2 - array_dim);
        *not_proved = update_result(*not_proved,
                        array_name, data.lost, array_dim, dim2 - array_dim);
        return 0;
}

/* The input arrays of the two programs are supposed to be the same,
 * so they should at least have the same dimension.  Make sure
 * this is true, because we depend on it later on.
 */
static int check_input_arrays(dependence_graph *dg1, dependence_graph *dg2)
{
        for (int i = 0; i < dg1->input_computations.size(); ++i)
                for (int j = 0; j < dg2->input_computations.size(); ++j) {
                        if (strcmp(dg1->input_computations[i]->operation,
                                   dg2->input_computations[j]->operation))
                                continue;
                        if (dg1->input_computations[i]->dim ==
                            dg2->input_computations[j]->dim)
                                continue;
                        fprintf(stderr,
                            "input arrays \"%s\" do not have the same dimension\n",
                            dg1->input_computations[i]->operation);
                        return -1;
                }

        return 0;
}

static __isl_give isl_union_set *add_array(__isl_take isl_union_set *only,
        dependence_graph *dg, int array)
{
        isl_space *dim;
        isl_set *set;

        dim = isl_union_set_get_space(only);
        dim = isl_space_add_dims(dim, isl_dim_set, dg->output_array_dims[array]);
        dim = isl_space_set_tuple_name(dim, isl_dim_set, dg->output_arrays[array]);
        set = isl_set_universe(dim);
        only = isl_union_set_add_set(only, set);

        return only;
}

static void print_results(const char *str, __isl_keep isl_union_set *only)
{
        isl_printer *prn;

        if (isl_union_set_is_empty(only))
                return;

        fprintf(stdout, "%s: '", str);
        prn = isl_printer_to_file(isl_union_set_get_ctx(only), stdout);
        prn = isl_printer_print_union_set(prn, only);
        isl_printer_free(prn);
        fprintf(stdout, "'\n");
}

static int check_equivalence(isl_ctx *ctx,
        dependence_graph *dg1, dependence_graph *dg2, options *options)
{
        isl_space *dim;
        isl_set *context;
        isl_union_set *only1, *only2, *proved, *not_proved;
        dg1->flatten_associative_operators(options->ops->associative);
        dg2->flatten_associative_operators(options->ops->associative);
        equivalence_checker ec(options);
        int i1 = 0, i2 = 0;

        if (check_input_arrays(dg1, dg2))
                return -1;

        dim = isl_space_set_alloc(ctx, 0, 0);
        proved = isl_union_set_empty(isl_space_copy(dim));
        not_proved = isl_union_set_empty(isl_space_copy(dim));
        only1 = isl_union_set_empty(isl_space_copy(dim));
        only2 = isl_union_set_empty(dim);

        while (i1 < dg1->output_arrays.size() || i2 < dg2->output_arrays.size()) {
                int cmp;
                cmp = i1 == dg1->output_arrays.size() ? 1 :
                      i2 == dg2->output_arrays.size() ? -1 :
                      strcmp(dg1->output_arrays[i1], dg2->output_arrays[i2]);
                if (cmp < 0) {
                        only1 = add_array(only1, dg1, i1);
                        ++i1;
                } else if (cmp > 0) {
                        only2 = add_array(only2, dg2, i2);
                        ++i2;
                } else {
                        check_equivalence_array(ctx, &ec, dg1, dg2, i1, i2,
                                                &proved, &not_proved);
                        ++i1;
                        ++i2;
                }
        }

        context = isl_set_union(isl_set_copy(dg1->context),
                                isl_set_copy(dg2->context));
        proved = isl_union_set_gist_params(proved, isl_set_copy(context));
        not_proved = isl_union_set_gist_params(not_proved, context);

        print_results("Equivalence proved", proved);
        print_results("Equivalence NOT proved", not_proved);
        print_results("Only in program 1", only1);
        print_results("Only in program 2", only2);

        isl_union_set_free(proved);
        isl_union_set_free(not_proved);
        isl_union_set_free(only1);
        isl_union_set_free(only2);

        if (options->print_stats) {
                fprintf(stderr, "widenings: %d\n", ec.widenings);
                if (options->narrowing)
                        fprintf(stderr, "narrowings: %d\n", ec.narrowings);
        }

        return 0;
}

static void dump_vertex(FILE *out, computation *comp)
{
        fprintf(out, "ND_%p [label = \"%d,%s/%d\"];\n",
                comp, comp->location, comp->operation, comp->arity);
        for (int i = 0; i < comp->edges.size(); ++i)
                fprintf(out, "ND_%p -> ND_%p%s;\n",
                        comp, comp->edges[i]->source,
                        comp->edges[i]->type == edge::expansion ?
                                                " [color=\"blue\"]" : "");
}

static void dump_graph(FILE *out, dependence_graph *dg)
{
        fprintf(out, "digraph dummy {\n");
        dump_vertex(out, dg->out);
        for (int i = 0; i < dg->vertices.size(); ++i)
                dump_vertex(out, dg->vertices[i]);
        fprintf(out, "}\n");
}

static void dump_graphs(dependence_graph **dg, struct options *options)
{
        int i;
        char path[PATH_MAX];

        if (!options->dump_graphs)
                return;

        for (i = 0; i < 2; ++i) {
                FILE *out;
                int s;
                s = snprintf(path, sizeof(path), "%s.dot", options->program[i]);
                assert(s < sizeof(path));
                out = fopen(path, "w");
                assert(out);
                dump_graph(out, dg[i]);
                fclose(out);
        }
}

void parse_ops(struct options *options) {
        char *tok;

        if (options->associative) {
                tok = strtok(options->associative, ",");
                options->ops->associative.push_back(strdup(tok));
                while ((tok = strtok(NULL, ",")) != NULL)
                        options->ops->associative.push_back(strdup(tok));
        }

        if (options->commutative) {
                tok = strtok(options->commutative, ",");
                options->ops->commutative.push_back(strdup(tok));
                while ((tok = strtok(NULL, ",")) != NULL)
                        options->ops->commutative.push_back(strdup(tok));
        }
}

int main(int argc, char *argv[])
{
        struct options *options = options_new_with_defaults();
        struct isl_ctx *ctx;
        isl_set *context = NULL;
        dependence_graph *dg[2];

        argc = options_parse(options, argc, argv, ISL_ARG_ALL);
        parse_ops(options);

        ctx = isl_ctx_alloc_with_options(&options_args, options);
        if (!ctx) {
                fprintf(stderr, "Unable to allocate ctx\n");
                return -1;
        }

        if (options->context)
                context = isl_set_read_from_str(ctx, options->context);

        pdg::PDG *pdg[2];
        unsigned out_dim = 0;
        for (int i = 0; i < 2; ++i) {
                FILE *in;
                in = fopen(options->program[i], "r");
                if (!in) {
                        fprintf(stderr, "Unable to open %s\n", options->program[i]);
                        return -1;
                }
                pdg[i] = yaml::Load<pdg::PDG>(in, ctx);
                fclose(in);
                if (!pdg[i]) {
                        fprintf(stderr, "Unable to read %s\n", options->program[i]);
                        return -1;
                }
                out_dim = update_out_dim(pdg[i], out_dim);

                if (context &&
                    pdg[i]->params.size() != isl_set_dim(context, isl_dim_param)) {
                        fprintf(stdout,
                            "Parameter dimension mismatch; context ignored\n");
                        isl_set_free(context);
                        context = NULL;
                }
        }

        for (int i = 0; i < 2; ++i) {
                dg[i] = pdg_to_dg(pdg[i], out_dim, isl_set_copy(context));
                pdg[i]->free();
                delete pdg[i];
        }

        dump_graphs(dg, options);

        int res = check_equivalence(ctx, dg[0], dg[1], options);

        for (int i = 0; i < 2; ++i)
                delete dg[i];
        isl_set_free(context);
        isl_ctx_free(ctx);

        return res;
}