update pet to version 0.11.8

[ppcg.git] / hybrid.c

/*
 * Copyright 2013      Ecole Normale Superieure
 * Copyright 2015      Sven Verdoolaege
 *
 * Use of this software is governed by the MIT license
 *
 * Written by Sven Verdoolaege,
 * Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
 */

#include <string.h>

#include <isl/space.h>
#include <isl/constraint.h>
#include <isl/val.h>
#include <isl/aff.h>
#include <isl/set.h>
#include <isl/map.h>
#include <isl/union_set.h>
#include <isl/union_map.h>

#include "hybrid.h"
#include "schedule.h"

/* The hybrid tiling implemented in this file is based on
 * Grosser et al., "Hybrid Hexagonal/Classical Tiling for GPUs".
 */

/* Bounds on relative dependence distances in input to hybrid tiling.
 * upper is an upper bound on the relative dependence distances
 * in the first space dimension
 * -lower is a lower bound on the relative dependence distances
 * in all space dimensions.
 *
 * In particular,
 *
 *      d_i >= -lower_i d_0
 * and
 *      d_1 <= upper d_0
 *
 * for each dependence distance vector d, where d_1 is the component
 * corresponding to the first space dimension.
 *
 * upper and lower are always non-negative.
 * Some of the values may be NaN if no bound could be found.
 */
struct ppcg_ht_bounds {
        isl_val *upper;
        isl_multi_val *lower;
};

/* Free "bounds" along with all its fields.
 */
__isl_null ppcg_ht_bounds *ppcg_ht_bounds_free(
        __isl_take ppcg_ht_bounds *bounds)
{
        if (!bounds)
                return NULL;
        isl_val_free(bounds->upper);
        isl_multi_val_free(bounds->lower);
        free(bounds);

        return NULL;
}

/* Create a ppcg_ht_bounds object for a band living in "space".
 * The bounds are initialized to NaN.
 */
__isl_give ppcg_ht_bounds *ppcg_ht_bounds_alloc(__isl_take isl_space *space)
{
        int i, n;
        isl_ctx *ctx;
        ppcg_ht_bounds *bounds;

        if (!space)
                return NULL;

        ctx = isl_space_get_ctx(space);
        bounds = isl_alloc_type(ctx, struct ppcg_ht_bounds);
        if (!bounds)
                goto error;
        bounds->upper = isl_val_nan(ctx);
        bounds->lower = isl_multi_val_zero(space);
        n = isl_multi_val_dim(bounds->lower, isl_dim_set);
        for (i = 0; i < n; ++i) {
                isl_val *v = isl_val_copy(bounds->upper);
                bounds->lower = isl_multi_val_set_val(bounds->lower, i, v);
        }

        if (!bounds->lower || !bounds->upper)
                return ppcg_ht_bounds_free(bounds);

        return bounds;
error:
        isl_space_free(space);
        return NULL;
}

void ppcg_ht_bounds_dump(__isl_keep ppcg_ht_bounds *bounds)
{
        if (!bounds)
                return;

        fprintf(stderr, "lower: ");
        isl_multi_val_dump(bounds->lower);
        fprintf(stderr, "upper: ");
        isl_val_dump(bounds->upper);
}

/* Return the upper bound on the relative dependence distances
 * in the first space dimension.
 */
__isl_give isl_val *ppcg_ht_bounds_get_upper(__isl_keep ppcg_ht_bounds *bounds)
{
        if (!bounds)
                return NULL;
        return isl_val_copy(bounds->upper);
}

/* Replace the upper bound on the relative dependence distances
 * in the first space dimension by "upper".
 */
__isl_give ppcg_ht_bounds *ppcg_ht_bounds_set_upper(
        __isl_take ppcg_ht_bounds *bounds, __isl_take isl_val *upper)
{
        if (!bounds || !upper)
                goto error;
        isl_val_free(bounds->upper);
        bounds->upper = upper;
        return bounds;
error:
        ppcg_ht_bounds_free(bounds);
        isl_val_free(upper);
        return NULL;
}

/* Return the lower bound on the relative dependence distances
 * in space dimension "pos".
 */
__isl_give isl_val *ppcg_ht_bounds_get_lower(__isl_keep ppcg_ht_bounds *bounds,
        int pos)
{
        if (!bounds)
                return NULL;
        return isl_multi_val_get_val(bounds->lower, pos);
}

/* Replace the lower bound on the relative dependence distances
 * in space dimension "pos" by "lower".
 */
__isl_give ppcg_ht_bounds *ppcg_ht_bounds_set_lower(
        __isl_take ppcg_ht_bounds *bounds, int pos, __isl_take isl_val *lower)
{
        if (!bounds || !lower)
                goto error;
        bounds->lower = isl_multi_val_set_val(bounds->lower, pos, lower);
        if (!bounds->lower)
                return ppcg_ht_bounds_free(bounds);
        return bounds;
error:
        ppcg_ht_bounds_free(bounds);
        isl_val_free(lower);
        return NULL;
}

/* Can the bounds on relative dependence distances recorded in "bounds"
 * be used to perform hybrid tiling?
 * In particular, have appropriate lower and upper bounds been found?
 * Any NaN indicates that no corresponding bound was found.
 */
isl_bool ppcg_ht_bounds_is_valid(__isl_keep ppcg_ht_bounds *bounds)
{
        isl_bool is_nan;
        int i, n;

        if (!bounds)
                return isl_bool_error;
        is_nan = isl_val_is_nan(bounds->upper);
        if (is_nan < 0)
                return isl_bool_error;
        if (is_nan)
                return isl_bool_false;

        n = isl_multi_val_dim(bounds->lower, isl_dim_set);
        for (i = 0; i < n; ++i) {
                isl_val *v;

                v = isl_multi_val_get_val(bounds->lower, i);
                is_nan = isl_val_is_nan(v);
                if (is_nan < 0)
                        return isl_bool_error;
                if (is_nan)
                        return isl_bool_false;
                isl_val_free(v);
        }

        return isl_bool_true;
}

/* Structure that represents the basic hexagonal tiling,
 * along with information that is needed to perform the hybrid tiling.
 *
 * "bounds" are the bounds on the dependence distances that
 * define the hexagonal shape and the required skewing in the remaining
 * space dimensions.
 *
 * "input_node" points to the input pair of band nodes.
 * "input_schedule" is the partial schedule of this input pair of band nodes.
 * The space of this schedule is [P -> C], where P is the space
 * of the parent node and C is the space of the child node.
 *
 * "space_sizes" represent the total size of a tile for the space
 * dimensions, i.e., those corresponding to the child node.
 * The space of "space_sizes" is C.
 * If S_0 is the original tile size in the first space dimension,
 * then the first entry of "space_sizes" is equal to
 * W = 2*S_0 + floor(d_l h) + floor(d_u h).
 * The remaining entries are the same as in the original tile sizes.
 *
 * The basic hexagonal tiling "hex" is defined
 * in a "ts" (time-space) space and corresponds to the phase-1 tiles.
 * "time_tile" maps the "ts" space to outer time tile.
 * Is is equal to ts[t, s] -> floor(t/(2 * S_t)), with S_t the original tile
 * size corresponding to the parent node.
 * "local_time" maps the "ts" space to the time dimension inside each tile.
 * It is equal to ts[t, s] -> t mod (2 S_t), with S_t the original tile
 * size corresponding to the parent node.
 * "shift_space" shifts the tiles at time tile T = floor(t/(2 S_t))
 * in the space dimension such that they align to a multiple of W.
 * It is equal to ts[t, s] -> s + (-(2 * shift_s)*T) % W,
 * with shift_s = S_0 + floor(d_u h).
 * "shift_phase" is the shift taken to go from phase 0 to phase 1
 * It is equal to ts[t, s] -> ts[t + S_t, s + shift_s],
 * with shift_s = S_0 + floor(d_u h).
 *
 * "project_ts" projects the space of the input schedule to the ts-space.
 * It is equal to [P[t] -> C[s_0, ...]] -> ts[t, s_0].
 */
struct ppcg_ht_tiling {
        int ref;

        ppcg_ht_bounds *bounds;
        isl_schedule_node *input_node;
        isl_multi_union_pw_aff *input_schedule;

        isl_multi_val *space_sizes;

        isl_aff *time_tile;
        isl_aff *local_time;
        isl_aff *shift_space;
        isl_multi_aff *shift_phase;
        isl_set *hex;

        isl_multi_aff *project_ts;
};
typedef struct ppcg_ht_tiling ppcg_ht_tiling;

/* Return the space of the pair of band nodes that form the input
 * to the hybrid tiling.
 * In particular, return the space [P -> C], where P is the space
 * of the parent node and C is the space of the child node.
 */
__isl_give isl_space *ppcg_ht_tiling_get_input_space(
        __isl_keep ppcg_ht_tiling *tile)
{
        if (!tile)
                return NULL;

        return isl_multi_union_pw_aff_get_space(tile->input_schedule);
}

/* Remove a reference to "tile" and free "tile" along with all its fields
 * as soon as the reference count drops to zero.
 */
static __isl_null ppcg_ht_tiling *ppcg_ht_tiling_free(
        __isl_take ppcg_ht_tiling *tiling)
{
        if (!tiling)
                return NULL;
        if (--tiling->ref > 0)
                return NULL;

        ppcg_ht_bounds_free(tiling->bounds);
        isl_schedule_node_free(tiling->input_node);
        isl_multi_union_pw_aff_free(tiling->input_schedule);
        isl_multi_val_free(tiling->space_sizes);
        isl_aff_free(tiling->time_tile);
        isl_aff_free(tiling->local_time);
        isl_aff_free(tiling->shift_space);
        isl_multi_aff_free(tiling->shift_phase);
        isl_set_free(tiling->hex);
        isl_multi_aff_free(tiling->project_ts);
        free(tiling);

        return NULL;
}

/* Return a new reference to "tiling".
 */
__isl_give ppcg_ht_tiling *ppcg_ht_tiling_copy(
        __isl_keep ppcg_ht_tiling *tiling)
{
        if (!tiling)
                return NULL;

        tiling->ref++;
        return tiling;
}

/* Return the isl_ctx to which "tiling" belongs.
 */
isl_ctx *ppcg_ht_tiling_get_ctx(__isl_keep ppcg_ht_tiling *tiling)
{
        if (!tiling)
                return NULL;

        return isl_multi_union_pw_aff_get_ctx(tiling->input_schedule);
}

/* Representation of one of the two phases of hybrid tiling.
 *
 * "tiling" points to the shared tiling data.
 *
 * "time_tile", "local_time" and "shift_space" are equal to the corresponding
 * fields of "tiling", pulled back to the input space.
 * In case of phase 0, these expressions have also been moved
 * from phase 1 to phase 0.
 *
 * "domain" contains the hexagonal tiling of this phase.
 *
 * "space_shift" is the shift that should be added to the space band
 * in order to be able to apply rectangular tiling to the space.
 * For phase 1, it is equal to
 *
 *      [P[t] -> C[s_0, s_i]] -> C[(-(2 * shift_s)*T) % W, dl_i * u]
 *
 * with shift_s = S_0 + floor(d_u h),
 * T equal to "time_tile" and u equal to "local_time".
 * For phase 0, it is equal to
 *
 *      [P[t] -> C[s_0, s_i]] -> C[shift_s + (-(2 * shift_s)*T) % W, dl_i * u]
 *
 * "space_tile" is the space tiling.  It is equal to
 *
 *      [P[t] -> C[s]] -> C[floor((s + space_shift)/space_size]
 */
struct ppcg_ht_phase {
        ppcg_ht_tiling *tiling;

        isl_aff *time_tile;
        isl_aff *local_time;
        isl_aff *shift_space;
        isl_set *domain;

        isl_multi_aff *space_shift;
        isl_multi_aff *space_tile;
};

/* Free "phase" along with all its fields.
 */
static __isl_null ppcg_ht_phase *ppcg_ht_phase_free(
        __isl_take ppcg_ht_phase *phase)
{
        if (!phase)
                return NULL;

        ppcg_ht_tiling_free(phase->tiling);
        isl_aff_free(phase->time_tile);
        isl_aff_free(phase->local_time);
        isl_aff_free(phase->shift_space);
        isl_set_free(phase->domain);
        isl_multi_aff_free(phase->space_shift);
        isl_multi_aff_free(phase->space_tile);
        free(phase);

        return NULL;
}

/* Wrapper around ppcg_ht_phase_free for use as an argument
 * to isl_id_set_free_user.
 */
static void ppcg_ht_phase_free_wrap(void *user)
{
        ppcg_ht_phase *phase = user;

        ppcg_ht_phase_free(phase);
}

/* Return the domain of hybrid tiling phase "phase".
 */
static __isl_give isl_set *ppcg_ht_phase_get_domain(ppcg_ht_phase *phase)
{
        if (!phase)
                return NULL;

        return isl_set_copy(phase->domain);
}

/* Return the space of the pair of band nodes that form the input
 * to the hybrid tiling of which "phase" is a phase.
 * In particular, return the space [P -> C], where P is the space
 * of the parent node and C is the space of the child node.
 */
static __isl_give isl_space *ppcg_ht_phase_get_input_space(
        __isl_keep ppcg_ht_phase *phase)
{
        if (!phase)
                return NULL;

        return ppcg_ht_tiling_get_input_space(phase->tiling);
}

/* Construct the lower left constraint of the hexagonal tile, i.e.,
 *
 *      du a - b <= (2h+1) du - duh
 *      -du a + b + (2h+1) du - duh >= 0
 *
 * where duh = floor(du * h).
 *
 * This constraint corresponds to (6) in
 * "Hybrid Hexagonal/Classical Tiling for GPUs".
 */
static __isl_give isl_constraint *hex_lower_left(__isl_take isl_local_space *ls,
        __isl_keep isl_val *h, __isl_keep isl_val *du, __isl_keep isl_val *duh)
{
        isl_val *v;
        isl_aff *aff;

        v = isl_val_add_ui(isl_val_mul_ui(isl_val_copy(h), 2), 1);
        v = isl_val_mul(v, isl_val_copy(du));
        v = isl_val_sub(v, isl_val_copy(duh));
        aff = isl_aff_val_on_domain(ls, v);
        v = isl_val_neg(isl_val_copy(du));
        aff = isl_aff_set_coefficient_val(aff, isl_dim_in, 0, v);
        aff = isl_aff_set_coefficient_si(aff, isl_dim_in, 1, 1);

        return isl_inequality_from_aff(aff);
}

/* Construct the lower constraint of the hexagonal tile, i.e.,
 *
 *      a <= 2h+1
 *      -a + 2h+1 >= 0
 *
 * This constraint corresponds to (7) in
 * "Hybrid Hexagonal/Classical Tiling for GPUs".
 */
static __isl_give isl_constraint *hex_lower(__isl_take isl_local_space *ls,
        __isl_keep isl_val *h)
{
        isl_val *v;
        isl_aff *aff;

        v = isl_val_add_ui(isl_val_mul_ui(isl_val_copy(h), 2), 1);
        aff = isl_aff_val_on_domain(ls, v);
        aff = isl_aff_set_coefficient_si(aff, isl_dim_in, 0, -1);

        return isl_inequality_from_aff(aff);
}

/* Construct the lower right constraint of the hexagonal tile, i.e.,
 *
 *      dl a + b <= (2h+1) dl + duh + (s0-1)
 *      -dl a - b + (2h+1) dl + duh + (s0-1) >= 0
 *
 * where duh = floor(du * h).
 *
 * This constraint corresponds to (8) in
 * "Hybrid Hexagonal/Classical Tiling for GPUs".
 */
static __isl_give isl_constraint *hex_lower_right(
        __isl_take isl_local_space *ls, __isl_keep isl_val *h,
        __isl_keep isl_val *s0, __isl_keep isl_val *dl, __isl_keep isl_val *duh)
{
        isl_val *v;
        isl_aff *aff;

        v = isl_val_add_ui(isl_val_mul_ui(isl_val_copy(h), 2), 1);
        v = isl_val_mul(v, isl_val_copy(dl));
        v = isl_val_add(v, isl_val_copy(duh));
        v = isl_val_add(v, isl_val_copy(s0));
        v = isl_val_sub_ui(v, 1);
        aff = isl_aff_val_on_domain(ls, v);
        v = isl_val_neg(isl_val_copy(dl));
        aff = isl_aff_set_coefficient_val(aff, isl_dim_in, 0, v);
        aff = isl_aff_set_coefficient_si(aff, isl_dim_in, 1, -1);

        return isl_inequality_from_aff(aff);
}

/* Construct the upper left constraint of the hexagonal tile, i.e.,
 *
 *      dl a + b >= h dl - (d - 1)/d                            with d = den(dl)
 *      dl a + b - h dl + (d - 1)/d >= 0
 *
 * This constraint corresponds to (10) in
 * "Hybrid Hexagonal/Classical Tiling for GPUs".
 */
static __isl_give isl_constraint *hex_upper_left(__isl_take isl_local_space *ls,
        __isl_keep isl_val *h, __isl_keep isl_val *dl)
{
        isl_val *v, *d;
        isl_aff *aff;

        d = isl_val_get_den_val(dl);
        v = isl_val_sub_ui(isl_val_copy(d), 1);
        v = isl_val_div(v, d);
        v = isl_val_sub(v, isl_val_mul(isl_val_copy(h), isl_val_copy(dl)));
        aff = isl_aff_val_on_domain(ls, v);
        aff = isl_aff_set_coefficient_val(aff, isl_dim_in, 0, isl_val_copy(dl));
        aff = isl_aff_set_coefficient_si(aff, isl_dim_in, 1, 1);

        return isl_inequality_from_aff(aff);
}

/* Construct the upper right constraint of the hexagonal tile, i.e.,
 *
 *      du a - b >= du h - duh - (s0-1) - dlh - (d - 1)/d       with d = den(du)
 *      du a - b - du h + duh + (s0-1) + dlh + (d - 1)/d >= 0
 *
 * where dlh = floor(dl * h) and duh = floor(du * h).
 *
 * This constraint corresponds to (12) in
 * "Hybrid Hexagonal/Classical Tiling for GPUs".
 */
static __isl_give isl_constraint *hex_upper_right(
        __isl_take isl_local_space *ls, __isl_keep isl_val *h,
        __isl_keep isl_val *s0, __isl_keep isl_val *du,
        __isl_keep isl_val *dlh, __isl_keep isl_val *duh)
{
        isl_val *v, *d;
        isl_aff *aff;

        d = isl_val_get_den_val(du);
        v = isl_val_sub_ui(isl_val_copy(d), 1);
        v = isl_val_div(v, d);
        v = isl_val_sub(v, isl_val_mul(isl_val_copy(h), isl_val_copy(du)));
        v = isl_val_add(v, isl_val_copy(duh));
        v = isl_val_add(v, isl_val_copy(dlh));
        v = isl_val_add(v, isl_val_copy(s0));
        v = isl_val_sub_ui(v, 1);
        aff = isl_aff_val_on_domain(ls, v);
        aff = isl_aff_set_coefficient_val(aff, isl_dim_in, 0, isl_val_copy(du));
        aff = isl_aff_set_coefficient_si(aff, isl_dim_in, 1, -1);

        return isl_inequality_from_aff(aff);
}

/* Construct the uppper constraint of the hexagonal tile, i.e.,
 *
 *      a >= 0
 *
 * This constraint corresponds to (13) in
 * "Hybrid Hexagonal/Classical Tiling for GPUs".
 */
static __isl_give isl_constraint *hex_upper(__isl_take isl_local_space *ls)
{
        isl_aff *aff;

        aff = isl_aff_var_on_domain(ls, isl_dim_set, 0);

        return isl_inequality_from_aff(aff);
}

/* Construct the basic hexagonal tile shape.
 * "space" is the 2D space in which the hexagon should be constructed.
 * h is st-1, with st the tile size in the time dimension
 * s0 is the tile size in the space dimension
 * dl is a bound on the negative relative dependence distances, i.e.,
 *
 *      d_s >= -dl d_t
 *
 * du is a bound on the positive relative dependence distances, i.e.,
 *
 *      d_s <= du d_t
 *
 * with (d_t,d_s) any dependence distance vector.
 * dlh = floor(dl * h)
 * duh = floor(du * h)
 *
 * The shape of the hexagon is as follows:
 *
 *              0 dlh   dlh+s0-1
 *                 ______                __
 * 0              /      \_             /
 *               /         \_          /
 * h            /            \ ______ /
 * h+1          \_           //      \\_
 *                \_        //         \\_
 * 2h+1             \______//            \\
 *              0   duh   duh+s0-1
 *                           duh+s0-1+dlh
 *                                duh+s0-1+dlh+1+s0+1
 *
 * The next hexagon is shifted by duh + dlh + 2 * s0.
 *
 * The slope of the "/" constraints is dl.
 * The slope of the "\_" constraints is du.
 */
static __isl_give isl_set *compute_hexagon(__isl_take isl_space *space,
        __isl_keep isl_val *h, __isl_keep isl_val *s0,
        __isl_keep isl_val *dl, __isl_keep isl_val *du,
        __isl_keep isl_val *dlh, __isl_keep isl_val *duh)
{
        isl_local_space *ls;
        isl_constraint *c;
        isl_basic_set *bset;

        ls = isl_local_space_from_space(space);

        c = hex_lower_left(isl_local_space_copy(ls), h, du, duh);
        bset = isl_basic_set_from_constraint(c);

        c = hex_lower(isl_local_space_copy(ls), h);
        bset = isl_basic_set_add_constraint(bset, c);

        c = hex_lower_right(isl_local_space_copy(ls), h, s0, dl, duh);
        bset = isl_basic_set_add_constraint(bset, c);

        c = hex_upper_left(isl_local_space_copy(ls), h, dl);
        bset = isl_basic_set_add_constraint(bset, c);

        c = hex_upper_right(isl_local_space_copy(ls), h, s0, du, dlh, duh);
        bset = isl_basic_set_add_constraint(bset, c);

        c = hex_upper(ls);
        bset = isl_basic_set_add_constraint(bset, c);

        return isl_set_from_basic_set(bset);
}

/* Name of the ts-space.
 */
static const char *ts_space_name = "ts";

/* Construct and return the space ts[t, s].
 */
static __isl_give isl_space *construct_ts_space(isl_ctx *ctx)
{
        isl_space *s;

        s = isl_space_set_alloc(ctx, 0, 2);
        s = isl_space_set_tuple_name(s, isl_dim_set, ts_space_name);

        return s;
}

/* Name of the local ts-space.
 */
static const char *local_ts_space_name = "local_ts";

/* Construct and return the space local_ts[t, s].
 */
static __isl_give isl_space *construct_local_ts_space(isl_ctx *ctx)
{
        isl_space *s;

        s = isl_space_set_alloc(ctx, 0, 2);
        s = isl_space_set_tuple_name(s, isl_dim_set, local_ts_space_name);

        return s;
}

/* Compute the total size of a tile for the space dimensions,
 * i.e., those corresponding to the child node
 * of the input pattern.
 * If S_0 is the original tile size in the first space dimension,
 * then the first entry of "space_sizes" is equal to
 * W = 2*S_0 + floor(d_l h) + floor(d_u h).
 * The remaining entries are the same as in the original tile sizes.
 * "tile_sizes" contains the original tile sizes, including
 * the tile size corresponding to the parent node.
 * "dlh" is equal to floor(d_l h).
 * "duh" is equal to floor(d_u h).
 */
static __isl_give isl_multi_val *compute_space_sizes(
        __isl_keep isl_multi_val *tile_sizes,
        __isl_keep isl_val *dlh, __isl_keep isl_val *duh)
{
        isl_val *size;
        isl_multi_val *space_sizes;

        space_sizes = isl_multi_val_copy(tile_sizes);
        space_sizes = isl_multi_val_factor_range(space_sizes);
        size = isl_multi_val_get_val(space_sizes, 0);
        size = isl_val_mul_ui(size, 2);
        size = isl_val_add(size, isl_val_copy(duh));
        size = isl_val_add(size, isl_val_copy(dlh));
        space_sizes = isl_multi_val_set_val(space_sizes, 0, size);

        return space_sizes;
}

/* Compute the offset of phase 1 with respect to phase 0
 * in the ts-space ("space").
 * In particular, return
 *
 *      ts[st, s0 + duh]
 */
static __isl_give isl_multi_val *compute_phase_shift(
        __isl_keep isl_space *space, __isl_keep isl_val *st,
        __isl_keep isl_val *s0, __isl_keep isl_val *duh)
{
        isl_val *v;
        isl_multi_val *phase_shift;

        phase_shift = isl_multi_val_zero(isl_space_copy(space));
        phase_shift = isl_multi_val_set_val(phase_shift, 0, isl_val_copy(st));
        v = isl_val_add(isl_val_copy(duh), isl_val_copy(s0));
        phase_shift = isl_multi_val_set_val(phase_shift, 1, v);

        return phase_shift;
}

/* Return the function
 *
 *      ts[t, s] -> floor(t/(2 * st))
 *
 * representing the time tile.
 * "space" is the space ts[t, s].
 */
static __isl_give isl_aff *compute_time_tile(__isl_keep isl_space *space,
        __isl_keep isl_val *st)
{
        isl_val *v;
        isl_aff *t;
        isl_local_space *ls;

        ls = isl_local_space_from_space(isl_space_copy(space));
        t = isl_aff_var_on_domain(ls, isl_dim_set, 0);
        v = isl_val_mul_ui(isl_val_copy(st), 2);
        t = isl_aff_floor(isl_aff_scale_down_val(t, v));

        return t;
}

/* Compute a shift in the space dimension for tiles
 * at time tile T = floor(t/(2 * S_t))
 * such that they align to a multiple of the total space tile dimension W.
 * In particular, compute
 *
 *      ts[t, s] -> s + (-(2 * shift_s)*T) % W
 *
 * where shift_s is the shift of phase 1 with respect to phase 0
 * in the space dimension (the first element of "phase_shift").
 * W is stored in the first element of "space_sizes".
 * "time_tile" is the function
 *
 *      ts[t, s] -> floor(t/(2 * S_T))
 *
 * Since phase 1 is shifted by shift_s with respect to phase 0,
 * the next line of phase 0 (at T+1) is shifted by 2*shift_s
 * with respect to the previous line (at T).
 * A shift of -(2 * shift_s)*T therefore allows the basic pattern
 * (which starts at 0) to be applied.
 * However, this shift will be used to obtain the tile coordinate
 * in the first space dimension and if the original values
 * in the space dimension are non-negative, then the shift should
 * not make them negative.  Moreover, the shift should be as minimal
 * as possible.
 * Since the pattern repeats itself with a period of W in the space
 * dimension, the shift can be replaced by (-(2 * shift_s)*T) % W.
 */
static __isl_give isl_aff *compute_shift_space(__isl_keep isl_aff *time_tile,
        __isl_keep isl_multi_val *space_sizes,
        __isl_keep isl_multi_val *phase_shift)
{
        isl_val *v;
        isl_aff *s, *t;
        isl_local_space *ls;

        ls = isl_local_space_from_space(isl_aff_get_domain_space(time_tile));
        t = isl_aff_copy(time_tile);
        v = isl_val_mul_ui(isl_multi_val_get_val(phase_shift, 1), 2);
        v = isl_val_neg(v);
        t = isl_aff_scale_val(t, v);
        v = isl_multi_val_get_val(space_sizes, 0);
        t = isl_aff_mod_val(t, v);
        s = isl_aff_var_on_domain(ls, isl_dim_set, 1);
        s = isl_aff_add(s, t);

        return s;
}

/* Give the phase_shift ts[S_t, S_0 + floor(d_u h)],
 * compute a function that applies the shift, i.e.,
 *
 *      ts[t, s] -> ts[t + S_t, s + S_0 + floor(d_u h)],
 */
static __isl_give isl_multi_aff *compute_shift_phase(
        __isl_keep isl_multi_val *phase_shift)
{
        isl_space *space;
        isl_multi_aff *shift;

        space = isl_multi_val_get_space(phase_shift);
        shift = isl_multi_aff_multi_val_on_space(space,
                                        isl_multi_val_copy(phase_shift));
        space = isl_multi_aff_get_space(shift);
        shift = isl_multi_aff_add(shift, isl_multi_aff_identity(space));

        return shift;
}

/* Compute a mapping from the ts-space to the local coordinates
 * within each tile.  In particular, compute
 *
 *      ts[t, s] -> local_ts[t % (2 S_t), (s + (-(2 * shift_s)*T) % W) % W]
 *
 * "ts" is the space ts[t, s]
 * "local_ts" is the space local_ts[t, s]
 * "shift_space" is equal to ts[t, s] -> s + (-(2 * shift_s)*T) % W
 * "st" is the tile size in the time dimension S_t.
 * The first element of "space_sizes" is equal to W.
 */
static __isl_give isl_multi_aff *compute_localize(
        __isl_keep isl_space *local_ts, __isl_keep isl_aff *shift_space,
        __isl_keep isl_val *st, __isl_keep isl_multi_val *space_sizes)
{
        isl_val *v;
        isl_space *space;
        isl_aff *s, *t;
        isl_multi_aff *localize;

        space = isl_aff_get_domain_space(shift_space);
        local_ts = isl_space_copy(local_ts);
        space = isl_space_map_from_domain_and_range(space, local_ts);
        localize = isl_multi_aff_identity(space);
        t = isl_multi_aff_get_aff(localize, 0);
        v = isl_val_mul_ui(isl_val_copy(st), 2);
        t = isl_aff_mod_val(t, v);
        localize = isl_multi_aff_set_aff(localize, 0, t);
        s = isl_aff_copy(shift_space);
        v = isl_multi_val_get_val(space_sizes, 0);
        s = isl_aff_mod_val(s, v);
        localize = isl_multi_aff_set_aff(localize, 1, s);

        return localize;
}

/* Set the project_ts field of "tiling".
 *
 * This field projects the space of the input schedule to the ts-space.
 * It is equal to [P[t] -> C[s_0, ...]] -> ts[t, s_0].
 */
static __isl_give ppcg_ht_tiling *ppcg_ht_tiling_set_project_ts(
        __isl_take ppcg_ht_tiling *tiling)
{
        int n;
        isl_space *space;
        isl_multi_aff *project;

        if (!tiling)
                return NULL;

        space = ppcg_ht_tiling_get_input_space(tiling);
        n = isl_space_dim(space, isl_dim_set);
        project = isl_multi_aff_project_out_map(space, isl_dim_set, 2, n - 2);
        project = isl_multi_aff_set_tuple_name(project,
                                                isl_dim_out, ts_space_name);
        if (!project)
                return ppcg_ht_tiling_free(tiling);

        tiling->project_ts = project;

        return tiling;
}

/* Construct a hybrid tiling description from bounds on the dependence
 * distances "bounds".
 * "input_node" points to the original parent node.
 * "input_schedule" is the combined schedule of the parent and child
 * node in the input.
 * "tile_sizes" are the original, user specified tile sizes.
 */
static __isl_give ppcg_ht_tiling *ppcg_ht_bounds_construct_tiling(
        __isl_take ppcg_ht_bounds *bounds,
        __isl_keep isl_schedule_node *input_node,
        __isl_keep isl_multi_union_pw_aff *input_schedule,
        __isl_keep isl_multi_val *tile_sizes)
{
        isl_ctx *ctx;
        ppcg_ht_tiling *tiling;
        isl_multi_val *space_sizes, *phase_shift;
        isl_aff *time_tile, *shift_space;
        isl_multi_aff *localize;
        isl_val *h, *duh, *dlh;
        isl_val *st, *s0, *du, *dl;
        isl_space *ts, *local_ts;

        if (!bounds || !input_node || !input_schedule || !tile_sizes)
                goto error;

        ctx = isl_multi_union_pw_aff_get_ctx(input_schedule);
        tiling = isl_calloc_type(ctx, struct ppcg_ht_tiling);
        if (!tiling)
                goto error;
        tiling->ref = 1;

        st = isl_multi_val_get_val(tile_sizes, 0);
        h = isl_val_sub_ui(isl_val_copy(st), 1);
        s0 = isl_multi_val_get_val(tile_sizes, 1);
        du = ppcg_ht_bounds_get_upper(bounds);
        dl = ppcg_ht_bounds_get_lower(bounds, 0);

        duh = isl_val_floor(isl_val_mul(isl_val_copy(du), isl_val_copy(h)));
        dlh = isl_val_floor(isl_val_mul(isl_val_copy(dl), isl_val_copy(h)));

        ts = construct_ts_space(ctx);
        local_ts = construct_local_ts_space(ctx);

        space_sizes = compute_space_sizes(tile_sizes, dlh, duh);
        phase_shift = compute_phase_shift(ts, st, s0, duh);
        time_tile = compute_time_tile(ts, st);
        shift_space = compute_shift_space(time_tile, space_sizes, phase_shift);
        localize = compute_localize(local_ts, shift_space, st, space_sizes);
        isl_space_free(ts);

        tiling->input_node = isl_schedule_node_copy(input_node);
        tiling->input_schedule = isl_multi_union_pw_aff_copy(input_schedule);
        tiling->space_sizes = space_sizes;
        tiling->bounds = bounds;
        tiling->local_time = isl_multi_aff_get_aff(localize, 0);
        tiling->hex = compute_hexagon(local_ts, h, s0, dl, du, dlh, duh);
        tiling->hex = isl_set_preimage_multi_aff(tiling->hex, localize);
        tiling->time_tile = time_tile;
        tiling->shift_space = shift_space;
        tiling->shift_phase = compute_shift_phase(phase_shift);
        isl_multi_val_free(phase_shift);

        isl_val_free(duh);
        isl_val_free(dlh);
        isl_val_free(du);
        isl_val_free(dl);
        isl_val_free(s0);
        isl_val_free(st);
        isl_val_free(h);

        if (!tiling->input_schedule || !tiling->local_time || !tiling->hex ||
            !tiling->shift_space || !tiling->shift_phase)
                return ppcg_ht_tiling_free(tiling);

        tiling = ppcg_ht_tiling_set_project_ts(tiling);

        return tiling;
error:
        ppcg_ht_bounds_free(bounds);
        return NULL;
}

/* Are all members of the band node "node" coincident?
 */
static isl_bool all_coincident(__isl_keep isl_schedule_node *node)
{
        int i, n;

        n = isl_schedule_node_band_n_member(node);
        for (i = 0; i < n; ++i) {
                isl_bool c;

                c = isl_schedule_node_band_member_get_coincident(node, i);
                if (c < 0 || !c)
                        return c;
        }

        return isl_bool_true;
}

/* Does "node" satisfy the properties of the inner node in the input
 * pattern for hybrid tiling?
 * That is, is it a band node with only coincident members, of which
 * there is at least one?
 */
static isl_bool has_child_properties(__isl_keep isl_schedule_node *node)
{
        if (!node)
                return isl_bool_error;
        if (isl_schedule_node_get_type(node) != isl_schedule_node_band)
                return isl_bool_false;
        if (isl_schedule_node_band_n_member(node) < 1)
                return isl_bool_false;
        return all_coincident(node);
}

/* Does "node" satisfy the properties of the outer node in the input
 * pattern for hybrid tiling?
 * That is, is it a band node with a single member?
 */
static isl_bool has_parent_properties(__isl_keep isl_schedule_node *node)
{
        if (!node)
                return isl_bool_error;
        if (isl_schedule_node_get_type(node) != isl_schedule_node_band)
                return isl_bool_false;
        if (isl_schedule_node_band_n_member(node) != 1)
                return isl_bool_false;
        return isl_bool_true;
}

/* Does the parent of "node" satisfy the input patttern for hybrid tiling?
 * That is, does "node" satisfy the properties of the inner node and
 * does the parent of "node" satisfy the properties of the outer node?
 */
isl_bool ppcg_ht_parent_has_input_pattern(__isl_keep isl_schedule_node *node)
{
        isl_bool has_pattern;

        has_pattern = has_child_properties(node);
        if (has_pattern < 0 || !has_pattern)
                return has_pattern;

        node = isl_schedule_node_copy(node);
        node = isl_schedule_node_parent(node);
        has_pattern = has_parent_properties(node);
        isl_schedule_node_free(node);

        return has_pattern;
}

/* Does "node" satisfy the input patttern for hybrid tiling?
 * That is, does "node" satisfy the properties of the outer node and
 * does the child of "node" satisfy the properties of the inner node?
 */
isl_bool ppcg_ht_has_input_pattern(__isl_keep isl_schedule_node *node)
{
        isl_bool has_pattern;

        has_pattern = has_parent_properties(node);
        if (has_pattern < 0 || !has_pattern)
                return has_pattern;

        node = isl_schedule_node_get_child(node, 0);
        has_pattern = has_child_properties(node);
        isl_schedule_node_free(node);

        return has_pattern;
}

/* Check that "node" satisfies the input pattern for hybrid tiling.
 * Error out if it does not.
 */
static isl_stat check_input_pattern(__isl_keep isl_schedule_node *node)
{
        isl_bool has_pattern;

        has_pattern = ppcg_ht_has_input_pattern(node);
        if (has_pattern < 0)
                return isl_stat_error;
        if (!has_pattern)
                isl_die(isl_schedule_node_get_ctx(node), isl_error_invalid,
                        "invalid input pattern for hybrid tiling",
                        return isl_stat_error);

        return isl_stat_ok;
}

/* Extract the input schedule from "node", i.e., the product
 * of the partial schedules of the parent and child nodes
 * in the input pattern.
 */
static __isl_give isl_multi_union_pw_aff *extract_input_schedule(
        __isl_keep isl_schedule_node *node)
{
        isl_multi_union_pw_aff *partial, *partial2;

        partial = isl_schedule_node_band_get_partial_schedule(node);
        node = isl_schedule_node_get_child(node, 0);
        partial2 = isl_schedule_node_band_get_partial_schedule(node);
        isl_schedule_node_free(node);

        return isl_multi_union_pw_aff_range_product(partial, partial2);
}

/* Collect all dependences from "scop" that are relevant for performing
 * hybrid tiling on "node" and its child and map them to the schedule
 * space of this pair of nodes.
 *
 * In case live range reordering is not used,
 * the flow and the false dependences are collected.
 * In case live range reordering is used,
 * the flow and the forced dependences are collected, as well
 * as the order dependences that are adjacent to non-local
 * flow dependences.
 *
 * In all cases, only dependences that map to the same instance
 * of the outer part of the schedule are considered.
 */
static __isl_give isl_map *collect_deps(struct ppcg_scop *scop,
        __isl_keep isl_schedule_node *node)
{
        isl_space *space;
        isl_multi_union_pw_aff *prefix, *partial;
        isl_union_map *flow, *other, *dep, *umap;
        isl_map *map;

        prefix = isl_schedule_node_get_prefix_schedule_multi_union_pw_aff(node);
        partial = extract_input_schedule(node);
        space = isl_multi_union_pw_aff_get_space(partial);

        flow = isl_union_map_copy(scop->dep_flow);
        flow = isl_union_map_eq_at_multi_union_pw_aff(flow,
                                        isl_multi_union_pw_aff_copy(prefix));
        if (!scop->options->live_range_reordering) {
                other = isl_union_map_copy(scop->dep_false);
                other = isl_union_map_eq_at_multi_union_pw_aff(other, prefix);
        } else {
                isl_union_map *local, *non_local, *order, *adj;
                isl_union_set *domain, *range;

                other = isl_union_map_copy(scop->dep_forced);
                other = isl_union_map_eq_at_multi_union_pw_aff(other,
                                        isl_multi_union_pw_aff_copy(prefix));
                local = isl_union_map_copy(flow);
                local = isl_union_map_eq_at_multi_union_pw_aff(local,
                                        isl_multi_union_pw_aff_copy(partial));
                non_local = isl_union_map_copy(flow);
                non_local = isl_union_map_subtract(non_local, local);

                order = isl_union_map_copy(scop->dep_order);
                order = isl_union_map_eq_at_multi_union_pw_aff(order, prefix);
                adj = isl_union_map_copy(order);
                domain = isl_union_map_domain(isl_union_map_copy(non_local));
                domain = isl_union_set_coalesce(domain);
                adj = isl_union_map_intersect_range(adj, domain);
                other = isl_union_map_union(other, adj);

                adj = order;
                range = isl_union_map_range(non_local);
                range = isl_union_set_coalesce(range);
                adj = isl_union_map_intersect_domain(adj, range);
                other = isl_union_map_union(other, adj);
        }
        dep = isl_union_map_union(flow, other);

        umap = isl_union_map_from_multi_union_pw_aff(partial);
        dep = isl_union_map_apply_domain(dep, isl_union_map_copy(umap));
        dep = isl_union_map_apply_range(dep, umap);

        space = isl_space_map_from_set(space);
        map = isl_union_map_extract_map(dep, space);
        isl_union_map_free(dep);

        map = isl_map_coalesce(map);

        return map;
}

/* Given a constraint of the form
 *
 *      a i_0 + b i_1 >= 0
 * or
 *      a i_0 + b i_1 = 0
 *
 * use it to update one or both of the non-negative bounds
 * in "list" = (min, max) such that
 *
 *      i_1 >= -min i_0
 * and
 *      i_1 <= max i_0
 *
 * If b = 0, then the constraint cannot be used.
 * Otherwise, the constraint is equivalent to
 *
 *      sgn(b) i_1 >= - a/abs(b) i_0
 * i.e.,
 *      i_1 >= - a/abs(b) i_0
 * or
 *      i_1 <= a/abs(b) i_0
 *
 * Set the first or second element of "list" to max(0, a/abs(b)),
 * according to the sign of "b".  Or set both in case the constraint
 * is an equality, taking into account the sign change.
 */
static __isl_give isl_val_list *list_set_min_max(__isl_take isl_val_list *list,
        __isl_keep isl_constraint *c)
{
        isl_val *a, *b;
        int sign;
        int pos;
        isl_bool eq, is_zero, is_neg;

        eq = isl_constraint_is_equality(c);
        if (eq < 0)
                return isl_val_list_free(list);

        b = isl_constraint_get_coefficient_val(c, isl_dim_set, 1);
        is_zero = isl_val_is_zero(b);
        if (is_zero == isl_bool_true) {
                isl_val_free(b);
                return list;
        }
        a = isl_constraint_get_coefficient_val(c, isl_dim_set, 0);
        sign = isl_val_sgn(b);
        b = isl_val_abs(b);
        a = isl_val_div(a, b);

        if (eq)
                b = isl_val_copy(a);

        pos = sign > 0 ? 0 : 1;
        is_neg = isl_val_is_neg(a);
        if (is_neg == isl_bool_true)
                a = isl_val_set_si(a, 0);
        list = isl_val_list_set_val(list, pos, a);

        if (!eq)
                return is_neg < 0 ? isl_val_list_free(list) : list;

        pos = 1 - pos;
        a = isl_val_neg(b);
        is_neg = isl_val_is_neg(a);
        if (is_neg == isl_bool_true)
                a = isl_val_set_si(a, 0);
        list = isl_val_list_set_val(list, pos, a);

        return is_neg < 0 ? isl_val_list_free(list) : list;
}

/* If constraint "c" passes through the origin, then try and use it
 * to update the non-negative bounds in "list" = (min, max) such that
 *
 *      i_1 >= -min i_0
 * and
 *      i_1 <= max i_0
 */
static isl_stat set_min_max(__isl_take isl_constraint *c, void *user)
{
        isl_val *v;
        isl_val_list **list = user;
        isl_bool is_zero;

        v = isl_constraint_get_constant_val(c);
        is_zero = isl_val_is_zero(v);
        isl_val_free(v);

        if (is_zero == isl_bool_true)
                *list = list_set_min_max(*list, c);

        isl_constraint_free(c);
        return is_zero < 0 ? isl_stat_error : isl_stat_ok;
}

/* Given a set of dependence distance vectors "dist", compute
 * pair of non-negative bounds min and max such that
 *
 *      d_pos >= -min d_0
 * and
 *      d_pos <= max d_0
 *
 * and return the pair (min, max).
 * If no bound can be found in either direction, then the bound
 * is replaced by NaN.
 *
 * The dependence distances are first projected onto the (d_0, d_pos).
 * Then the zero dependence distance is added and the convex hull is computed.
 * Finally, the bounds are extracted from the constraints of the convex hull
 * that pass through the origin.
 */
static __isl_give isl_val_list *min_max_dist(__isl_keep isl_set *dist, int pos)
{
        isl_space *space;
        isl_basic_set *hull;
        int dim;
        isl_ctx *ctx;
        isl_val *nan;
        isl_val_list *list;

        ctx = isl_set_get_ctx(dist);
        nan = isl_val_nan(ctx);
        list = isl_val_list_alloc(ctx, 2);
        list = isl_val_list_add(list, isl_val_copy(nan));
        list = isl_val_list_add(list, nan);

        dist = isl_set_copy(dist);
        dim = isl_set_dim(dist, isl_dim_set);
        if (dist && pos >= dim)
                isl_die(ctx, isl_error_internal, "position out of bounds",
                        dist = isl_set_free(dist));
        dist = isl_set_project_out(dist, isl_dim_set, pos + 1, dim - (pos + 1));
        dist = isl_set_project_out(dist, isl_dim_set, 1, pos - 1);

        space = isl_set_get_space(dist);
        dist = isl_set_union(dist, isl_set_from_point(isl_point_zero(space)));
        dist = isl_set_remove_divs(dist);
        hull = isl_set_convex_hull(dist);

        if (isl_basic_set_foreach_constraint(hull, &set_min_max, &list) < 0)
                list = isl_val_list_free(list);
        isl_basic_set_free(hull);

        return list;
}

/* Given a schedule node "node" that, together with its child,
 * satisfies the input pattern for hybrid tiling, compute bounds
 * on the relative dependence distances of the child node with
 * respect to the parent node.  These bounds are needed to
 * construct a hybrid tiling.
 *
 * First all relevant dependences are collected and mapped
 * to the schedule space of the pair of nodes.  Then, the
 * dependence distances are computed in this space.
 *
 * These dependence distances are then projected onto a two-dimensional
 * space consisting of the single schedule dimension of the outer node
 * and one of the schedule dimensions of the inner node.
 * The maximal and minimal relative dependence distances are extracted
 * from these projections.
 * This process is repeated for each of the schedule dimensions
 * of the inner node.  For the first dimension, both minimal and
 * maximal relative dependence distances are stored in the result.
 * For the other dimensions, only the minimal relative dependence
 * distance is stored.
 */
__isl_give ppcg_ht_bounds *ppcg_ht_compute_bounds(struct ppcg_scop *scop,
        __isl_keep isl_schedule_node *node)
{
        ppcg_ht_bounds *bnd;
        isl_space *space;
        isl_map *map;
        isl_set *dist;
        isl_val_list *pair;
        isl_schedule_node *child;
        int n;
        int i, dim;

        if (!scop || !node || check_input_pattern(node) < 0)
                return NULL;

        child = isl_schedule_node_get_child(node, 0);
        space = isl_schedule_node_band_get_space(child);
        dim = isl_schedule_node_band_n_member(child);
        isl_schedule_node_free(child);
        bnd = ppcg_ht_bounds_alloc(space);
        if (!bnd)
                return NULL;

        map = collect_deps(scop, node);

        dist = isl_map_deltas(map);
        n = isl_set_dim(dist, isl_dim_param);
        dist = isl_set_project_out(dist, isl_dim_param, 0, n);

        pair = min_max_dist(dist, 1);
        bnd = ppcg_ht_bounds_set_lower(bnd, 0, isl_val_list_get_val(pair, 0));
        bnd = ppcg_ht_bounds_set_upper(bnd, isl_val_list_get_val(pair, 1));
        isl_val_list_free(pair);

        for (i = 1; i < dim; ++i) {
                pair = min_max_dist(dist, 1 + i);
                bnd = ppcg_ht_bounds_set_lower(bnd, i,
                                                isl_val_list_get_val(pair, 0));
                isl_val_list_free(pair);
        }

        isl_set_free(dist);

        return bnd;
}

/* Check if all the fields of "phase" are valid, freeing "phase"
 * if they are not.
 */
static __isl_give ppcg_ht_phase *check_phase(__isl_take ppcg_ht_phase *phase)
{
        if (!phase)
                return NULL;

        if (!phase->tiling || !phase->local_time ||
            !phase->shift_space || !phase->domain)
                return ppcg_ht_phase_free(phase);

        return phase;
}

/* Construct a ppcg_ht_phase object, that simply copies
 * information from "tiling".
 * That is, the result is defined over the "ts" space and
 * corresponds to phase 1.
 */
static __isl_give ppcg_ht_phase *construct_phase(
        __isl_keep ppcg_ht_tiling *tiling)
{
        isl_ctx *ctx;
        ppcg_ht_phase *phase;

        if (!tiling)
                return NULL;

        ctx = ppcg_ht_tiling_get_ctx(tiling);
        phase = isl_calloc_type(ctx, struct ppcg_ht_phase);
        if (!phase)
                return NULL;
        phase->tiling = ppcg_ht_tiling_copy(tiling);
        phase->time_tile = isl_aff_copy(tiling->time_tile);
        phase->local_time = isl_aff_copy(tiling->local_time);
        phase->shift_space = isl_aff_copy(tiling->shift_space);
        phase->domain = isl_set_copy(tiling->hex);

        return check_phase(phase);
}

/* Align the parameters of the elements of "phase" to those of "space".
 */
static __isl_give ppcg_ht_phase *phase_align_params(
        __isl_take ppcg_ht_phase *phase, __isl_take isl_space *space)
{
        if (!phase)
                goto error;

        phase->time_tile = isl_aff_align_params(phase->time_tile,
                                                        isl_space_copy(space));
        phase->local_time = isl_aff_align_params(phase->local_time,
                                                        isl_space_copy(space));
        phase->shift_space = isl_aff_align_params(phase->shift_space,
                                                        isl_space_copy(space));
        phase->domain = isl_set_align_params(phase->domain, space);

        return check_phase(phase);
error:
        isl_space_free(space);
        return NULL;
}

/* Pull back "phase" over "ma".
 * That is, take a phase defined over the range of "ma" and
 * turn it into a phase defined over the domain of "ma".
 */
static __isl_give ppcg_ht_phase *pullback_phase(__isl_take ppcg_ht_phase *phase,
        __isl_take isl_multi_aff *ma)
{
        phase = phase_align_params(phase, isl_multi_aff_get_space(ma));
        if (!phase)
                goto error;

        phase->time_tile = isl_aff_pullback_multi_aff(phase->time_tile,
                                                        isl_multi_aff_copy(ma));
        phase->local_time = isl_aff_pullback_multi_aff(phase->local_time,
                                                        isl_multi_aff_copy(ma));
        phase->shift_space = isl_aff_pullback_multi_aff(phase->shift_space,
                                                        isl_multi_aff_copy(ma));
        phase->domain = isl_set_preimage_multi_aff(phase->domain, ma);

        return check_phase(phase);
error:
        isl_multi_aff_free(ma);
        return NULL;
}

/* Pullback "phase" over phase->tiling->shift_phase, which shifts
 * phase 0 to phase 1.  The pullback therefore takes a phase 1
 * description and turns it into a phase 0 description.
 */
static __isl_give ppcg_ht_phase *shift_phase(__isl_take ppcg_ht_phase *phase)
{
        ppcg_ht_tiling *tiling;

        if (!phase)
                return NULL;

        tiling = phase->tiling;
        return pullback_phase(phase, isl_multi_aff_copy(tiling->shift_phase));
}

/* Take a "phase" defined over the ts-space and plug in the projection
 * from the input schedule space to the ts-space.
 * The result is then defined over this input schedule space.
 */
static __isl_give ppcg_ht_phase *lift_phase(__isl_take ppcg_ht_phase *phase)
{
        ppcg_ht_tiling *tiling;

        if (!phase)
                return NULL;

        tiling = phase->tiling;
        return pullback_phase(phase, isl_multi_aff_copy(tiling->project_ts));
}

/* Compute the shift that should be added to the space band
 * in order to be able to apply rectangular tiling to the space.
 * Store the shift in phase->space_shift.
 *
 * In the first dimension, it is equal to shift_space - s.
 * For phase 1, this results in
 *
 *      (-(2 * shift_s)*T) % W
 *
 * In phase 0, the "s" in shift_space has been replaced by "s + shift_s",
 * so the result is
 *
 *      shift_s + (-(2 * shift_s)*T) % W
 *
 * In the other dimensions, the shift is equal to
 *
 *      dl_i * local_time.
 */
static __isl_give ppcg_ht_phase *compute_space_shift(
        __isl_take ppcg_ht_phase *phase)
{
        int i, n;
        isl_space *space;
        isl_local_space *ls;
        isl_aff *aff, *s;
        isl_multi_aff *space_shift;

        if (!phase)
                return NULL;

        space = ppcg_ht_phase_get_input_space(phase);
        space = isl_space_unwrap(space);
        space = isl_space_range_map(space);

        space_shift = isl_multi_aff_zero(space);
        aff = isl_aff_copy(phase->shift_space);
        ls = isl_local_space_from_space(isl_aff_get_domain_space(aff));
        s = isl_aff_var_on_domain(ls, isl_dim_set, 1);
        aff = isl_aff_sub(aff, s);
        space_shift = isl_multi_aff_set_aff(space_shift, 0, aff);

        n = isl_multi_aff_dim(space_shift, isl_dim_out);
        for (i = 1; i < n; ++i) {
                isl_val *v;
                isl_aff *time;

                v = ppcg_ht_bounds_get_lower(phase->tiling->bounds, i);
                time = isl_aff_copy(phase->local_time);
                time = isl_aff_scale_val(time, v);
                space_shift = isl_multi_aff_set_aff(space_shift, i, time);
        }

        if (!space_shift)
                return ppcg_ht_phase_free(phase);
        phase->space_shift = space_shift;
        return phase;
}

/* Compute the space tiling and store the result in phase->space_tile.
 * The space tiling is of the form
 *
 *      [P[t] -> C[s]] -> C[floor((s + space_shift)/space_size]
 */
static __isl_give ppcg_ht_phase *compute_space_tile(
        __isl_take ppcg_ht_phase *phase)
{
        isl_space *space;
        isl_multi_val *space_sizes;
        isl_multi_aff *space_shift;
        isl_multi_aff *tile;

        if (!phase)
                return NULL;

        space = ppcg_ht_phase_get_input_space(phase);
        space = isl_space_unwrap(space);
        tile = isl_multi_aff_range_map(space);
        space_shift = isl_multi_aff_copy(phase->space_shift);
        tile = isl_multi_aff_add(space_shift, tile);
        space_sizes = isl_multi_val_copy(phase->tiling->space_sizes);
        tile = isl_multi_aff_scale_down_multi_val(tile, space_sizes);
        tile = isl_multi_aff_floor(tile);

        if (!tile)
                return ppcg_ht_phase_free(phase);
        phase->space_tile = tile;
        return phase;
}

/* Construct a representation for one of the two phase for hybrid tiling
 * "tiling".  If "shift" is not set, then the phase is constructed
 * directly from the hexagonal tile shape in "tiling", which represents
 * the phase-1 tiles.  If "shift" is set, then this tile shape is shifted
 * back over tiling->shift_phase to obtain the phase-0 tiles.
 *
 * First copy data from "tiling", then optionally shift the phase and
 * finally move the tiling from the "ts" space of "tiling" to
 * the space of the input pattern.
 *
 * After the basic phase has been computed, also compute
 * the corresponding space shift.
 */
static __isl_give ppcg_ht_phase *ppcg_ht_tiling_compute_phase(
        __isl_keep ppcg_ht_tiling *tiling, int shift)
{
        ppcg_ht_phase *phase;

        phase = construct_phase(tiling);
        if (shift)
                phase = shift_phase(phase);
        phase = lift_phase(phase);

        phase = compute_space_shift(phase);
        phase = compute_space_tile(phase);

        return phase;
}

/* Consruct a function that is equal to the time tile of "phase0"
 * on the domain of "phase0" and equal to the time tile of "phase1"
 * on the domain of "phase1".
 * The two domains are assumed to form a partition of the input
 * schedule space.
 */
static __isl_give isl_pw_multi_aff *combine_time_tile(
        __isl_keep ppcg_ht_phase *phase0, __isl_keep ppcg_ht_phase *phase1)
{
        isl_aff *T;
        isl_pw_aff *time, *time1;

        if (!phase0 || !phase1)
                return NULL;

        T = isl_aff_copy(phase0->time_tile);
        time = isl_pw_aff_alloc(ppcg_ht_phase_get_domain(phase0), T);

        T = isl_aff_copy(phase1->time_tile);
        time1 = isl_pw_aff_alloc(ppcg_ht_phase_get_domain(phase1), T);

        time = isl_pw_aff_union_add(time, time1);

        return isl_pw_multi_aff_from_pw_aff(time);
}

/* Name used in mark nodes that contain a pointer to a ppcg_ht_phase.
 */
static char *ppcg_phase_name = "phase";

/* Does "id" contain a pointer to a ppcg_ht_phase?
 * That is, is it called "phase"?
 */
static isl_bool is_phase_id(__isl_keep isl_id *id)
{
        const char *name;

        name = isl_id_get_name(id);
        if (!name)
                return isl_bool_error;

        return !strcmp(name, ppcg_phase_name);
}

/* Given a mark node with an identifier that points to a ppcg_ht_phase,
 * extract this ppcg_ht_phase pointer.
 */
__isl_keep ppcg_ht_phase *ppcg_ht_phase_extract_from_mark(
        __isl_keep isl_schedule_node *node)
{
        isl_bool is_phase;
        isl_id *id;
        void *p;

        if (!node)
                return NULL;
        if (isl_schedule_node_get_type(node) != isl_schedule_node_mark)
                isl_die(isl_schedule_node_get_ctx(node), isl_error_internal,
                        "not a phase mark", return NULL);

        id = isl_schedule_node_mark_get_id(node);
        is_phase = is_phase_id(id);
        p = isl_id_get_user(id);
        isl_id_free(id);

        if (is_phase < 0)
                return NULL;
        if (!is_phase)
                isl_die(isl_schedule_node_get_ctx(node), isl_error_internal,
                        "not a phase mark", return NULL);

        return p;
}

/* Insert a mark node at "node" holding a pointer to "phase".
 */
static __isl_give isl_schedule_node *insert_phase(
        __isl_take isl_schedule_node *node, __isl_take ppcg_ht_phase *phase)
{
        isl_ctx *ctx;
        isl_id *id;

        if (!node)
                goto error;
        ctx = isl_schedule_node_get_ctx(node);
        id = isl_id_alloc(ctx, ppcg_phase_name, phase);
        if (!id)
                goto error;
        id = isl_id_set_free_user(id, &ppcg_ht_phase_free_wrap);
        node = isl_schedule_node_insert_mark(node, id);

        return node;
error:
        ppcg_ht_phase_free(phase);
        isl_schedule_node_free(node);
        return NULL;
}

/* Construct a mapping from the elements of the original pair of bands
 * to which tiling was applied that belong to a tile of "phase"
 * to that tile, preserving the values for the outer bands.
 *
 * The mapping is of the form
 *
 *      [[outer] -> [P -> C]] -> [[outer] -> [tile]]
 *
 * where tile is defined by a concatenation of the time_tile and
 * the space_tile.
 */
static __isl_give isl_map *construct_tile_map(__isl_keep ppcg_ht_phase *phase)
{
        int depth;
        isl_space *space;
        isl_multi_aff *ma;
        isl_multi_aff *tiling;
        isl_map *el2tile;

        depth = isl_schedule_node_get_schedule_depth(
                                                phase->tiling->input_node);
        space = isl_aff_get_space(phase->time_tile);
        space = isl_space_params(space);
        space = isl_space_set_from_params(space);
        space = isl_space_add_dims(space, isl_dim_set, depth);
        space = isl_space_map_from_set(space);
        ma = isl_multi_aff_identity(space);

        tiling = isl_multi_aff_flat_range_product(
                isl_multi_aff_from_aff(isl_aff_copy(phase->time_tile)),
                isl_multi_aff_copy(phase->space_tile));
        el2tile = isl_map_from_multi_aff(tiling);
        el2tile = isl_map_intersect_domain(el2tile,
                                                isl_set_copy(phase->domain));
        el2tile = isl_map_product(isl_map_from_multi_aff(ma), el2tile);

        return el2tile;
}

/* Return a description of the full tiles of "phase" at the point
 * in the original schedule tree where the tiling was applied.
 *
 * First construct a mapping from the input schedule dimensions
 * up to and including the original pair of bands to which hybrid tiling
 * was applied to schedule dimensions in which this original pair
 * has been replaced by the tiles.
 * This mapping is of the form
 *
 *      [[outer] -> [P -> C]] -> [[outer] -> [tile]]
 *
 * Apply this mapping to the set of all values for the input
 * schedule dimensions and then apply its inverse.
 * The result is the set of values for the input schedule dimensions
 * that would map to any of the tiles.  Subtracting from this set
 * the set of values that are actually executed produces the set
 * of values that belong to a tile but that are not executed.
 * Mapping these back to the tiles produces a description of
 * the partial tiles.  Subtracting these from the set of all tiles
 * produces a description of the full tiles in the form
 *
 *      [[outer] -> [tile]]
 */
static __isl_give isl_set *compute_full_tile(__isl_keep ppcg_ht_phase *phase)
{
        isl_schedule_node *node;
        isl_union_set *domain;
        isl_union_map *prefix, *schedule;
        isl_set *all, *partial, *all_el;
        isl_map *tile2el, *el2tile;
        isl_multi_union_pw_aff *mupa;

        el2tile = construct_tile_map(phase);
        tile2el = isl_map_reverse(isl_map_copy(el2tile));

        node = phase->tiling->input_node;
        prefix = isl_schedule_node_get_prefix_schedule_union_map(node);
        domain = isl_schedule_node_get_domain(node);
        mupa = isl_multi_union_pw_aff_copy(phase->tiling->input_schedule);
        schedule = isl_union_map_from_multi_union_pw_aff(mupa);
        schedule = isl_union_map_range_product(prefix, schedule);
        all_el = isl_set_from_union_set(isl_union_set_apply(domain, schedule));
        all_el = isl_set_coalesce(all_el);

        all = isl_set_apply(isl_set_copy(all_el), isl_map_copy(el2tile));

        partial = isl_set_copy(all);
        partial = isl_set_apply(partial, tile2el);
        partial = isl_set_subtract(partial, all_el);
        partial = isl_set_apply(partial, el2tile);

        return isl_set_subtract(all, partial);
}

/* Copy the AST loop types of the non-isolated part to those
 * of the isolated part.
 */
static __isl_give isl_schedule_node *set_isolate_loop_type(
        __isl_take isl_schedule_node *node)
{
        int i, n;

        n = isl_schedule_node_band_n_member(node);
        for (i = 0; i < n; ++i) {
                enum isl_ast_loop_type type;

                type = isl_schedule_node_band_member_get_ast_loop_type(node, i);
                node = isl_schedule_node_band_member_set_isolate_ast_loop_type(
                                                                node, i, type);
        }

        return node;
}

/* If options->isolate_full_tiles is set, then mark the full tiles
 * in "node" for isolation.  The full tiles are derived from "phase".
 * "node" may point to a part of the tiling, e.g., the space tiling.
 *
 * The full tiles are originally computed in the form
 *
 *      [[outer] -> [tile]]
 *
 * However, the band that "node" points to may only contain
 * subset of the tile dimensions.
 * The description above is therefore treated as
 *
 *      [[outer] -> [before; this; after]]
 *
 * before is of size "pos"; this is of size "dim"; and
 * after is of size "out - pos - dim".
 * The after part is first project out.  Then the range is split
 * into a before and this part and finally the before part is moved
 * to the domain, resulting in
 *
 *      [[outer; before] -> [this]]
 *
 * This description is then used as the isolate option.
 *
 * The AST loop type for the isolated part is set to be the same
 * as that of the non-isolated part.
 */
static __isl_give isl_schedule_node *ppcg_ht_phase_isolate_full_tile_node(
        __isl_keep ppcg_ht_phase *phase, __isl_take isl_schedule_node *node,
        struct ppcg_options *options)
{
        int in, out, pos, depth, dim;
        isl_space *space;
        isl_multi_aff *ma1, *ma2;
        isl_set *tile;
        isl_map *map;
        isl_set *set;
        isl_union_set *opt;

        if (!options->isolate_full_tiles)
                return node;

        depth = isl_schedule_node_get_schedule_depth(node);
        dim = isl_schedule_node_band_n_member(node);

        tile = compute_full_tile(phase);
        map = isl_set_unwrap(tile);
        in = isl_map_dim(map, isl_dim_in);
        out = isl_map_dim(map, isl_dim_out);
        pos = depth - in;
        map = isl_map_project_out(map, isl_dim_out, pos + dim,
                                out - (pos + dim));
        space = isl_space_range(isl_map_get_space(map));
        ma1 = isl_multi_aff_project_out_map(isl_space_copy(space),
                                           isl_dim_set, pos, dim);
        ma2 = isl_multi_aff_project_out_map(space, isl_dim_set, 0, pos);
        ma1 = isl_multi_aff_range_product(ma1, ma2);
        map = isl_map_apply_range(map, isl_map_from_multi_aff(ma1));
        map = isl_map_uncurry(map);
        map = isl_map_flatten_domain(map);
        set = isl_map_wrap(map);
        set = isl_set_set_tuple_name(set, "isolate");

        opt = isl_schedule_node_band_get_ast_build_options(node);
        opt = isl_union_set_add_set(opt, set);
        node = isl_schedule_node_band_set_ast_build_options(node, opt);
        node = set_isolate_loop_type(node);

        return node;
}

/* Insert a band node for performing the space tiling for "phase" at "node".
 * In particular, insert a band node with partial schedule
 *
 *      [P[t] -> C[s]] -> C[floor((s + space_shift)/space_size)]
 *
 * pulled back over the input schedule.
 * "options" determines whether full tiles should be separated
 * from partial tiles.
 *
 * The first tile dimension iterates over the hexagons in the same
 * phase, which are independent by construction.  The first dimension
 * is therefore marked coincident.
 * All dimensions are also marked for being generated as atomic loops
 * because separation is usually not desirable on tile loops.
 */
static __isl_give isl_schedule_node *insert_space_tiling(
        __isl_keep ppcg_ht_phase *phase, __isl_take isl_schedule_node *node,
        struct ppcg_options *options)
{
        isl_multi_aff *space_tile;
        isl_multi_union_pw_aff *mupa;

        if (!phase)
                return isl_schedule_node_free(node);

        space_tile = isl_multi_aff_copy(phase->space_tile);
        mupa = isl_multi_union_pw_aff_copy(phase->tiling->input_schedule);
        mupa = isl_multi_union_pw_aff_apply_multi_aff(mupa, space_tile);
        node = isl_schedule_node_insert_partial_schedule(node, mupa);
        node = ppcg_set_schedule_node_type(node, isl_ast_loop_atomic);
        node = ppcg_ht_phase_isolate_full_tile_node(phase, node, options);
        node = isl_schedule_node_band_member_set_coincident(node, 0, 1);

        return node;
}

/* Given a pointer "node" to (a copy of) the original child node
 * in the input pattern, adjust its partial schedule such that
 * it starts at zero within each tile.
 *
 * That is, replace "s" by (s + space_shift) % space_sizes.
 */
__isl_give isl_schedule_node *ppcg_ht_phase_shift_space_point(
        __isl_keep ppcg_ht_phase *phase, __isl_take isl_schedule_node *node)
{
        isl_multi_val *space_sizes;
        isl_multi_aff *space_shift;
        isl_multi_union_pw_aff *mupa;

        space_shift = isl_multi_aff_copy(phase->space_shift);
        mupa = isl_multi_union_pw_aff_copy(phase->tiling->input_schedule);
        mupa = isl_multi_union_pw_aff_apply_multi_aff(mupa, space_shift);
        node = isl_schedule_node_band_shift(node, mupa);
        space_sizes = isl_multi_val_copy(phase->tiling->space_sizes);
        node = isl_schedule_node_band_mod(node, space_sizes);

        return node;
}

/* Does
 *
 *      s0 > delta + 2 * {delta * h} - 1
 *
 * hold?
 */
static isl_bool wide_enough(__isl_keep isl_val *s0, __isl_keep isl_val *delta,
        __isl_keep isl_val *h)
{
        isl_val *v, *v2;
        isl_bool ok;

        v = isl_val_mul(isl_val_copy(delta), isl_val_copy(h));
        v2 = isl_val_floor(isl_val_copy(v));
        v = isl_val_sub(v, v2);
        v = isl_val_mul_ui(v, 2);
        v = isl_val_add(v, isl_val_copy(delta));
        v = isl_val_sub_ui(v, 1);
        ok = isl_val_gt(s0, v);
        isl_val_free(v);

        return ok;
}

/* Is the tile size specified by "sizes" wide enough in the first space
 * dimension, i.e., the base of the hexagon?  This ensures that,
 * after hybrid tiling using "bounds" and these sizes,
 * neighboring hexagons in the same phase are far enough apart
 * that they do not depend on each other.
 * The test is only meaningful if the bounds are valid.
 *
 * Let st be (half) the size in the time dimension and s0 the base
 * size in the first space dimension.  Let delta be the dependence
 * distance in either positive or negative direction.  In principle,
 * it should be enough to have s0 + 1 > delta, i.e., s0 >= delta.
 * However, in case of fractional delta, the tile is not extended
 * with delta * (st - 1), but instead with floor(delta * (st - 1)).
 * The condition therefore needs to be adjusted to
 *
 *      s0 + 1 > delta + 2 {delta * (st - 1)}
 *
 * (with {} the fractional part) to account for the two slanted sides.
 * The condition in the paper "Hybrid Hexagonal/Classical Tiling for GPUs"
 * translates to
 *
 *      s0 >= delta + {delta * (st - 1)}
 *
 * Since 1 > frac(delta * (st - 1)), this condition implies
 * the condition above.
 *
 * The condition is checked for both directions.
 */
isl_bool ppcg_ht_bounds_supports_sizes(__isl_keep ppcg_ht_bounds *bounds,
        __isl_keep isl_multi_val *sizes)
{
        isl_val *s0, *h;
        isl_val *delta;
        isl_bool ok;

        ok = ppcg_ht_bounds_is_valid(bounds);
        if (ok < 0 || !ok)
                return ok;

        h = isl_val_sub_ui(isl_multi_val_get_val(sizes, 0), 1);
        s0 = isl_multi_val_get_val(sizes, 1);

        delta = ppcg_ht_bounds_get_lower(bounds, 0);
        ok = wide_enough(s0, delta, h);
        isl_val_free(delta);

        delta = ppcg_ht_bounds_get_upper(bounds);
        if (ok == isl_bool_true)
                ok = wide_enough(s0, delta, h);
        isl_val_free(delta);

        isl_val_free(s0);
        isl_val_free(h);

        return ok;
}

/* Check that the tile will be wide enough in the first space
 * dimension, i.e., the base of the hexagon.  This ensures that
 * neighboring hexagons in the same phase are far enough apart
 * that they do not depend on each other.
 *
 * Error out if the condition fails to hold.
 */
static isl_stat check_width(__isl_keep ppcg_ht_bounds *bounds,
        __isl_keep isl_multi_val *sizes)
{
        isl_bool ok;

        ok = ppcg_ht_bounds_supports_sizes(bounds, sizes);

        if (ok < 0)
                return isl_stat_error;
        if (!ok)
                isl_die(isl_multi_val_get_ctx(sizes), isl_error_invalid,
                        "base of hybrid tiling hexagon not sufficiently wide",
                        return isl_stat_error);

        return isl_stat_ok;
}

/* Given valid bounds on the relative dependence distances for
 * the pair of nested nodes that "node" point to, as well as sufficiently
 * wide tile sizes "sizes", insert the corresponding time and space tiling
 * at "node", along with a pair of phase nodes that can be used
 * to make further changes.
 * The space of "sizes" should be the product of the spaces
 * of the schedules of the pair of parent and child nodes.
 * "options" determines whether full tiles should be separated
 * from partial tiles.
 *
 * In particular, given an input of the form
 *
 *      P - C - ...
 *
 * the output has the form
 *
 *              /- F0 - M0 - CT0 - P - C - ...
 *      PT - seq
 *              \- F1 - M1 - CT1 - P - C - ...
 *
 * PT is the global time tiling.  Within each of these tiles,
 * two phases are executed in order.  Within each phase, the schedule
 * space is further subdivided into tiles through CT0 and CT1.
 * The first dimension of each of these iterates over the hexagons
 * within a phase and these are independent by construction.
 * The F0 and F1 filters filter the statement instances that belong
 * to the corresponding phase.  The M0 and M1 marks contain a pointer
 * to a ppcg_ht_phase object that can be used to perform further changes.
 *
 * After checking that input satisfies the requirements,
 * a data structure is constructed that represents the tiling and
 * two additional data structures are constructed for the two phases
 * of the tiling.  These are then used to define the filters F0 and F1 and
 * combined to construct the time tiling PT.
 * Then the time tiling node PT is inserted, followed by
 * the sequence with the two filters, the CT space tiling nodes and
 * the phase markers M.
 */
__isl_give isl_schedule_node *ppcg_ht_bounds_insert_tiling(
        __isl_take ppcg_ht_bounds *bounds, __isl_take isl_multi_val *sizes,
        __isl_take isl_schedule_node *node, struct ppcg_options *options)
{
        isl_ctx *ctx;
        isl_union_set *phase0;
        isl_union_set *phase1;
        isl_multi_union_pw_aff *input, *dom_time;
        isl_union_pw_multi_aff *upma;
        isl_pw_multi_aff *time;
        isl_union_set_list *phases;
        ppcg_ht_tiling *tiling;
        ppcg_ht_phase *phase_0;
        ppcg_ht_phase *phase_1;

        if (!node || !sizes || !bounds)
                goto error;
        if (check_input_pattern(node) < 0 || check_width(bounds, sizes) < 0)
                goto error;

        ctx = isl_schedule_node_get_ctx(node);

        input = extract_input_schedule(node);

        tiling = ppcg_ht_bounds_construct_tiling(bounds, node, input, sizes);
        phase_0 = ppcg_ht_tiling_compute_phase(tiling, 1);
        phase_1 = ppcg_ht_tiling_compute_phase(tiling, 0);
        time = combine_time_tile(phase_0, phase_1);
        ppcg_ht_tiling_free(tiling);

        upma = isl_union_pw_multi_aff_from_multi_union_pw_aff(
                                        isl_multi_union_pw_aff_copy(input));
        phase0 = isl_union_set_from_set(ppcg_ht_phase_get_domain(phase_0));
        phase0 = isl_union_set_preimage_union_pw_multi_aff(phase0,
                                        isl_union_pw_multi_aff_copy(upma));
        phase1 = isl_union_set_from_set(ppcg_ht_phase_get_domain(phase_1));
        phase1 = isl_union_set_preimage_union_pw_multi_aff(phase1, upma);

        phases = isl_union_set_list_alloc(ctx, 2);
        phases = isl_union_set_list_add(phases, phase0);
        phases = isl_union_set_list_add(phases, phase1);

        dom_time = isl_multi_union_pw_aff_apply_pw_multi_aff(input, time);
        node = isl_schedule_node_insert_partial_schedule(node, dom_time);

        node = isl_schedule_node_child(node, 0);

        node = isl_schedule_node_insert_sequence(node, phases);
        node = isl_schedule_node_child(node, 0);
        node = isl_schedule_node_child(node, 0);
        node = insert_space_tiling(phase_0, node, options);
        node = insert_phase(node, phase_0);
        node = isl_schedule_node_parent(node);
        node = isl_schedule_node_next_sibling(node);
        node = isl_schedule_node_child(node, 0);
        node = insert_space_tiling(phase_1, node, options);
        node = insert_phase(node, phase_1);
        node = isl_schedule_node_parent(node);
        node = isl_schedule_node_parent(node);

        node = isl_schedule_node_parent(node);

        isl_multi_val_free(sizes);
        return node;
error:
        isl_multi_val_free(sizes);
        isl_schedule_node_free(node);
        ppcg_ht_bounds_free(bounds);
        return NULL;
}

/* Given a branch "node" that contains a sequence node with two phases
 * of hybrid tiling as input, call "fn" on each of the two phase marker
 * nodes.
 *
 * That is, the input is as follows
 *
 *               /- F0 - M0 - ...
 *      ... - seq
 *               \- F1 - M1 - ...
 *
 * and "fn" is called on M0 and on M1.
 */
__isl_give isl_schedule_node *hybrid_tile_foreach_phase(
        __isl_take isl_schedule_node *node,
        __isl_give isl_schedule_node *(*fn)(__isl_take isl_schedule_node *node,
                void *user), void *user)
{
        int depth0, depth;

        depth0 = isl_schedule_node_get_tree_depth(node);

        while (node &&
            isl_schedule_node_get_type(node) != isl_schedule_node_sequence)
                node = isl_schedule_node_child(node, 0);

        node = isl_schedule_node_child(node, 0);
        node = isl_schedule_node_child(node, 0);
        if (!node)
                return NULL;
        node = fn(node, user);
        node = isl_schedule_node_parent(node);
        node = isl_schedule_node_next_sibling(node);
        node = isl_schedule_node_child(node, 0);
        if (!node)
                return NULL;
        node = fn(node, user);
        node = isl_schedule_node_parent(node);
        node = isl_schedule_node_parent(node);

        depth = isl_schedule_node_get_tree_depth(node);
        node = isl_schedule_node_ancestor(node, depth - depth0);

        return node;
}

/* This function is called on each of the two phase marks
 * in a hybrid tiling tree.
 * Drop the phase mark at "node".
 */
static __isl_give isl_schedule_node *drop_phase_mark(
        __isl_take isl_schedule_node *node, void *user)
{
        isl_id *id;
        isl_bool is_phase;

        if (isl_schedule_node_get_type(node) != isl_schedule_node_mark)
                return node;

        id = isl_schedule_node_mark_get_id(node);
        is_phase = is_phase_id(id);
        isl_id_free(id);

        if (is_phase < 0)
                return isl_schedule_node_free(node);
        if (is_phase)
                node = isl_schedule_node_delete(node);

        return node;
}

/* Given a branch "node" that contains a sequence node with two phases
 * of hybrid tiling as input, remove the two phase marker nodes.
 *
 * That is, the input is as follows
 *
 *               /- F0 - M0 - ...
 *      ... - seq
 *               \- F1 - M1 - ...
 *
 * and the output is
 *
 *               /- F0 - ...
 *      ... - seq
 *               \- F1 - ...
 */
__isl_give isl_schedule_node *hybrid_tile_drop_phase_marks(
        __isl_take isl_schedule_node *node)
{
        return hybrid_tile_foreach_phase(node, &drop_phase_mark, NULL);
}