add isl_aff_mod_val
[platform/upstream/isl.git] / isl_schedule.c
index 7837295..dd7cac8 100644 (file)
@@ -1,16 +1,19 @@
 /*
  * Copyright 2011      INRIA Saclay
+ * Copyright 2012-2013 Ecole Normale Superieure
  *
- * Use of this software is governed by the GNU LGPLv2.1 license
+ * Use of this software is governed by the MIT license
  *
  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
  * 91893 Orsay, France
+ * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
  */
 
 #include <isl_ctx_private.h>
 #include <isl_map_private.h>
 #include <isl_space_private.h>
+#include <isl/aff.h>
 #include <isl/hash.h>
 #include <isl/constraint.h>
 #include <isl/schedule.h>
 #include <isl_tab.h>
 #include <isl_dim_map.h>
 #include <isl_hmap_map_basic_set.h>
-#include <isl_qsort.h>
+#include <isl_sort.h>
 #include <isl_schedule_private.h>
 #include <isl_band_private.h>
-#include <isl_list_private.h>
 #include <isl_options_private.h>
+#include <isl_tarjan.h>
 
 /*
  * The scheduling algorithm implemented in this file was inspired by
  * indicating whether the corresponding scheduling dimension results
  * in zero dependence distances within its band and with respect
  * to the proximity edges.
- *
- * index, min_index and on_stack are used during the SCC detection
- * index represents the order in which nodes are visited.
- * min_index is the index of the root of a (sub)component.
- * on_stack indicates whether the node is currently on the stack.
  */
 struct isl_sched_node {
        isl_space *dim;
@@ -82,11 +80,6 @@ struct isl_sched_node {
        int     *band;
        int     *band_id;
        int     *zero;
-
-       /* scc detection */
-       int      index;
-       int      min_index;
-       int      on_stack;
 };
 
 static int node_has_dim(const void *entry, const void *val)
@@ -126,6 +119,7 @@ struct isl_sched_edge {
 
 enum isl_edge_type {
        isl_edge_validity = 0,
+       isl_edge_first = isl_edge_validity,
        isl_edge_proximity,
        isl_edge_last = isl_edge_proximity
 };
@@ -141,6 +135,7 @@ enum isl_edge_type {
  * n is the number of nodes
  * node is the list of nodes
  * maxvar is the maximal number of variables over all nodes
+ * max_row is the allocated number of rows in the schedule
  * n_row is the current (maximal) number of linearly independent
  *     rows in the node schedules
  * n_total_row is the current number of rows in the node schedules
@@ -171,11 +166,7 @@ enum isl_edge_type {
  * src_scc and dst_scc are the source and sink SCCs of an edge with
  *     conflicting constraints
  *
- * scc, sp, index and stack are used during the detection of SCCs
- * scc is the number of the next SCC
- * stack contains the nodes on the path from the root to the current node
- * sp is the stack pointer
- * index is the index of the last node visited
+ * scc represents the number of components
  */
 struct isl_sched_graph {
        isl_hmap_map_basic_set *intra_hmap;
@@ -184,6 +175,7 @@ struct isl_sched_graph {
        struct isl_sched_node *node;
        int n;
        int maxvar;
+       int max_row;
        int n_row;
 
        int *sorted;
@@ -207,11 +199,7 @@ struct isl_sched_graph {
        int src_scc;
        int dst_scc;
 
-       /* scc detection */
        int scc;
-       int sp;
-       int index;
-       int *stack;
 };
 
 /* Initialize node_table based on the list of nodes.
@@ -339,7 +327,7 @@ static struct isl_sched_edge *graph_find_edge(struct isl_sched_graph *graph,
        return entry->data;
 }
 
-/* Check whether the dependence graph has an edge of the give type
+/* Check whether the dependence graph has an edge of the given type
  * between the given two nodes.
  */
 static int graph_has_edge(struct isl_sched_graph *graph,
@@ -367,10 +355,10 @@ static int graph_has_edge(struct isl_sched_graph *graph,
 static struct isl_sched_edge *graph_find_any_edge(struct isl_sched_graph *graph,
        struct isl_sched_node *src, struct isl_sched_node *dst)
 {
-       int i;
+       enum isl_edge_type i;
        struct isl_sched_edge *edge;
 
-       for (i = 0; i <= isl_edge_last; ++i) {
+       for (i = isl_edge_first; i <= isl_edge_last; ++i) {
                edge = graph_find_edge(graph, i, src, dst);
                if (edge)
                        return edge;
@@ -385,9 +373,9 @@ static void graph_remove_edge(struct isl_sched_graph *graph,
        struct isl_sched_edge *edge)
 {
        isl_ctx *ctx = isl_map_get_ctx(edge->map);
-       int i;
+       enum isl_edge_type i;
 
-       for (i = 0; i <= isl_edge_last; ++i) {
+       for (i = isl_edge_first; i <= isl_edge_last; ++i) {
                struct isl_hash_table_entry *entry;
 
                entry = graph_find_edge_entry(graph, i, edge->src, edge->dst);
@@ -405,10 +393,10 @@ static void graph_remove_edge(struct isl_sched_graph *graph,
 static int graph_has_any_edge(struct isl_sched_graph *graph,
        struct isl_sched_node *src, struct isl_sched_node *dst)
 {
-       int i;
+       enum isl_edge_type i;
        int r;
 
-       for (i = 0; i <= isl_edge_last; ++i) {
+       for (i = isl_edge_first; i <= isl_edge_last; ++i) {
                r = graph_has_edge(graph, i, src, dst);
                if (r < 0 || r)
                        return r;
@@ -436,15 +424,13 @@ static int graph_alloc(isl_ctx *ctx, struct isl_sched_graph *graph,
        graph->node = isl_calloc_array(ctx, struct isl_sched_node, graph->n);
        graph->sorted = isl_calloc_array(ctx, int, graph->n);
        graph->region = isl_alloc_array(ctx, struct isl_region, graph->n);
-       graph->stack = isl_alloc_array(ctx, int, graph->n);
        graph->edge = isl_calloc_array(ctx,
                                        struct isl_sched_edge, graph->n_edge);
 
        graph->intra_hmap = isl_hmap_map_basic_set_alloc(ctx, 2 * n_edge);
        graph->inter_hmap = isl_hmap_map_basic_set_alloc(ctx, 2 * n_edge);
 
-       if (!graph->node || !graph->region || !graph->stack || !graph->edge ||
-           !graph->sorted)
+       if (!graph->node || !graph->region || !graph->edge || !graph->sorted)
                return -1;
 
        for(i = 0; i < graph->n; ++i)
@@ -477,13 +463,49 @@ static void graph_free(isl_ctx *ctx, struct isl_sched_graph *graph)
                isl_map_free(graph->edge[i].map);
        free(graph->edge);
        free(graph->region);
-       free(graph->stack);
        for (i = 0; i <= isl_edge_last; ++i)
                isl_hash_table_free(ctx, graph->edge_table[i]);
        isl_hash_table_free(ctx, graph->node_table);
        isl_basic_set_free(graph->lp);
 }
 
+/* For each "set" on which this function is called, increment
+ * graph->n by one and update graph->maxvar.
+ */
+static int init_n_maxvar(__isl_take isl_set *set, void *user)
+{
+       struct isl_sched_graph *graph = user;
+       int nvar = isl_set_dim(set, isl_dim_set);
+
+       graph->n++;
+       if (nvar > graph->maxvar)
+               graph->maxvar = nvar;
+
+       isl_set_free(set);
+
+       return 0;
+}
+
+/* Compute the number of rows that should be allocated for the schedule.
+ * The graph can be split at most "n - 1" times, there can be at most
+ * two rows for each dimension in the iteration domains (in particular,
+ * we usually have one row, but it may be split by split_scaled),
+ * and there can be one extra row for ordering the statements.
+ * Note that if we have actually split "n - 1" times, then no ordering
+ * is needed, so in principle we could use "graph->n + 2 * graph->maxvar - 1".
+ */
+static int compute_max_row(struct isl_sched_graph *graph,
+       __isl_keep isl_union_set *domain)
+{
+       graph->n = 0;
+       graph->maxvar = 0;
+       if (isl_union_set_foreach_set(domain, &init_n_maxvar, graph) < 0)
+               return -1;
+       graph->max_row = graph->n + 2 * graph->maxvar;
+
+       return 0;
+}
+
 /* Add a new node to the graph representing the given set.
  */
 static int extract_node(__isl_take isl_set *set, void *user)
@@ -508,11 +530,11 @@ static int extract_node(__isl_take isl_set *set, void *user)
        graph->node[graph->n].nparam = nparam;
        graph->node[graph->n].sched = sched;
        graph->node[graph->n].sched_map = NULL;
-       band = isl_alloc_array(ctx, int, graph->n_edge + nvar);
+       band = isl_alloc_array(ctx, int, graph->max_row);
        graph->node[graph->n].band = band;
-       band_id = isl_calloc_array(ctx, int, graph->n_edge + nvar);
+       band_id = isl_calloc_array(ctx, int, graph->max_row);
        graph->node[graph->n].band_id = band_id;
-       zero = isl_calloc_array(ctx, int, graph->n_edge + nvar);
+       zero = isl_calloc_array(ctx, int, graph->max_row);
        graph->node[graph->n].zero = zero;
        graph->n++;
 
@@ -588,92 +610,60 @@ static int extract_edge(__isl_take isl_map *map, void *user)
        return graph_edge_table_add(ctx, graph, data->type, edge);
 }
 
-/* Check whether there is a validity dependence from src to dst,
- * forcing dst to follow src (if weak is not set).
- * If weak is set, then check if there is any dependence from src to dst.
+/* Check whether there is any dependence from node[j] to node[i]
+ * or from node[i] to node[j].
  */
-static int node_follows(struct isl_sched_graph *graph, 
-       struct isl_sched_node *dst, struct isl_sched_node *src, int weak)
+static int node_follows_weak(int i, int j, void *user)
 {
-       if (weak)
-               return graph_has_any_edge(graph, src, dst);
-       else
-               return graph_has_validity_edge(graph, src, dst);
+       int f;
+       struct isl_sched_graph *graph = user;
+
+       f = graph_has_any_edge(graph, &graph->node[j], &graph->node[i]);
+       if (f < 0 || f)
+               return f;
+       return graph_has_any_edge(graph, &graph->node[i], &graph->node[j]);
 }
 
-/* Perform Tarjan's algorithm for computing the strongly connected components
+/* Check whether there is a validity dependence from node[j] to node[i],
+ * forcing node[i] to follow node[j].
+ */
+static int node_follows_strong(int i, int j, void *user)
+{
+       struct isl_sched_graph *graph = user;
+
+       return graph_has_validity_edge(graph, &graph->node[j], &graph->node[i]);
+}
+
+/* Use Tarjan's algorithm for computing the strongly connected components
  * in the dependence graph (only validity edges).
  * If weak is set, we consider the graph to be undirected and
  * we effectively compute the (weakly) connected components.
  * Additionally, we also consider other edges when weak is set.
  */
-static int detect_sccs_tarjan(struct isl_sched_graph *g, int i, int weak)
+static int detect_ccs(isl_ctx *ctx, struct isl_sched_graph *graph, int weak)
 {
-       int j;
-
-       g->node[i].index = g->index;
-       g->node[i].min_index = g->index;
-       g->node[i].on_stack = 1;
-       g->index++;
-       g->stack[g->sp++] = i;
+       int i, n;
+       struct isl_tarjan_graph *g = NULL;
 
-       for (j = g->n - 1; j >= 0; --j) {
-               int f;
+       g = isl_tarjan_graph_init(ctx, graph->n,
+               weak ? &node_follows_weak : &node_follows_strong, graph);
+       if (!g)
+               return -1;
 
-               if (j == i)
-                       continue;
-               if (g->node[j].index >= 0 &&
-                       (!g->node[j].on_stack ||
-                        g->node[j].index > g->node[i].min_index))
-                       continue;
-               
-               f = node_follows(g, &g->node[i], &g->node[j], weak);
-               if (f < 0)
-                       return -1;
-               if (!f && weak) {
-                       f = node_follows(g, &g->node[j], &g->node[i], weak);
-                       if (f < 0)
-                               return -1;
+       graph->scc = 0;
+       i = 0;
+       n = graph->n;
+       while (n) {
+               while (g->order[i] != -1) {
+                       graph->node[g->order[i]].scc = graph->scc;
+                       --n;
+                       ++i;
                }
-               if (!f)
-                       continue;
-               if (g->node[j].index < 0) {
-                       detect_sccs_tarjan(g, j, weak);
-                       if (g->node[j].min_index < g->node[i].min_index)
-                               g->node[i].min_index = g->node[j].min_index;
-               } else if (g->node[j].index < g->node[i].min_index)
-                       g->node[i].min_index = g->node[j].index;
+               ++i;
+               graph->scc++;
        }
 
-       if (g->node[i].index != g->node[i].min_index)
-               return 0;
-
-       do {
-               j = g->stack[--g->sp];
-               g->node[j].on_stack = 0;
-               g->node[j].scc = g->scc;
-       } while (j != i);
-       g->scc++;
-
-       return 0;
-}
-
-static int detect_ccs(struct isl_sched_graph *graph, int weak)
-{
-       int i;
-
-       graph->index = 0;
-       graph->sp = 0;
-       graph->scc = 0;
-       for (i = graph->n - 1; i >= 0; --i)
-               graph->node[i].index = -1;
-
-       for (i = graph->n - 1; i >= 0; --i) {
-               if (graph->node[i].index >= 0)
-                       continue;
-               if (detect_sccs_tarjan(graph, i, weak) < 0)
-                       return -1;
-       }
+       isl_tarjan_graph_free(g);
 
        return 0;
 }
@@ -681,17 +671,17 @@ static int detect_ccs(struct isl_sched_graph *graph, int weak)
 /* Apply Tarjan's algorithm to detect the strongly connected components
  * in the dependence graph.
  */
-static int detect_sccs(struct isl_sched_graph *graph)
+static int detect_sccs(isl_ctx *ctx, struct isl_sched_graph *graph)
 {
-       return detect_ccs(graph, 0);
+       return detect_ccs(ctx, graph, 0);
 }
 
 /* Apply Tarjan's algorithm to detect the (weakly) connected components
  * in the dependence graph.
  */
-static int detect_wccs(struct isl_sched_graph *graph)
+static int detect_wccs(isl_ctx *ctx, struct isl_sched_graph *graph)
 {
-       return detect_ccs(graph, 1);
+       return detect_ccs(ctx, graph, 1);
 }
 
 static int cmp_scc(const void *a, const void *b, void *data)
@@ -705,9 +695,9 @@ static int cmp_scc(const void *a, const void *b, void *data)
 
 /* Sort the elements of graph->sorted according to the corresponding SCCs.
  */
-static void sort_sccs(struct isl_sched_graph *graph)
+static int sort_sccs(struct isl_sched_graph *graph)
 {
-       isl_quicksort(graph->sorted, graph->n, sizeof(int), &cmp_scc, graph);
+       return isl_sort(graph->sorted, graph->n, sizeof(int), &cmp_scc, graph);
 }
 
 /* Given a dependence relation R from a node to itself,
@@ -805,6 +795,8 @@ static int add_intra_validity_constraints(struct isl_sched_graph *graph,
 
        coef = isl_basic_set_transform_dims(coef, isl_dim_set,
                    isl_space_dim(dim, isl_dim_set), isl_mat_copy(node->cmap));
+       if (!coef)
+               goto error;
 
        total = isl_basic_set_total_dim(graph->lp);
        dim_map = isl_dim_map_alloc(ctx, total);
@@ -821,6 +813,9 @@ static int add_intra_validity_constraints(struct isl_sched_graph *graph,
        isl_space_free(dim);
 
        return 0;
+error:
+       isl_space_free(dim);
+       return -1;
 }
 
 /* Add constraints to graph->lp that force validity for the given
@@ -862,6 +857,8 @@ static int add_inter_validity_constraints(struct isl_sched_graph *graph,
        coef = isl_basic_set_transform_dims(coef, isl_dim_set,
                    isl_space_dim(dim, isl_dim_set) + src->nvar,
                    isl_mat_copy(dst->cmap));
+       if (!coef)
+               goto error;
 
        total = isl_basic_set_total_dim(graph->lp);
        dim_map = isl_dim_map_alloc(ctx, total);
@@ -891,10 +888,15 @@ static int add_inter_validity_constraints(struct isl_sched_graph *graph,
                        coef->n_eq, coef->n_ineq);
        graph->lp = isl_basic_set_add_constraints_dim_map(graph->lp,
                                                           coef, dim_map);
+       if (!graph->lp)
+               goto error;
        isl_space_free(dim);
        edge->end = graph->lp->n_ineq;
 
        return 0;
+error:
+       isl_space_free(dim);
+       return -1;
 }
 
 /* Add constraints to graph->lp that bound the dependence distance for the given
@@ -944,6 +946,8 @@ static int add_intra_proximity_constraints(struct isl_sched_graph *graph,
 
        coef = isl_basic_set_transform_dims(coef, isl_dim_set,
                    isl_space_dim(dim, isl_dim_set), isl_mat_copy(node->cmap));
+       if (!coef)
+               goto error;
 
        nparam = isl_space_dim(node->dim, isl_dim_param);
        total = isl_basic_set_total_dim(graph->lp);
@@ -964,6 +968,9 @@ static int add_intra_proximity_constraints(struct isl_sched_graph *graph,
        isl_space_free(dim);
 
        return 0;
+error:
+       isl_space_free(dim);
+       return -1;
 }
 
 /* Add constraints to graph->lp that bound the dependence distance for the given
@@ -1023,6 +1030,8 @@ static int add_inter_proximity_constraints(struct isl_sched_graph *graph,
        coef = isl_basic_set_transform_dims(coef, isl_dim_set,
                    isl_space_dim(dim, isl_dim_set) + src->nvar,
                    isl_mat_copy(dst->cmap));
+       if (!coef)
+               goto error;
 
        nparam = isl_space_dim(src->dim, isl_dim_param);
        total = isl_basic_set_total_dim(graph->lp);
@@ -1059,6 +1068,9 @@ static int add_inter_proximity_constraints(struct isl_sched_graph *graph,
        isl_space_free(dim);
 
        return 0;
+error:
+       isl_space_free(dim);
+       return -1;
 }
 
 static int add_all_validity_constraints(struct isl_sched_graph *graph)
@@ -1507,6 +1519,9 @@ static int update_schedule(struct isl_sched_graph *graph,
        if (sol->size == 0)
                isl_die(sol->ctx, isl_error_internal,
                        "no solution found", goto error);
+       if (graph->n_total_row >= graph->max_row)
+               isl_die(sol->ctx, isl_error_internal,
+                       "too many schedule rows", goto error);
 
        if (check_zero)
                zero = isl_int_is_zero(sol->el[1]) &&
@@ -1563,58 +1578,65 @@ error:
        return -1;
 }
 
-/* Convert node->sched into a map and return this map.
- * We simply add equality constraints that express each output variable
- * as the affine combination of parameters and input variables specified
- * by the schedule matrix.
- *
- * The result is cached in node->sched_map, which needs to be released
- * whenever node->sched is updated.
+/* Convert node->sched into a multi_aff and return this multi_aff.
  */
-static __isl_give isl_map *node_extract_schedule(struct isl_sched_node *node)
+static __isl_give isl_multi_aff *node_extract_schedule_multi_aff(
+       struct isl_sched_node *node)
 {
        int i, j;
-       isl_space *dim;
+       isl_space *space;
        isl_local_space *ls;
-       isl_basic_map *bmap;
-       isl_constraint *c;
+       isl_aff *aff;
+       isl_multi_aff *ma;
        int nrow, ncol;
        isl_int v;
 
-       if (node->sched_map)
-               return isl_map_copy(node->sched_map);
-
        nrow = isl_mat_rows(node->sched);
        ncol = isl_mat_cols(node->sched) - 1;
-       dim = isl_space_from_domain(isl_space_copy(node->dim));
-       dim = isl_space_add_dims(dim, isl_dim_out, nrow);
-       bmap = isl_basic_map_universe(isl_space_copy(dim));
-       ls = isl_local_space_from_space(dim);
+       space = isl_space_from_domain(isl_space_copy(node->dim));
+       space = isl_space_add_dims(space, isl_dim_out, nrow);
+       ma = isl_multi_aff_zero(space);
+       ls = isl_local_space_from_space(isl_space_copy(node->dim));
 
        isl_int_init(v);
 
        for (i = 0; i < nrow; ++i) {
-               c = isl_equality_alloc(isl_local_space_copy(ls));
-               isl_constraint_set_coefficient_si(c, isl_dim_out, i, -1);
+               aff = isl_aff_zero_on_domain(isl_local_space_copy(ls));
                isl_mat_get_element(node->sched, i, 0, &v);
-               isl_constraint_set_constant(c, v);
+               aff = isl_aff_set_constant(aff, v);
                for (j = 0; j < node->nparam; ++j) {
                        isl_mat_get_element(node->sched, i, 1 + j, &v);
-                       isl_constraint_set_coefficient(c, isl_dim_param, j, v);
+                       aff = isl_aff_set_coefficient(aff, isl_dim_param, j, v);
                }
                for (j = 0; j < node->nvar; ++j) {
                        isl_mat_get_element(node->sched,
                                            i, 1 + node->nparam + j, &v);
-                       isl_constraint_set_coefficient(c, isl_dim_in, j, v);
+                       aff = isl_aff_set_coefficient(aff, isl_dim_in, j, v);
                }
-               bmap = isl_basic_map_add_constraint(bmap, c);
+               ma = isl_multi_aff_set_aff(ma, i, aff);
        }
 
        isl_int_clear(v);
 
        isl_local_space_free(ls);
 
-       node->sched_map = isl_map_from_basic_map(bmap);
+       return ma;
+}
+
+/* Convert node->sched into a map and return this map.
+ *
+ * The result is cached in node->sched_map, which needs to be released
+ * whenever node->sched is updated.
+ */
+static __isl_give isl_map *node_extract_schedule(struct isl_sched_node *node)
+{
+       if (!node->sched_map) {
+               isl_multi_aff *ma;
+
+               ma = node_extract_schedule_multi_aff(node);
+               node->sched_map = isl_map_from_multi_aff(ma);
+       }
+
        return isl_map_copy(node->sched_map);
 }
 
@@ -1680,9 +1702,13 @@ static int sort_statements(isl_ctx *ctx, struct isl_sched_graph *graph)
        if (graph->n_edge == 0)
                return 0;
 
-       if (detect_sccs(graph) < 0)
+       if (detect_sccs(ctx, graph) < 0)
                return -1;
 
+       if (graph->n_total_row >= graph->max_row)
+               isl_die(ctx, isl_error_internal,
+                       "too many schedule rows", return -1);
+
        for (i = 0; i < graph->n; ++i) {
                struct isl_sched_node *node = &graph->node[i];
                int row = isl_mat_rows(node->sched);
@@ -1736,10 +1762,18 @@ static __isl_give isl_schedule *extract_schedule(struct isl_sched_graph *graph,
                int r, b;
                int *band_end, *band_id, *zero;
 
+               sched->node[i].sched =
+                       node_extract_schedule_multi_aff(&graph->node[i]);
+               if (!sched->node[i].sched)
+                       goto error;
+
+               sched->node[i].n_band = graph->n_band;
+               if (graph->n_band == 0)
+                       continue;
+
                band_end = isl_alloc_array(ctx, int, graph->n_band);
                band_id = isl_alloc_array(ctx, int, graph->n_band);
                zero = isl_alloc_array(ctx, int, graph->n_total_row);
-               sched->node[i].sched = node_extract_schedule(&graph->node[i]);
                sched->node[i].band_end = band_end;
                sched->node[i].band_id = band_id;
                sched->node[i].zero = zero;
@@ -1809,7 +1843,7 @@ static int copy_edges(isl_ctx *ctx, struct isl_sched_graph *dst,
        int (*edge_pred)(struct isl_sched_edge *edge, int data), int data)
 {
        int i;
-       int t;
+       enum isl_edge_type t;
 
        dst->n_edge = 0;
        for (i = 0; i < src->n_edge; ++i) {
@@ -1841,7 +1875,7 @@ static int copy_edges(isl_ctx *ctx, struct isl_sched_graph *dst,
                dst->edge[dst->n_edge].proximity = edge->proximity;
                dst->n_edge++;
 
-               for (t = 0; t <= isl_edge_last; ++t) {
+               for (t = isl_edge_first; t <= isl_edge_last; ++t) {
                        if (edge !=
                            graph_find_edge(src, t, edge->src, edge->dst))
                                continue;
@@ -1876,6 +1910,7 @@ static int copy_schedule(struct isl_sched_graph *dst,
                src->n++;
        }
 
+       dst->max_row = src->max_row;
        dst->n_total_row = src->n_total_row;
        dst->n_band = src->n_band;
 
@@ -1943,6 +1978,7 @@ static int compute_sub_schedule(isl_ctx *ctx,
        if (copy_edges(ctx, &split, graph, edge_pred, data) < 0)
                goto error;
        split.n_row = graph->n_row;
+       split.max_row = graph->max_row;
        split.n_total_row = graph->n_total_row;
        split.n_band = graph->n_band;
        split.band_start = graph->band_start;
@@ -2044,6 +2080,10 @@ static int compute_split_schedule(isl_ctx *ctx, struct isl_sched_graph *graph)
        int n_band, orig_band;
        int drop;
 
+       if (graph->n_total_row >= graph->max_row)
+               isl_die(ctx, isl_error_internal,
+                       "too many schedule rows", return -1);
+
        drop = graph->n_total_row - graph->band_start;
        graph->n_total_row -= drop;
        graph->n_row -= drop;
@@ -2151,6 +2191,8 @@ static int add_intra_constraints(struct isl_sched_graph *graph,
        struct isl_sched_node *node = edge->src;
 
        coef = intra_coefficients(graph, map);
+       if (!coef)
+               return -1;
 
        dim = isl_space_domain(isl_space_unwrap(isl_basic_set_get_space(coef)));
 
@@ -2199,6 +2241,8 @@ static int add_inter_constraints(struct isl_sched_graph *graph,
        struct isl_sched_node *dst = edge->dst;
 
        coef = inter_coefficients(graph, map);
+       if (!coef)
+               return -1;
 
        dim = isl_space_domain(isl_space_unwrap(isl_basic_set_get_space(coef)));
 
@@ -2428,6 +2472,10 @@ static int split_scaled(isl_ctx *ctx, struct isl_sched_graph *graph)
        if (graph->n <= 1)
                return 0;
 
+       if (graph->n_total_row >= graph->max_row)
+               isl_die(ctx, isl_error_internal,
+                       "too many schedule rows", return -1);
+
        isl_int_init(gcd);
        isl_int_init(gcd_i);
 
@@ -2481,8 +2529,61 @@ error:
        return -1;
 }
 
+static int compute_component_schedule(isl_ctx *ctx,
+       struct isl_sched_graph *graph);
+
+/* Is the schedule row "sol" trivial on node "node"?
+ * That is, is the solution zero on the dimensions orthogonal to
+ * the previously found solutions?
+ * Each coefficient is represented as the difference between
+ * two non-negative values in "sol".  The coefficient is then
+ * zero if those two values are equal to each other.
+ */
+static int is_trivial(struct isl_sched_node *node, __isl_keep isl_vec *sol)
+{
+       int i;
+       int pos;
+       int len;
+
+       pos = 1 + node->start + 1 + 2 * (node->nparam + node->rank);
+       len = 2 * (node->nvar - node->rank);
+
+       if (len == 0)
+               return 0;
+
+       for (i = 0; i < len; i += 2)
+               if (isl_int_ne(sol->el[pos + i], sol->el[pos + i + 1]))
+                       return 0;
+
+       return 1;
+}
+
+/* Is the schedule row "sol" trivial on any node where it should
+ * not be trivial?
+ */
+static int is_any_trivial(struct isl_sched_graph *graph,
+       __isl_keep isl_vec *sol)
+{
+       int i;
+
+       for (i = 0; i < graph->n; ++i) {
+               struct isl_sched_node *node = &graph->node[i];
+
+               if (!needs_row(graph, node))
+                       continue;
+               if (is_trivial(node, sol))
+                       return 1;
+       }
+
+       return 0;
+}
+
 /* Construct a schedule row for each node such that as many dependences
  * as possible are carried and then continue with the next band.
+ *
+ * If the computed schedule row turns out to be trivial on one or
+ * more nodes where it should not be trivial, then we throw it away
+ * and try again on each component separately.
  */
 static int carry_dependences(isl_ctx *ctx, struct isl_sched_graph *graph)
 {
@@ -2509,12 +2610,21 @@ static int carry_dependences(isl_ctx *ctx, struct isl_sched_graph *graph)
                        "error in schedule construction", return -1);
        }
 
+       isl_int_divexact(sol->el[1], sol->el[1], sol->el[0]);
        if (isl_int_cmp_si(sol->el[1], n_edge) >= 0) {
                isl_vec_free(sol);
                isl_die(ctx, isl_error_unknown,
                        "unable to carry dependences", return -1);
        }
 
+       if (is_any_trivial(graph, sol)) {
+               isl_vec_free(sol);
+               if (graph->scc > 1)
+                       return compute_component_schedule(ctx, graph);
+               isl_die(ctx, isl_error_unknown,
+                       "unable to construct non-trivial solution", return -1);
+       }
+
        if (update_schedule(graph, sol, 0, 0) < 0)
                return -1;
 
@@ -2601,9 +2711,10 @@ static int compute_schedule_wcc(isl_ctx *ctx, struct isl_sched_graph *graph)
 {
        int force_zero = 0;
 
-       if (detect_sccs(graph) < 0)
+       if (detect_sccs(ctx, graph) < 0)
+               return -1;
+       if (sort_sccs(graph) < 0)
                return -1;
-       sort_sccs(graph);
 
        if (compute_maxvar(graph) < 0)
                return -1;
@@ -2649,10 +2760,14 @@ static int compute_schedule_wcc(isl_ctx *ctx, struct isl_sched_graph *graph)
 /* Add a row to the schedules that separates the SCCs and move
  * to the next band.
  */
-static int split_on_scc(struct isl_sched_graph *graph)
+static int split_on_scc(isl_ctx *ctx, struct isl_sched_graph *graph)
 {
        int i;
 
+       if (graph->n_total_row >= graph->max_row)
+               isl_die(ctx, isl_error_internal,
+                       "too many schedule rows", return -1);
+
        for (i = 0; i < graph->n; ++i) {
                struct isl_sched_node *node = &graph->node[i];
                int row = isl_mat_rows(node->sched);
@@ -2690,8 +2805,10 @@ static int compute_component_schedule(isl_ctx *ctx,
        int n_total_row, orig_total_row;
        int n_band, orig_band;
 
-       if (ctx->opt->schedule_fuse == ISL_SCHEDULE_FUSE_MIN)
-               split_on_scc(graph);
+       if (ctx->opt->schedule_fuse == ISL_SCHEDULE_FUSE_MIN ||
+           ctx->opt->schedule_separate_components)
+               if (split_on_scc(ctx, graph) < 0)
+                       return -1;
 
        n_total_row = 0;
        orig_total_row = graph->n_total_row;
@@ -2738,10 +2855,10 @@ static int compute_component_schedule(isl_ctx *ctx,
 static int compute_schedule(isl_ctx *ctx, struct isl_sched_graph *graph)
 {
        if (ctx->opt->schedule_fuse == ISL_SCHEDULE_FUSE_MIN) {
-               if (detect_sccs(graph) < 0)
+               if (detect_sccs(ctx, graph) < 0)
                        return -1;
        } else {
-               if (detect_wccs(graph) < 0)
+               if (detect_wccs(ctx, graph) < 0)
                        return -1;
        }
 
@@ -2787,6 +2904,8 @@ __isl_give isl_schedule *isl_union_set_compute_schedule(
        if (graph_alloc(ctx, &graph, graph.n,
            isl_union_map_n_map(validity) + isl_union_map_n_map(proximity)) < 0)
                goto error;
+       if (compute_max_row(&graph, domain) < 0)
+               goto error;
        graph.root = 1;
        graph.n = 0;
        if (isl_union_set_foreach_set(domain, &extract_node, &graph) < 0)
@@ -2836,7 +2955,7 @@ void *isl_schedule_free(__isl_take isl_schedule *sched)
                return NULL;
 
        for (i = 0; i < sched->n; ++i) {
-               isl_map_free(sched->node[i].sched);
+               isl_multi_aff_free(sched->node[i].sched);
                free(sched->node[i].band_end);
                free(sched->node[i].band_id);
                free(sched->node[i].zero);
@@ -2852,6 +2971,84 @@ isl_ctx *isl_schedule_get_ctx(__isl_keep isl_schedule *schedule)
        return schedule ? isl_space_get_ctx(schedule->dim) : NULL;
 }
 
+/* Set max_out to the maximal number of output dimensions over
+ * all maps.
+ */
+static int update_max_out(__isl_take isl_map *map, void *user)
+{
+       int *max_out = user;
+       int n_out = isl_map_dim(map, isl_dim_out);
+
+       if (n_out > *max_out)
+               *max_out = n_out;
+
+       isl_map_free(map);
+       return 0;
+}
+
+/* Internal data structure for map_pad_range.
+ *
+ * "max_out" is the maximal schedule dimension.
+ * "res" collects the results.
+ */
+struct isl_pad_schedule_map_data {
+       int max_out;
+       isl_union_map *res;
+};
+
+/* Pad the range of the given map with zeros to data->max_out and
+ * then add the result to data->res.
+ */
+static int map_pad_range(__isl_take isl_map *map, void *user)
+{
+       struct isl_pad_schedule_map_data *data = user;
+       int i;
+       int n_out = isl_map_dim(map, isl_dim_out);
+
+       map = isl_map_add_dims(map, isl_dim_out, data->max_out - n_out);
+       for (i = n_out; i < data->max_out; ++i)
+               map = isl_map_fix_si(map, isl_dim_out, i, 0);
+
+       data->res = isl_union_map_add_map(data->res, map);
+       if (!data->res)
+               return -1;
+
+       return 0;
+}
+
+/* Pad the ranges of the maps in the union map with zeros such they all have
+ * the same dimension.
+ */
+static __isl_give isl_union_map *pad_schedule_map(
+       __isl_take isl_union_map *umap)
+{
+       struct isl_pad_schedule_map_data data;
+
+       if (!umap)
+               return NULL;
+       if (isl_union_map_n_map(umap) <= 1)
+               return umap;
+
+       data.max_out = 0;
+       if (isl_union_map_foreach_map(umap, &update_max_out, &data.max_out) < 0)
+               return isl_union_map_free(umap);
+
+       data.res = isl_union_map_empty(isl_union_map_get_space(umap));
+       if (isl_union_map_foreach_map(umap, &map_pad_range, &data) < 0)
+               data.res = isl_union_map_free(data.res);
+
+       isl_union_map_free(umap);
+       return data.res;
+}
+
+/* Return an isl_union_map of the schedule.  If we have already constructed
+ * a band forest, then this band forest may have been modified so we need
+ * to extract the isl_union_map from the forest rather than from
+ * the originally computed schedule.  This reconstructed schedule map
+ * then needs to be padded with zeros to unify the schedule space
+ * since the result of isl_band_list_get_suffix_schedule may not have
+ * a unified schedule space.
+ */
 __isl_give isl_union_map *isl_schedule_get_map(__isl_keep isl_schedule *sched)
 {
        int i;
@@ -2860,10 +3057,18 @@ __isl_give isl_union_map *isl_schedule_get_map(__isl_keep isl_schedule *sched)
        if (!sched)
                return NULL;
 
+       if (sched->band_forest) {
+               umap = isl_band_list_get_suffix_schedule(sched->band_forest);
+               return pad_schedule_map(umap);
+       }
+
        umap = isl_union_map_empty(isl_space_copy(sched->dim));
-       for (i = 0; i < sched->n; ++i)
-               umap = isl_union_map_add_map(umap,
-                                           isl_map_copy(sched->node[i].sched));
+       for (i = 0; i < sched->n; ++i) {
+               isl_multi_aff *ma;
+
+               ma = isl_multi_aff_copy(sched->node[i].sched);
+               umap = isl_union_map_add_map(umap, isl_map_from_multi_aff(ma));
+       }
 
        return umap;
 }
@@ -2880,6 +3085,10 @@ static __isl_give isl_band_list *construct_band_list(
  * Because of the way the schedule is constructed, we know that
  * the position of the band inside the schedule of a node is the same
  * for all active nodes.
+ *
+ * The partial schedule for the band is created before the children
+ * are created to that construct_band_list can refer to the partial
+ * schedule of the parent.
  */
 static __isl_give isl_band *construct_band(__isl_keep isl_schedule *schedule,
        __isl_keep isl_band *parent,
@@ -2890,26 +3099,14 @@ static __isl_give isl_band *construct_band(__isl_keep isl_schedule *schedule,
        isl_band *band;
        unsigned start, end;
 
-       band = isl_calloc_type(ctx, isl_band);
+       band = isl_band_alloc(ctx);
        if (!band)
                return NULL;
 
-       band->ref = 1;
        band->schedule = schedule;
        band->parent = parent;
 
        for (i = 0; i < schedule->n; ++i)
-               if (active[i] && schedule->node[i].n_band > band_nr + 1)
-                       break;
-
-       if (i < schedule->n) {
-               band->children = construct_band_list(schedule, band,
-                                               band_nr + 1, active, n_active);
-               if (!band->children)
-                       goto error;
-       }
-
-       for (i = 0; i < schedule->n; ++i)
                if (active[i])
                        break;
 
@@ -2929,30 +3126,265 @@ static __isl_give isl_band *construct_band(__isl_keep isl_schedule *schedule,
        for (j = 0; j < band->n; ++j)
                band->zero[j] = schedule->node[i].zero[start + j];
 
-       band->map = isl_union_map_empty(isl_space_copy(schedule->dim));
+       band->pma = isl_union_pw_multi_aff_empty(isl_space_copy(schedule->dim));
        for (i = 0; i < schedule->n; ++i) {
-               isl_map *map;
+               isl_multi_aff *ma;
+               isl_pw_multi_aff *pma;
                unsigned n_out;
 
                if (!active[i])
                        continue;
 
-               map = isl_map_copy(schedule->node[i].sched);
-               n_out = isl_map_dim(map, isl_dim_out);
-               map = isl_map_project_out(map, isl_dim_out, end, n_out - end);
-               map = isl_map_project_out(map, isl_dim_out, 0, start);
-               band->map = isl_union_map_union(band->map,
-                                               isl_union_map_from_map(map));
+               ma = isl_multi_aff_copy(schedule->node[i].sched);
+               n_out = isl_multi_aff_dim(ma, isl_dim_out);
+               ma = isl_multi_aff_drop_dims(ma, isl_dim_out, end, n_out - end);
+               ma = isl_multi_aff_drop_dims(ma, isl_dim_out, 0, start);
+               pma = isl_pw_multi_aff_from_multi_aff(ma);
+               band->pma = isl_union_pw_multi_aff_add_pw_multi_aff(band->pma,
+                                                                   pma);
        }
-       if (!band->map)
+       if (!band->pma)
                goto error;
 
+       for (i = 0; i < schedule->n; ++i)
+               if (active[i] && schedule->node[i].n_band > band_nr + 1)
+                       break;
+
+       if (i < schedule->n) {
+               band->children = construct_band_list(schedule, band,
+                                               band_nr + 1, active, n_active);
+               if (!band->children)
+                       goto error;
+       }
+
        return band;
 error:
        isl_band_free(band);
        return NULL;
 }
 
+/* Internal data structure used inside cmp_band and pw_multi_aff_extract_int.
+ *
+ * r is set to a negative value if anything goes wrong.
+ *
+ * c1 stores the result of extract_int.
+ * c2 is a temporary value used inside cmp_band_in_ancestor.
+ * t is a temporary value used inside extract_int.
+ *
+ * first and equal are used inside extract_int.
+ * first is set if we are looking at the first isl_multi_aff inside
+ * the isl_union_pw_multi_aff.
+ * equal is set if all the isl_multi_affs have been equal so far.
+ */
+struct isl_cmp_band_data {
+       int r;
+
+       int first;
+       int equal;
+
+       isl_int t;
+       isl_int c1;
+       isl_int c2;
+};
+
+/* Check if "ma" assigns a constant value.
+ * Note that this function is only called on isl_multi_affs
+ * with a single output dimension.
+ *
+ * If "ma" assigns a constant value then we compare it to data->c1
+ * or assign it to data->c1 if this is the first isl_multi_aff we consider.
+ * If "ma" does not assign a constant value or if it assigns a value
+ * that is different from data->c1, then we set data->equal to zero
+ * and terminate the check.
+ */
+static int multi_aff_extract_int(__isl_take isl_set *set,
+       __isl_take isl_multi_aff *ma, void *user)
+{
+       isl_aff *aff;
+       struct isl_cmp_band_data *data = user;
+
+       aff = isl_multi_aff_get_aff(ma, 0);
+       data->r = isl_aff_is_cst(aff);
+       if (data->r >= 0 && data->r) {
+               isl_aff_get_constant(aff, &data->t);
+               if (data->first) {
+                       isl_int_set(data->c1, data->t);
+                       data->first = 0;
+               } else if (!isl_int_eq(data->c1, data->t))
+                       data->equal = 0;
+       } else if (data->r >= 0 && !data->r)
+               data->equal = 0;
+
+       isl_aff_free(aff);
+       isl_set_free(set);
+       isl_multi_aff_free(ma);
+
+       if (data->r < 0)
+               return -1;
+       if (!data->equal)
+               return -1;
+       return 0;
+}
+
+/* This function is called for each isl_pw_multi_aff in
+ * the isl_union_pw_multi_aff checked by extract_int.
+ * Check all the isl_multi_affs inside "pma".
+ */
+static int pw_multi_aff_extract_int(__isl_take isl_pw_multi_aff *pma,
+       void *user)
+{
+       int r;
+
+       r = isl_pw_multi_aff_foreach_piece(pma, &multi_aff_extract_int, user);
+       isl_pw_multi_aff_free(pma);
+
+       return r;
+}
+
+/* Check if "upma" assigns a single constant value to its domain.
+ * If so, return 1 and store the result in data->c1.
+ * If not, return 0.
+ *
+ * A negative return value from isl_union_pw_multi_aff_foreach_pw_multi_aff
+ * means that either an error occurred or that we have broken off the check
+ * because we already know the result is going to be negative.
+ * In the latter case, data->equal is set to zero.
+ */
+static int extract_int(__isl_keep isl_union_pw_multi_aff *upma,
+       struct isl_cmp_band_data *data)
+{
+       data->first = 1;
+       data->equal = 1;
+
+       if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
+                                       &pw_multi_aff_extract_int, data) < 0) {
+               if (!data->equal)
+                       return 0;
+               return -1;
+       }
+
+       return !data->first && data->equal;
+}
+
+/* Compare "b1" and "b2" based on the parent schedule of their ancestor
+ * "ancestor".
+ *
+ * If the parent of "ancestor" also has a single member, then we
+ * first try to compare the two band based on the partial schedule
+ * of this parent.
+ *
+ * Otherwise, or if the result is inconclusive, we look at the partial schedule
+ * of "ancestor" itself.
+ * In particular, we specialize the parent schedule based
+ * on the domains of the child schedules, check if both assign
+ * a single constant value and, if so, compare the two constant values.
+ * If the specialized parent schedules do not assign a constant value,
+ * then they cannot be used to order the two bands and so in this case
+ * we return 0.
+ */
+static int cmp_band_in_ancestor(__isl_keep isl_band *b1,
+       __isl_keep isl_band *b2, struct isl_cmp_band_data *data,
+       __isl_keep isl_band *ancestor)
+{
+       isl_union_pw_multi_aff *upma;
+       isl_union_set *domain;
+       int r;
+
+       if (data->r < 0)
+               return 0;
+
+       if (ancestor->parent && ancestor->parent->n == 1) {
+               r = cmp_band_in_ancestor(b1, b2, data, ancestor->parent);
+               if (data->r < 0)
+                       return 0;
+               if (r)
+                       return r;
+       }
+
+       upma = isl_union_pw_multi_aff_copy(b1->pma);
+       domain = isl_union_pw_multi_aff_domain(upma);
+       upma = isl_union_pw_multi_aff_copy(ancestor->pma);
+       upma = isl_union_pw_multi_aff_intersect_domain(upma, domain);
+       r = extract_int(upma, data);
+       isl_union_pw_multi_aff_free(upma);
+
+       if (r < 0)
+               data->r = -1;
+       if (r < 0 || !r)
+               return 0;
+
+       isl_int_set(data->c2, data->c1);
+
+       upma = isl_union_pw_multi_aff_copy(b2->pma);
+       domain = isl_union_pw_multi_aff_domain(upma);
+       upma = isl_union_pw_multi_aff_copy(ancestor->pma);
+       upma = isl_union_pw_multi_aff_intersect_domain(upma, domain);
+       r = extract_int(upma, data);
+       isl_union_pw_multi_aff_free(upma);
+
+       if (r < 0)
+               data->r = -1;
+       if (r < 0 || !r)
+               return 0;
+
+       return isl_int_cmp(data->c2, data->c1);
+}
+
+/* Compare "a" and "b" based on the parent schedule of their parent.
+ */
+static int cmp_band(const void *a, const void *b, void *user)
+{
+       isl_band *b1 = *(isl_band * const *) a;
+       isl_band *b2 = *(isl_band * const *) b;
+       struct isl_cmp_band_data *data = user;
+
+       return cmp_band_in_ancestor(b1, b2, data, b1->parent);
+}
+
+/* Sort the elements in "list" based on the partial schedules of its parent
+ * (and ancestors).  In particular if the parent assigns constant values
+ * to the domains of the bands in "list", then the elements are sorted
+ * according to that order.
+ * This order should be a more "natural" order for the user, but otherwise
+ * shouldn't have any effect.
+ * If we would be constructing an isl_band forest directly in
+ * isl_union_set_compute_schedule then there wouldn't be any need
+ * for a reordering, since the children would be added to the list
+ * in their natural order automatically.
+ *
+ * If there is only one element in the list, then there is no need to sort
+ * anything.
+ * If the partial schedule of the parent has more than one member
+ * (or if there is no parent), then it's
+ * defnitely not assigning constant values to the different children in
+ * the list and so we wouldn't be able to use it to sort the list.
+ */
+static __isl_give isl_band_list *sort_band_list(__isl_take isl_band_list *list,
+       __isl_keep isl_band *parent)
+{
+       struct isl_cmp_band_data data;
+
+       if (!list)
+               return NULL;
+       if (list->n <= 1)
+               return list;
+       if (!parent || parent->n != 1)
+               return list;
+
+       data.r = 0;
+       isl_int_init(data.c1);
+       isl_int_init(data.c2);
+       isl_int_init(data.t);
+       isl_sort(list->p, list->n, sizeof(list->p[0]), &cmp_band, &data);
+       if (data.r < 0)
+               list = isl_band_list_free(list);
+       isl_int_clear(data.c1);
+       isl_int_clear(data.c2);
+       isl_int_clear(data.t);
+
+       return list;
+}
+
 /* Construct a list of bands that start at the same position (with
  * sequence number band_nr) in the schedules of the nodes that
  * were active in the parent band.
@@ -3037,6 +3469,8 @@ static __isl_give isl_band_list *construct_band_list(
 
        free(active);
 
+       list = sort_band_list(list, parent);
+
        return list;
 }
 
@@ -3077,6 +3511,24 @@ __isl_give isl_band_list *isl_schedule_get_band_forest(
        return isl_band_list_dup(schedule->band_forest);
 }
 
+/* Call "fn" on each band in the schedule in depth-first post-order.
+ */
+int isl_schedule_foreach_band(__isl_keep isl_schedule *sched,
+       int (*fn)(__isl_keep isl_band *band, void *user), void *user)
+{
+       int r;
+       isl_band_list *forest;
+
+       if (!sched)
+               return -1;
+
+       forest = isl_schedule_get_band_forest(sched);
+       r = isl_band_list_foreach_band(forest, fn, user);
+       isl_band_list_free(forest);
+
+       return r;
+}
+
 static __isl_give isl_printer *print_band_list(__isl_take isl_printer *p,
        __isl_keep isl_band_list *list);
 
@@ -3086,7 +3538,7 @@ static __isl_give isl_printer *print_band(__isl_take isl_printer *p,
        isl_band_list *children;
 
        p = isl_printer_start_line(p);
-       p = isl_printer_print_union_map(p, band->map);
+       p = isl_printer_print_union_pw_multi_aff(p, band->pma);
        p = isl_printer_end_line(p);
 
        if (!isl_band_has_children(band))