Lesson 18.7 — Dependence-driven transformations: fusion, fission, interchange, tiling, the polyhedral model¶

Techniques: loop fusion and fission (distribution), loop interchange and permutation, strip-mining, skewing and tiling, the polyhedral model (Feautrier's dataflow analysis and scheduling; Polly, isl, GCC Graphite) · Pebble implements: — · Lab: — (the drill interchange-legal checks your legality verdicts by executing the permuted nest) · Prerequisites: Lesson 18.6 (distance and direction vectors) · Time: 6–7 hours

Lesson 18.6 computed which iterations depend on which. This lesson uses that information to reorder iterations: to split a loop into several (fission) or merge several into one (fusion), to swap the loops of a nest (interchange), to cut the iteration space into cache-sized blocks (tiling), and — in the polyhedral model — to compute an arbitrary new affine execution order at once. All of them rest on one principle: a reordering is legal if it executes the source of every dependence before its sink. Each technique turns that principle into a concrete test on distance vectors.

1. Problem and motivation¶

The problem. A loop nest written for clarity is rarely in the best order for the machine. for j … for i … A[i][j] walks a row-major array column by column, touching a new cache line in every iteration; two loops over the same array read it from memory twice; a loop mixing a recurrence with independent work cannot be vectorized as a whole. The fix is to reorder the iterations — and the question is always the same: which reorderings preserve the program's result? Dependence distances (Lesson 18.6) answer it.

Loop fusion and fission¶

Fission (loop distribution) splits one loop into several, each containing some of the statements; fusion merges adjacent loops with the same iteration space. Allen and Kennedy's vectorizing compiler PFC distributed loops around the strongly connected components of the statement dependence graph, so that every statement not on a dependence cycle could be vectorized [AK02, Ch. 2]. Kennedy and McKinley studied the reverse, fusion for locality and parallelism, and showed that fusion to maximize data reuse is NP-hard [KM93]. LLVM has both: loop-distribute (driven by LoopAccessAnalysis, to isolate what blocks vectorization) and loop-fusion [LLVM-Distribute; LLVM-Fuse]; GCC distributes loops and recognizes memset/memcpy partitions [GCC-Distribution].

Loop interchange and permutation¶

Swapping two loops of a perfect nest changes the order in which the iteration space is swept. Allen and Kennedy characterized its legality with direction vectors: interchange is illegal exactly when some dependence has direction (<, >) [AK84]. The general rule for any permutation (Theorem 18.7.4) follows. LLVM's loop-interchange is in the default -O2 pipeline of LLVM 23 [LLVM-Interchange]; GCC's -floop-interchange is enabled at -O3 (gcc 14.2.0's -Q --help=optimizers lists it as disabled at -O2) [GCC-Interchange].

Strip-mining, skewing and tiling¶

Strip-mining splits one loop into a loop over blocks and a loop within a block — the shape of every vectorized loop (Lesson 18.8). Tiling (blocking) strip-mines several loops of a nest and moves the block loops outward, so that the iterations of one tile reuse data while it is still in cache. Irigoin and Triolet gave the legality condition for tiling ("supernode partitioning") [IT88]; Wolf and Lam unified interchange, reversal and skewing as unimodular transformations and used skewing to make tiling legal [WL91].

The polyhedral model¶

Instead of choosing a sequence of loop transformations, the polyhedral model represents each statement's iterations as the integer points of a polyhedron, computes exact dependences as relations between polyhedra — Feautrier's array dataflow analysis even finds the unique last write that produced each read value [Fea91] — and then solves for a new affine schedule that respects all dependences at once [Fea92]. Pluto chooses schedules that make tiling and parallelization possible [BHRS08]; isl implements the integer-set operations and a Pluto-style scheduler [Ver10]; Polly [GGL12] and GCC's Graphite bring the model into LLVM and GCC.

2. Definitions and algorithms¶

Throughout, nests are perfect with constant bounds and step 1 (Definition 18.6.1), a distance vector \(\mathbf{d}\) is sink minus source (Definition 18.6.3), and \(\mathbf{d} \succ \mathbf{0}\) means lexicographically positive.

Loop fusion and fission¶

Definition 18.7.1 (Statement dependence graph; fission; fusion)

For a loop with body statements \(S_1, \dots, S_m\) (in textual order), the statement dependence graph has an edge \(S_p \to S_q\) labelled with distance \(d\) for each dependence from an instance of \(S_p\) to an instance of \(S_q\) \(d\) iterations later (\(d = 0\): in the same iteration, then \(p < q\) or the edge is a read-before-write inside one statement). Fission (distribution) replaces the loop by several loops over the same range, each executing a subsequence of the statements; fusion replaces two adjacent loops \(L_1; L_2\) over the same range by one loop whose body is \(L_1\)'s body followed by \(L_2\)'s. A dependence from \(L_1\)'s instance \(i\) to \(L_2\)'s instance \(i'\) is fusion-preventing if \(i' < i\): after fusion the sink would run before the source.

Algorithm 18.7.2 (Loop distribution by strongly connected components — Allen and Kennedy)

Input: a loop \(L\) with statements \(S_1..S_m\), its statement dependence graph \(G\).
Output: a sequence of loops, one per strongly connected component (SCC) of \(G\), each marked parallel (acyclic component) or sequential.
Precondition: all dependences of \(L\) are in \(G\) (sound dependence analysis, Lesson 18.6); the body is straight-line (or if-converted).
Postcondition: executing the loops in the output order gives the same result as \(L\) (Theorem 18.7.14).
Invariant: every edge of \(G\) between different components goes from an earlier emitted loop to a later one.

function Distribute(L, G):
    C ← the SCCs of G (Tarjan, Ch 15); form the condensation DAG
    order ← a topological order of the DAG, ties broken by textual order
    for each component c in order:
        emit  for i in range(L): the statements of c, in textual order
        mark it parallel if c is a single statement without a self-edge of distance ≠ 0

Algorithm 18.7.3 (Loop fusion)

Input: adjacent loops \(L_1; L_2\) with the same trip count and no code between them (or code independent of both).
Output: the fused loop, or "not fused".
Precondition: both are in loop-simplify form; dependence analysis can relate their accesses.
Postcondition: the fused loop computes the same result (Theorem 18.7.15).
Invariant: —

function Fuse(L1, L2):
    if trip counts differ or L1 and L2 are not control-flow equivalent: return not fused
    for each pair (a in L1, b in L2) of memory accesses, at least one a write:
        D ← dependence distances from a (iteration i of L1) to b (iteration i' of L2), as i' − i
        if some distance in D is negative: return not fused       # fusion-preventing
    body ← L1.body ; L2.body  with L2's induction variable replaced by L1's
    return the loop with this body

LLVM's loop-distribute isolates a recurrence so that the rest vectorizes

Reproduce (clang 23.1.2, opt 23.1.2):

cat > dist.c <<'X'
void dist(long *restrict a, long *restrict b, long *restrict c,
          long *restrict d, long n) {
  for (long i = 0; i < n; i++) {
    a[i + 1] = a[i] * b[i];   /* S1: a recurrence through memory */
    c[i] = d[i] * 3;          /* S2: independent of S1 */
  }
}
X
clang-23 -O0 -Xclang -disable-O0-optnone -fno-discard-value-names -S -emit-llvm dist.c -o - |
  opt -passes='mem2reg,instcombine,loop-simplify,loop-rotate' -S -o dist.ll
opt -passes='loop-distribute' -enable-loop-distribute -pass-remarks=loop-distribute \
  -S dist.ll -o dist.out.ll 2>&1
awk '/^[a-z.0-9_]+:/ { b = $1 } /store/ { print b, $0 }' dist.out.ll
opt -passes='loop-distribute,loop-vectorize' -enable-loop-distribute \
  -pass-remarks=loop-vectorize -pass-remarks-analysis=loop-vectorize -disable-output dist.ll 2>&1

Output (complete):

remark: <unknown>:0:0: distributed loop
for.body.ldist1:   store i64 %mul.ldist1, ptr %arrayidx2.ldist1, align 8
for.body:   store i64 %mul4, ptr %arrayidx5, align 8
remark: <unknown>:0:0: loop not vectorized: unsafe dependent memory operations in loop. Use #pragma clang loop distribute(enable) to allow loop distribution to attempt to isolate the offending operations into a separate loop
Backward loop carried data dependence.
remark: <unknown>:0:0: vectorized loop (vectorization width: 2, interleaved count: 2)

What to notice: S1 is on a dependence cycle (its store feeds the next iteration's load, distance 1); S2 is on none. Algorithm 18.7.2 gives two loops: for.body.ldist1 holds S1's store, the original loop keeps S2's. The vectorizer then rejects the first ("Backward loop carried data dependence") and vectorizes the second. LLVM's distribution is off by default (-enable-loop-distribute, or #pragma clang loop distribute(enable) per loop, as the vectorizer's remark suggests) and only splits loops for this purpose: it keeps the textual order of memory operations, because LoopAccessAnalysis presumes it [LLVM-Distribute].

LLVM's loop-fusion: fused, and prevented by a dependence

Reproduce (clang 23.1.2, opt 23.1.2):

cat > fuse.c <<'X'
void fuse(long *restrict a, long *restrict b, long *restrict c) {
  for (long i = 0; i < 1000; i++)
    a[i] = b[i] + 1;
  for (long i = 0; i < 1000; i++)
    c[i] = a[i] * 2;
}
void nofuse(long *restrict a, long *restrict b, long *restrict c) {
  for (long i = 0; i < 1000; i++)
    a[i] = b[i] + 1;
  for (long i = 0; i < 1000; i++)
    c[i] = a[i + 1] * 2;
}
X
clang-23 -O0 -Xclang -disable-O0-optnone -fno-discard-value-names -S -emit-llvm fuse.c -o - |
  opt -passes='mem2reg,instcombine,loop-simplify,loop-rotate' -S -o fuse.ll
opt -passes='loop-fusion,print<loops>' -pass-remarks=loop-fusion \
  -pass-remarks-analysis=loop-fusion -disable-output fuse.ll 2>&1

Output (complete):

remark: <unknown>:0:0: [fuse]: entry and for.end: Loops fused
Loop info for function 'fuse':
Loop at depth 1 containing: %for.body<header>,%for.inc,%for.inc8<latch><exiting>
Loop info for function 'nofuse':
Loop at depth 1 containing: %for.body5<header>,%for.inc9<latch><exiting>
Loop at depth 1 containing: %for.body<header>,%for.inc<latch><exiting>

What to notice: in fuse the second loop reads a[i], written by the first loop in the same iteration \(i\): distance 0, fusion legal, one loop remains. In nofuse the second loop's iteration \(i\) reads a[i+1], which the first loop writes in iteration \(i + 1\): distance \(-1\), fusion-preventing (Definition 18.7.1) — after fusion iteration \(i\) would read the old a[i+1]. Both loops remain.

GCC distributes a loop and turns one partition into memset

Reproduce (gcc 14.2.0):

cat > gdist.c <<'X'
void gdist(long *restrict a, long *restrict b, long *restrict c, long n) {
  for (long i = 0; i < n; i++) {
    a[i + 1] = a[i] * b[i];
    c[i] = 0;
  }
}
X
gcc-14 -O2 -ftree-loop-distribution -fopt-info-loop-optimized -c gdist.c -o /dev/null

Output (complete):

gdist.c:2:22: optimized: Loop 1 distributed: split to 1 loops and 1 library calls.

What to notice: GCC's distribution builds a reduced dependence graph (RDG) of the loop's statements, partitions it, and classifies each partition; a partition that only stores zero is replaced by a library call. One loop (the recurrence) and one memset remain [GCC-Distribution].

Loop interchange and permutation¶

A loop permutation \(\pi\) of a depth-\(d\) nest makes old loop \(\pi(k)\) the new \(k\)-th loop (outermost first); interchange of loops \(k\) and \(k + 1\) is the transposition. The new nest executes the same set of iterations, ordered lexicographically by the permuted vector \(\pi(\mathbf{i}) = (i_{\pi(1)}, \dots, i_{\pi(d)})\).

Theorem 18.7.4 (Legality of a loop permutation)

A loop permutation \(\pi\) of a perfect nest (the body executed as a unit) preserves the order of every dependence — and therefore the result of the nest — if and only if for every distance vector \(\mathbf{d} \ne \mathbf{0}\) of a dependence between different iterations, the permuted vector \(\pi(\mathbf{d})\) is lexicographically positive. (Zero distances are loop-independent: they stay inside one execution of the body, whose statement order is unchanged.)

Proof

If. Let a dependence go from iteration \(\mathbf{x}\) to iteration \(\mathbf{y} = \mathbf{x} + \mathbf{d}\), \(\mathbf{d} \ne \mathbf{0}\). The permuted nest runs \(\mathbf{x}\) before \(\mathbf{y}\) iff \(\pi(\mathbf{x}) \prec \pi(\mathbf{y})\), i.e. iff \(\pi(\mathbf{y}) - \pi(\mathbf{x}) = \pi(\mathbf{d}) \succ \mathbf{0}\) (permutation is linear and lexicographic order is translation-invariant). So every dependent pair of iterations keeps its order; independent pairs may swap, which cannot change any value read or the final contents of memory (they touch disjoint locations or only read). By induction over the executed instances, every read sees the same last write as before, so every value and the final state are the same.

Only if. If some dependence \(\mathbf{x} \to \mathbf{x} + \mathbf{d}\) has \(\pi(\mathbf{d}) \prec \mathbf{0}\) (\(\pi(\mathbf{d}) \ne \mathbf{0}\) since \(\mathbf{d} \ne \mathbf{0}\)), the permuted nest runs the sink before the source. For a flow dependence the read then sees an older value; for an anti dependence it sees the newer one; for an output dependence the final value comes from the other write. So the dependence is not preserved. Whether the result changes then depends on the values: it does for some array contents whenever the two instances can carry different values (as in A[f] = A[g] + 1), but not, for example, when both writes of an output dependence store the same constant. That is why the "only if" is stated for the order of dependences; for the result it holds for every dependence that can carry distinct values. ∎

With direction vectors the condition reads: after permuting, the leftmost non-= direction of every vector must be <. Interchange of two loops is illegal exactly when a dependence has (<, >) in those positions and = in all outer ones [AK84].

Algorithm 18.7.5 (Loop interchange — legality and profitability, as in LLVM)

Input: a perfect nest of depth \(d\) (innermost to outermost), its dependences.
Output: the nest with some adjacent loops swapped.
Precondition: tightly nested loops (no code between the headers except the IV updates), a computable trip count for each, reductions and LCSSA phis recognized.
Postcondition: each swap satisfies Theorem 18.7.4; the result computes the same values.
Invariant: the direction matrix (one row per dependence, one column per loop) is permuted along with the loops.

function Interchange(nest):
    M ← direction matrix from DependenceAnalysis; give up on any unknown (*) row that cannot be refined
    for inner ← innermost loop up to the second outermost:           # bubble outward
        outer ← the loop enclosing inner
        if some row r has M[r] lexicographically ≺ 0 after swapping columns inner and outer:
            continue                                                   # illegal (Theorem 18.7.4)
        if not profitable(inner, outer): continue                      # cache-line reuse, vectorizability
        swap the loops in the IR (move headers and latches; fix phis); swap columns in M

LLVM's loop-interchange fixes a column-major sweep and refuses the illegal one

Reproduce (clang 23.1.2, opt 23.1.2):

cat > ic.c <<'X'
long A[256][256], B[256][256];
void colmajor(void) {           /* inner loop strides by a whole row */
  for (int j = 0; j < 256; j++)
    for (int i = 0; i < 256; i++)
      A[i][j] = A[i][j] + B[i][j];
}
void skewed(void) {             /* distance (1, -1): interchange illegal */
  for (int i = 1; i < 256; i++)
    for (int j = 0; j < 255; j++)
      A[i][j] = A[i - 1][j + 1] + 1;
}
X
clang-23 -O1 -gline-tables-only -fno-vectorize -fno-unroll-loops -S -emit-llvm ic.c -o ic.ll
opt -passes='loop-interchange' -pass-remarks=loop-interchange \
  -pass-remarks-missed=loop-interchange -disable-output ic.ll 2>&1
opt -passes='print<da>' -disable-output ic.ll 2>&1 | grep -A1 'load.*--> Dst:  store'

Output (complete):

remark: ic.c:4:5: Loop interchanged with enclosing loop.
remark: ic.c:9:5: Cannot interchange loops due to dependences.
Src:  %12 = load i64, ptr %11, align 8, !dbg !28, !tbaa !29 --> Dst:  store i64 %15, ptr %11, align 8, !dbg !33, !tbaa !29
  da analyze - anti [0 0|<]!
--
Src:  %14 = load i64, ptr %13, align 8, !dbg !31, !tbaa !29 --> Dst:  store i64 %15, ptr %11, align 8, !dbg !33, !tbaa !29
  da analyze - none!
--
Src:  %13 = load i64, ptr %12, align 8, !dbg !28, !tbaa !29 --> Dst:  store i64 %14, ptr %15, align 8, !dbg !33, !tbaa !29
  da analyze - anti [-1 1]!

What to notice: colmajor has only a loop-independent dependence (0 0|<: the load and store of A[i][j] in one iteration): interchange is legal and profitable (the new inner loop walks consecutive addresses). In skewed, DA prints the pair in textual order (load, then store) and so reports the anti view [-1 1]; the underlying flow dependence goes from the store in iteration \((i, j)\) to the load in \((i + 1, j - 1)\): distance \((1, -1)\), direction (<, >). Interchanged it would be \((-1, 1) \prec \mathbf{0}\): illegal, as the remark says.

GCC interchanges the same loop (and then vectorizes it)

Reproduce (gcc 14.2.0):

cat > ic.c <<'X'
long A[256][256], B[256][256];
void colmajor(void) {
  for (int j = 0; j < 256; j++)
    for (int i = 0; i < 256; i++)
      A[i][j] = A[i][j] + B[i][j];
}
X
gcc-14 -O2 -floop-interchange -fopt-info-loop-optimized -c ic.c -o /dev/null

Output (complete):

ic.c:3:21: optimized: loops interchanged in loop nest
ic.c:4:23: optimized: loop vectorized using 16 byte vectors

What to notice: GCC's linterchange pass checks the data dependences of the nest (valid_data_dependences) and a stride-based cost model (should_interchange_loops) [GCC-Interchange]; after interchange the inner loop is unit-stride, which is what made it vectorizable (Lesson 18.8).

Strip-mining, skewing and tiling¶

Definition 18.7.6 (Strip-mining, tiling, skewing, fully permutable band)

Strip-mining loop \(i \in [L, U]\) by \(B\) replaces it by for \(t = \lfloor L/B \rfloor\) to \(\lfloor U/B \rfloor\) : for \(i = \max(L, tB)\) to \(\min(U, tB + B - 1)\): a tile loop \(t\) and a point loop \(i\). Tiling a band of consecutive loops \(k..m\) strip-mines each of them and orders all tile loops outside all point loops: iteration \(\mathbf{i}\) runs at position \((\lfloor i_k/B_k \rfloor, \dots, \lfloor i_m/B_m \rfloor, i_k, \dots, i_m)\) inside the outer loops. Skewing loop \(j\) by factor \(f\) with respect to an outer loop \(i\) renames \(j' = j + f i\) (bounds adjusted); a distance \((d_i, d_j)\) becomes \((d_i, d_j + f d_i)\). A band of loops \(k..m\) is fully permutable if every dependence not already carried by a loop outside the band (that is, whose distance is zero in all outer positions) has \(d_k, \dots, d_m \ge 0\).

Algorithm 18.7.7 (Tiling a band, with skewing to make it fully permutable)

Input: a perfect nest, its (uniform) distance vectors, a band \(k..m\), tile sizes \(B_k..B_m\).
Output: the tiled nest.
Precondition: the dependences are uniform (constant distances) or summarized by distance vectors with known signs.
Postcondition: the tiled nest computes the same result (Theorem 18.7.16).
Invariant: skewing keeps every distance lexicographically positive (Lemma 18.7.8).

function Tile(nest, D, band, B):
    for q in band (outer to inner), for each earlier p in band:
        f ← max over d in D with d_q < 0 and d_p > 0 of ⌈-d_q / d_p⌉   (0 if none)
        if f > 0: skew loop q by f with respect to p; update every d in D (d_q += f·d_p)
    if some d in D still has a negative band component: return nest    # cannot tile
    strip-mine each band loop by its B; move the tile loops outside the point loops
    return the new nest (bounds from Fourier–Motzkin on the tile/point constraints)

Lemma 18.7.8 (Skewing preserves legality)

Skewing an inner loop \(j\) by \(f\) with respect to an outer loop \(i\) keeps every lexicographically positive distance vector lexicographically positive, so it is always legal; with \(f \ge \max \lceil -d_j / d_i \rceil\) over the dependences with \(d_i > 0\), all their \(j\)-components become non-negative.

Proof

The components before position \(j\) are unchanged. If one of them is non-zero, the first non-zero one is still first and still positive. Otherwise all components before \(j\) are zero, in particular \(d_i = 0\), so \(d_j + f d_i = d_j\) is unchanged too, and the vector is unchanged. For the second claim: if \(d_i > 0\) then \(d_j + f d_i \ge d_j + \lceil -d_j/d_i \rceil d_i \ge 0\); if \(d_i = 0\) the \(j\)-component is unchanged, and if it was negative the vector's first non-zero component lies before \(i\) (it is positive), so that dependence is carried outside and does not constrain a band beginning at \(i\).

Clang's #pragma omp tile: strip-mining and tiling as loop structure

Reproduce (clang 23.1.2, opt 23.1.2):

cat > tile.c <<'X'
void mm(int n, double C[n][n], double A[n][n], double B[n][n]) {
  #pragma omp tile sizes(32, 32)
  for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
      for (int k = 0; k < n; k++)
        C[i][j] += A[i][k] * B[k][j];
}
void strip(int n, double *a) {
  #pragma omp tile sizes(4)
  for (int i = 0; i < n; i++)
    a[i] = 2 * a[i];
}
X
clang-23 -fopenmp -O1 -fno-vectorize -fno-unroll-loops -fno-discard-value-names \
  -S -emit-llvm tile.c -o tile.ll
opt -passes='print<loops>' -disable-output tile.ll 2>&1 |
  sed -E 's/(<header>).*/\1/; s/containing: %[^<]*,%/containing: %/'

Output (complete):

Loop info for function 'mm':
Loop at depth 1 containing: %for.cond8.preheader<header>
    Loop at depth 2 containing: %for.cond13.preheader<header>
        Loop at depth 3 containing: %for.cond23.preheader<header>
            Loop at depth 4 containing: %for.cond38.preheader<header>
                Loop at depth 5 containing: %for.body41<header>
Loop info for function 'strip':
Loop at depth 1 containing: %for.cond3.preheader<header>
    Loop at depth 2 containing: %for.body11<header>

What to notice: tiling the \(i, j\) band by \(32 \times 32\) turns the 3-deep nest into 5 loops (two tile loops, two point loops, then \(k\)); sizes(4) on a single loop is strip-mining: 2 loops. The directive is a programmer's assertion — clang tiles without checking dependences (OpenMP makes the programmer responsible), which is exactly why a compiler that tiles on its own needs Theorem 18.7.16. LLVM itself has no automatic tiling pass outside Polly. (With -fopenmp-simd instead of -fopenmp, this clang 23.1.2 build crashed on the tile directive; -fopenmp works.)

The polyhedral model¶

Definition 18.7.9 (Static control part: domains, accesses, schedules)

A static control part (SCoP) is a region whose loop bounds and conditions are affine in outer loop indices and parameters, and whose array subscripts are affine. Each statement \(S\) has an iteration domain \(\mathcal{D}_S = \{\mathbf{i} \in \mathbb{Z}^{d_S} : A_S \mathbf{i} + \mathbf{b}_S \ge 0\}\), access relations \(\{S[\mathbf{i}] \to X[F\mathbf{i} + \mathbf{f}]\}\) for its reads and writes, and a schedule \(\theta_S : \mathcal{D}_S \to \mathbb{Z}^p\), an affine map; instances execute in lexicographic order of their schedule values (ties broken by textual order). The original program's schedule interleaves loop indices with textual positions: \(S[\mathbf{i}] \mapsto (0, i_1, 0, i_2, 1)\) and so on.

Definition 18.7.10 (Dependence relation; dataflow)

The dependence relation between statements \(S\) and \(T\) through array \(X\) is \(\{S[\mathbf{x}] \to T[\mathbf{y}] : \mathbf{x} \in \mathcal{D}_S, \mathbf{y} \in \mathcal{D}_T,\ F_S\mathbf{x} + \mathbf{f}_S = F_T\mathbf{y} + \mathbf{f}_T,\ \theta_S(\mathbf{x}) \prec \theta_T(\mathbf{y})\}\), at least one access a write — a union of integer polyhedra, which is where the name comes from. The dataflow (value-based, or exact flow) relation keeps, for each read \(T[\mathbf{y}]\), only the last write executing before it: \(\mathrm{src}(\mathbf{y}) = \operatorname{lexmax}_{\theta} \{S[\mathbf{x}] : \text{writes the element } T[\mathbf{y}] \text{ reads},\ \theta_S(\mathbf{x}) \prec \theta_T(\mathbf{y})\}\).

Algorithm 18.7.11 (Array dataflow analysis — Feautrier)

Input: a SCoP; a read access of statement \(T\).
Output: for each instance \(\mathbf{y}\) of \(T\), the write instance that produced the value it reads (or "from before the SCoP"), as a piecewise affine function of \(\mathbf{y}\) (a quast, quasi-affine selection tree).
Precondition: affine domains and subscripts.
Postcondition: the output is exact: it is the dataflow relation of Definition 18.7.10.
Invariant: the candidate source for each region of \(\mathbf{y}\) is the latest among the writes considered so far.

function Dataflow(T, read):
    result ← ⊥ (no source) for all y in D_T
    for each write access W of statement S to the same array:
        for each depth ℓ at which S and T share loops (and the textual case):
            P_ℓ(y) ← lexmax{ x ∈ D_S : F_S x + f_S = F_T y + f_T, x precedes y at depth ℓ }
                    # a parametric integer program in the parameters y (PIP)
        cand(y) ← the latest of the P_ℓ(y) that exist
        result(y) ← the later of result(y) and cand(y)                # piecewise, by case splits on y
    return result

Theorem 18.7.12 (Affine form of Farkas' lemma)

Let \(\mathcal{D} = \{\mathbf{x} \in \mathbb{R}^n : A\mathbf{x} + \mathbf{b} \ge \mathbf{0}\}\) be non-empty. An affine function \(\phi(\mathbf{x}) = \mathbf{c}^T\mathbf{x} + c_0\) is non-negative everywhere on \(\mathcal{D}\) if and only if \(\phi(\mathbf{x}) \equiv \lambda_0 + \boldsymbol{\lambda}^T (A\mathbf{x} + \mathbf{b})\) for some \(\lambda_0 \ge 0\), \(\boldsymbol{\lambda} \ge \mathbf{0}\) (identically in \(\mathbf{x}\)) [Sch86, Cor. 7.1h].

Proof

If: on \(\mathcal{D}\) every term is a product of non-negative numbers. Only if: \(\phi \ge 0\) on \(\mathcal{D}\) means the linear program \(\min\{\mathbf{c}^T\mathbf{x} : A\mathbf{x} \ge -\mathbf{b}\}\) is feasible and bounded below by \(-c_0\). By LP strong duality its dual \(\max\{-\mathbf{b}^T\boldsymbol{\lambda} : A^T\boldsymbol{\lambda} = \mathbf{c},\ \boldsymbol{\lambda} \ge \mathbf{0}\}\) has an optimal solution \(\boldsymbol{\lambda}^*\) with \(-\mathbf{b}^T\boldsymbol{\lambda}^* = \min \mathbf{c}^T\mathbf{x} \ge -c_0\). Then \(\phi(\mathbf{x}) = \boldsymbol{\lambda}^{*T} A\mathbf{x} + c_0 = \boldsymbol{\lambda}^{*T}(A\mathbf{x} + \mathbf{b}) + (c_0 - \mathbf{b}^T\boldsymbol{\lambda}^*)\), and \(\lambda_0 = c_0 - \mathbf{b}^T\boldsymbol{\lambda}^* \ge 0\).

Algorithm 18.7.13 (One-dimensional affine scheduling — Feautrier)

Input: statements with domains \(\mathcal{D}_S\); dependence polyhedra \(\mathcal{P}_e = \{(\mathbf{x}, \mathbf{y})\}\) for each dependence edge \(e: S \to T\).
Output: affine schedules \(\theta_S(\mathbf{x}) = \mathbf{a}_S^T\mathbf{x} + a_{S,0}\) with \(\theta_T(\mathbf{y}) - \theta_S(\mathbf{x}) \ge 1\) on every \(\mathcal{P}_e\), minimizing latency; or "no one-dimensional schedule" (then a multidimensional one, [Fea92] Part II).
Precondition: non-empty polyhedra given by affine inequalities.
Postcondition: executing all instances with equal \(\theta\) in parallel, in increasing \(\theta\), respects every dependence (Theorem 18.7.17).
Invariant: every constraint added is linear in the unknown coefficients (thanks to Theorem 18.7.12).

function Schedule(statements, edges):
    unknowns ← the coefficients a_S of every θ_S
    for each statement S: require θ_S ≥ 0 on D_S       # Farkas: θ_S ≡ μ_0 + μ^T(A_S x + b_S), μ ≥ 0
    for each edge e: S → T with polyhedron P_e = {(x, y) : C(x, y) + c ≥ 0}:
        require θ_T(y) − θ_S(x) − 1 ≡ λ_0 + λ^T (C(x, y) + c), λ_0, λ ≥ 0
        # equate the coefficients of every x_k, y_k and the constant: linear equalities
    solve the LP (minimize a latency bound); if infeasible: no one-dimensional schedule
    return the θ_S

An exact dependence analysis and a new schedule with isl (the library under Polly and Graphite)

Reproduce (islpy 2026.2.2, isl bundled; uv 0.x installs it on the fly):

cat > poly.py <<'X'
import islpy as isl

# for (i = 1; i <= 8; i++) for (j = 1; j <= 8; j++)
#   S: A[i][j] = A[i-1][j+1] + A[i][j-1];
dom = isl.UnionSet("{ S[i, j] : 1 <= i <= 8 and 1 <= j <= 8 }")
write = isl.UnionMap("{ S[i, j] -> A[i, j] }").intersect_domain(dom)
read = isl.UnionMap("{ S[i, j] -> A[i - 1, j + 1]; S[i, j] -> A[i, j - 1] }").intersect_domain(dom)
orig = isl.UnionMap("{ S[i, j] -> [i, j] }")          # the original execution order

# exact array dataflow (Feautrier): for each read, the last write that produced its value
ai = isl.UnionAccessInfo.from_sink(read).set_must_source(write).set_schedule_map(orig)
flow = ai.compute_flow().get_must_dependence()
print("flow dependences:", flow)
print("distances:", flow.deltas())

# a new schedule: legal (validity) and with short dependences (proximity), as Pluto/isl do
sc = isl.ScheduleConstraints.on_domain(dom).set_validity(flow).set_proximity(flow)
sched = sc.compute_schedule()
m = sched.get_map()
print("new schedule:", m)
print("distances under it:", flow.apply_domain(m).apply_range(m).deltas())
print("interchanged (j, i) legal:",
      flow.apply_domain(isl.UnionMap("{ S[i, j] -> [j, i] }"))
          .apply_range(isl.UnionMap("{ S[i, j] -> [j, i] }"))
          .deltas().is_subset(isl.UnionSet("{ [a, b] : a > 0 or (a = 0 and b > 0) }")))

# tiling 4 x 4: legal iff every dependence stays lexicographically positive in the new order
lexpos = isl.UnionSet("{ [a, b, c, d] : a > 0 or (a = 0 and b > 0) or (a = 0 and b = 0 and c > 0)"
                      " or (a = 0 and b = 0 and c = 0 and d > 0) }")
def legal(s):
    return flow.apply_domain(s).apply_range(s).deltas().is_subset(lexpos)
plain = isl.UnionMap("{ S[i, j] -> [floor(i/4), floor(j/4), i, j] }").intersect_domain(dom)
tiled = isl.UnionMap("{ S[i, j] -> [floor(i/4), floor((i + j)/4), i, i + j] }").intersect_domain(dom)
print("tiling (i, j) legal:", legal(plain))
print("tiling the skewed (i, i + j) legal:", legal(tiled))
ast = isl.AstBuild.from_context(isl.Set("{ : }")).node_from_schedule_map(tiled)
print(ast.to_C_str())
X
uv run --quiet --no-project --with islpy==2026.2.2 python poly.py

Output (complete):

flow dependences: { S[i, j] -> S[i' = i, j' = 1 + j] : 0 < i <= 8 and 0 < j <= 7; S[i, j] -> S[i' = 1 + i, j' = -1 + j] : 0 < i <= 7 and 2 <= j <= 8 }
distances: { S[i = 0, j = 1]; S[i = 1, j = -1] }
new schedule: { S[i, j] -> [i, i + j] }
distances under it: { [0, 1]; [1, 0] }
interchanged (j, i) legal: False
tiling (i, j) legal: False
tiling the skewed (i, i + j) legal: True
for (int c0 = 0; c0 <= 2; c0 += 1)
  for (int c1 = c0; c1 <= c0 + 2; c1 += 1)
    for (int c2 = max(1, 4 * c0); c2 <= min(min(8, 4 * c0 + 3), 4 * c1 + 2); c2 += 1)
      for (int c3 = max(4 * c1, c2 + 1); c3 <= min(4 * c1 + 3, c2 + 8); c3 += 1)
        S(c2, -c2 + c3);

What to notice: compute_flow is Feautrier's dataflow analysis (Algorithm 18.7.11): each read gets its unique last writer, with exact domains (the A[i-1][j+1] value comes from iteration \((i - 1, j + 1)\) only for \(j \ge 2\) and \(i \le 7\) at the source). The distances \((0, 1)\) and \((1, -1)\) forbid interchange (Theorem 18.7.4) and plain tiling (the band is not fully permutable). isl's scheduler (Pluto's algorithm [BHRS08]) finds the skew \((i, i + j)\) of Lemma 18.7.8, after which both distances are non-negative, and \(4 \times 4\) tiling is legal; the AST generator then produces the tiled loops with exact bounds. Polly is not part of the LLVM 23.1.2 build this course uses (opt: unknown pass name 'polly-scops'), so Polly's own output is not shown; it performs the same steps on LLVM IR [GGL12; Polly-Docs].

GCC Graphite: a SCoP, its schedule, and isl's tiled schedule

Reproduce (gcc 14.2.0, built with isl):

cat > gr.c <<'X'
long A[256][256], B[256][256];
void colmajor(void) {
  for (int j = 0; j < 256; j++)
    for (int i = 0; i < 256; i++)
      A[i][j] = A[i][j] + B[i][j];
}
X
gcc-14 -O2 -floop-nest-optimize -fno-loop-interchange -fdump-tree-graphite-details \
  -fopt-info-loop -c gr.c -o /dev/null 2>&1
sed -n '/\[scheduler\] original schedule/,/^$/p;/isl transformed schedule/,/^$/p;/AST generated by isl/,/^$/p' \
  gr.c.*graphite

Output (complete):

gr.c:3:21: optimized: loop nest optimized
[scheduler] original schedule:
domain: "{ S_3[i1, i2] : 0 <= i1 <= 255 and 0 <= i2 <= 255 }"
child:
  schedule: "L_1[{ S_3[i1, i2] -> [(i1)] }]"
  child:
    schedule: "L_2[{ S_3[i1, i2] -> [(i2)] }]"

[scheduler] isl transformed schedule:
domain: "{ S_3[i1, i2] : 0 <= i1 <= 255 and 0 <= i2 <= 255 }"
child:
  schedule: "[{ S_3[i1, i2] -> [(floor((i1)/51))] }, { S_3[i1, i2] -> [(floor((i2)/51))] }]"
  permutable: 1
  coincident: [ 1, 1 ]
  child:
    schedule: "[{ S_3[i1, i2] -> [((i1) mod 51)] }, { S_3[i1, i2] -> [((i2) mod 51)] }]"
    permutable: 1
    coincident: [ 1, 1 ]

[scheduler] AST generated by isl:
for (int c0 = 0; c0 <= 5; c0 += 1)
  for (int c1 = 0; c1 <= 5; c1 += 1)
    for (int c2 = 0; c2 <= min(50, -51 * c0 + 255); c2 += 1)
      for (int c3 = 0; c3 <= min(50, -51 * c1 + 255); c3 += 1)
        S_3(51 * c0 + c2, 51 * c1 + c3);

What to notice: Graphite extracted the nest as a SCoP (domain \([0, 255]^2\), schedule = the two loops), let isl reschedule it, and tiled the permutable band with GCC's default --param loop-block-tile-size=51 (hence floor(i1/51) and bounds -51 * c0 + 255). permutable: 1 is Definition 18.7.6's "fully permutable"; coincident marks dimensions without loop-carried dependences (parallel). GCC's own interchange pass is disabled here to show Graphite's work alone [GCC-Graphite].

3. Worked example¶

Loop fusion and fission on the running example¶

Take the loop

for (i = 1; i <= n; i++) {
    S1: a[i] = b[i] + 1;
    S2: c[i] = a[i - 1] * 2 + a[i + 1];
    S3: d[i] = d[i - 1] + c[i];
}

Dependences (Lesson 18.6's strong SIV test on each pair): S1 → S2 flow, distance 1 (a[i] written in \(i\), read as a[i-1] in \(i + 1\)); S2 → S1 anti, distance 1 (a[i+1] read in \(i\), written by S1 in \(i + 1\)); S2 → S3 flow, distance 0 (c[i]); S3 → S3 flow, distance 1 (d). The graph has SCCs \(\{S1, S2\}\) (the cycle S1 → S2 → S1) and \(\{S3\}\) (a self-cycle). Algorithm 18.7.2 emits two loops: for i: S1; S2 and for i: S3 — both sequential. Without the a[i + 1] term, the SCCs would be \(\{S1\}, \{S2\}, \{S3\}\): three loops, the first two parallel (vectorizable), which is the point of distribution. Emitting S2's loop before S1's would be wrong: its reads of a[i - 1] would see old values — the topological order is mandatory.

Fusing back for i: S1 and for i: S2 (in the three-loop version): the only dependence from the first to the second has distance \(+1\) (\(i' - i = 1\)), non-negative: fusion legal. With S2 reading a[i + 1] instead, the distance is \(-1\): fusion-preventing — the nofuse box.

Loop interchange on the running example¶

The nest for i, j, k = 1..4: A[i][j][k] = A[i-1][j+1][k-1] + 1 has one flow dependence, distance \((1, -1, 1)\) (the course oracle brute_dependences). Theorem 18.7.4 for all six orders:

order	permuted distance	lexicographically positive?	oracle: same result when executed
`ijk`	\((1, -1, 1)\)	yes	yes
`ikj`	\((1, 1, -1)\)	yes	yes
`jik`	\((-1, 1, 1)\)	no	no
`jki`	\((-1, 1, 1)\)	no	no
`kij`	\((1, 1, -1)\)	yes	yes
`kji`	\((1, -1, 1)\)	yes	yes

The last column comes from run_nest in tools/course/lib/loops.py, which executes the permuted nest on concrete data: it agrees with the theorem in all six cases. Any order that puts \(j\) outermost is illegal, because \(j\) is the only loop whose component is negative.

Strip-mining, skewing and tiling on the running example¶

The stencil of the isl box, A[i][j] = A[i-1][j+1] + A[i][j-1], has distances \((1, -1)\) and \((0, 1)\). The band \((i, j)\) is not fully permutable (\(-1\)). Algorithm 18.7.7: for \(q = j\), \(p = i\): the dependences with \(d_j < 0\), \(d_i > 0\) are \(\{(1, -1)\}\), \(f = \lceil 1/1 \rceil = 1\). Skewing \(j' = j + i\) gives \((1, 0)\) and \((0, 1)\): fully permutable. Tiling \(4 \times 4\) in \((i, j')\): a dependence \(\mathbf{x} \to \mathbf{x} + (1, 0)\) has tile distance \((\lfloor (x_i + 1)/4 \rfloor - \lfloor x_i/4 \rfloor, 0) \in \{(0, 0), (1, 0)\}\), never negative — the isl box's "tiling the skewed (i, i + j) legal: True". Without skewing, \((1, -1)\) from \(\mathbf{x} = (1, 4)\) (tile \((0, 1)\)) to \((2, 3)\) (tile \((0, 0)\)) has tile distance \((0, -1)\): the sink's tile runs first — "tiling (i, j) legal: False".

The polyhedral model on the running example¶

Feautrier's one-dimensional schedule for the same stencil: \(\theta(i, j) = a i + b j + c\). The dependences are uniform, so Farkas reduces to the distance vectors: \(\theta(\mathbf{x} + \mathbf{d}) - \theta(\mathbf{x}) = a d_i + b d_j \ge 1\) for \((0, 1)\) and \((1, -1)\): \(b \ge 1\) and \(a - b \ge 1\). Minimizing latency picks \(b = 1, a = 2\): \(\theta(i, j) = 2i + j\) — the classic wavefront. On \([1, 8]^2\), \(\theta\) ranges over \([3, 24]\): 22 sequential steps instead of 64, and all iterations with the same \(2i + j\) can run in parallel.

Try it

./course drill interchange-legal --seed 3 --difficulty hard --solution shows the six-order table for a 3-deep nest; --difficulty medium includes non-uniform references like A[j][i], whose distances are not constant (the worked solution enumerates them).

4. Invariants and correctness¶

Loop fusion and fission¶

Theorem 18.7.14 (Distribution by SCCs preserves behavior)

The loops emitted by Algorithm 18.7.2 compute the same result as \(L\).

Proof

Take a dependence from instance \(S_p(i)\) to \(S_q(i')\) (source first in \(L\)). If \(S_p\) and \(S_q\) are in the same component, they are in the same emitted loop, which executes that component's statements in the original textual order for every \(i\): the relative order of any two instances is as in \(L\). If they are in different components, the edge \(S_p \to S_q\) exists in \(G\), so \(S_p\)'s component precedes \(S_q\)'s in the topological order; all instances of \(S_p\) then run before all instances of \(S_q\) — in particular the source before the sink. Every dependence is preserved, and independent instances commute; as in Theorem 18.7.4, every read sees the same last write.

Theorem 18.7.15 (Fusion legality)

Fusing \(L_1; L_2\) (same iteration range) preserves behavior if and only if no dependence from an instance of \(L_1\) in iteration \(i\) to an instance of \(L_2\) in iteration \(i'\) has \(i' < i\).

Proof

In the fused loop, \(L_1\)'s body in iteration \(i\) runs before \(L_2\)'s body in iteration \(i'\) iff \(i \le i'\); dependences inside \(L_1\) or inside \(L_2\) keep their order (each body's iterations run in the same order as before). Dependences from \(L_2\) to \(L_1\) cannot exist in the original (all of \(L_1\) ran first). So all dependences are preserved iff every cross dependence has \(i \le i'\). Necessity as in Theorem 18.7.4.

Loop interchange and permutation¶

Theorem 18.7.4 (§2) is the correctness argument; Algorithm 18.7.5 checks its condition for every swap on the direction matrix, which over-approximates the distance vectors (Lesson 18.6), so a swap it accepts is legal.

Strip-mining, skewing and tiling¶

Theorem 18.7.16 (Tiling legality)

Strip-mining alone is always legal. Tiling a band \(k..m\) is legal if the band is fully permutable (Definition 18.7.6).

Proof

Strip-mining maps \(i\) to \((\lfloor i/B \rfloor, i)\), which is order-preserving: \(i < i'\) implies \(\lfloor i/B \rfloor \le \lfloor i'/B \rfloor\) and, if equal, \(i < i'\). So the order of all iterations is unchanged. For tiling, take a dependence \(\mathbf{x} \to \mathbf{x} + \mathbf{d}\) with \(\mathbf{d} \succ \mathbf{0}\). If it is carried outside the band, the outer loops (unchanged) order it. Otherwise its outer components are zero and \(d_k, \dots, d_m \ge 0\). Since \(t \mapsto \lfloor t/B \rfloor\) is monotone, each tile component \(\lfloor (x_q + d_q)/B_q \rfloor - \lfloor x_q/B_q \rfloor \ge 0\). If some tile component is positive, the new position vector difference is lexicographically positive at the tile level. If all are zero, both instances are in the same tile, and the point loops run in the original order of the band, where \(\mathbf{d} \succ \mathbf{0}\). Either way the source runs first; conclude as in Theorem 18.7.4.

The polyhedral model¶

Theorem 18.7.17 (Legality of an affine schedule)

A schedule \(\theta\) preserves behavior if \(\theta_S(\mathbf{x}) \prec \theta_T(\mathbf{y})\) for every pair \((S[\mathbf{x}], T[\mathbf{y}])\) in every dependence relation (for one-dimensional \(\theta\): \(\theta_T(\mathbf{y}) - \theta_S(\mathbf{x}) \ge 1\)), and instances with equal schedule values are independent. The coefficients found by Algorithm 18.7.13 satisfy this.

Proof

The first part is Theorem 18.7.4's argument with \(\pi\) replaced by \(\theta\): every dependent pair keeps its order, independent pairs may be reordered or run in parallel. For the second: Algorithm 18.7.13 requires \(\theta_T(\mathbf{y}) - \theta_S(\mathbf{x}) - 1 \equiv \lambda_0 + \boldsymbol{\lambda}^T(C(\mathbf{x}, \mathbf{y}) + \mathbf{c})\) with non-negative multipliers; on \(\mathcal{P}_e\) the right-hand side is \(\ge 0\) (the "if" direction of Theorem 18.7.12), so \(\theta_T(\mathbf{y}) - \theta_S(\mathbf{x}) \ge 1\) on every dependence. Two instances with equal \(\theta\) therefore cannot be dependent. (Conversely, by the "only if" direction, every legal one-dimensional affine schedule is found this way, which is why the LP is exact for one-dimensional schedules [Fea92].)

5. Complexity¶

\(m\) = statements, \(E\) = dependence edges, \(d\) = nest depth, \(D\) = number of distance vectors.

Technique	Legality test	Transformation	Worst case	Variables
Loop fusion and fission	Fission: SCCs \(O(m + E)\); fusion: one pass over cross dependences \(O(E)\)	\(O(\text{size})\)	Fusion for maximal locality is NP-hard [KM93]	\(m\), \(E\)
Loop interchange and permutation	\(O(D \cdot d)\) per permutation; \(d!\) permutations	\(O(\text{size})\)	Choosing the best order: \(d!\) candidates, pruned by cost	\(D\), \(d\)
Strip-mining, skewing and tiling	\(O(D \cdot d)\) for full permutability; skewing factors \(O(D d^2)\)	New bounds by Fourier–Motzkin	Bounds can blow up doubly exponentially in \(d\)	\(D\), \(d\)
The polyhedral model	Dependences: parametric integer programming	LP / ILP per schedule dimension; code generation	Exponential (integer programming, Fourier–Motzkin in code generation)	constraints, statements

Proposition 18.7.18 (Number of legal permutations)

For a depth-\(d\) nest whose distance vectors are all \(\ge \mathbf{0}\) componentwise (a fully permutable nest), all \(d!\) permutations are legal; a single distance vector with exactly one negative component, in position \(q\), and positive components elsewhere forbids exactly the \((d - 1)!\) permutations that put \(q\) outermost.

Proof

If all components of \(\mathbf{d} \ne \mathbf{0}\) are \(\ge 0\), every permutation of them is \(\ge \mathbf{0}\) componentwise and non-zero, hence lexicographically positive. For the second claim, \(\pi(\mathbf{d})\)'s first component is \(d_{\pi(1)}\), which is non-zero; the vector is lexicographically negative iff \(d_{\pi(1)} < 0\), i.e. \(\pi(1) = q\), and there are \((d - 1)!\) such permutations. (The worked example: \(d = 3\), \(q = j\), the \(2! = 2\) orders jik, jki are illegal.)

Pathological input. Tiling a \(d\)-deep nest with skewed (non-rectangular) bounds: each tile loop's bounds come from projecting the tile constraints, and Fourier–Motzkin elimination can double the number of constraints per eliminated variable — code generators (isl's AST builder, CLooG) spend most of their time here, and polyhedral compilers limit the number of statements or parameters in a SCoP for that reason (Polly's runtime-check limits such as polly-rtc-max-parameters and polly-rtc-max-array-disjuncts; Graphite's --param graphite-max-nb-scop-params, default 10, and graphite-max-arrays-per-scop, default 100).

At scale. Production compilers restrict themselves to the cheap cases: LLVM's interchange handles perfect nests with direction matrices, distribution only isolates LAA-unsafe partitions, and fusion requires adjacent loops with equal trip counts. Full polyhedral optimization is opt-in (-polly, -floop-nest-optimize).

Loop fusion and fission¶

Fusion with peeling or shifting: when trip counts differ by a constant or a dependence has distance \(-1\), peel or shift one loop first (LLVM's loop-fusion-peel-max-count, default 0) — trade-off: extra code for the peeled iterations.
Partial distribution for vectorization (LLVM): only split out what LAA calls unsafe — trade-off: fewer loops, but no parallelism beyond vectorization.
Library-call partitions (GCC memset/memcpy recognition, LLVM's separate loop-idiom pass) — trade-off: only exact idioms.

Loop interchange and permutation¶

Unimodular transformations [WL91]: interchange, reversal and skewing as integer matrices with determinant \(\pm 1\); legality is \(T\mathbf{d} \succ \mathbf{0}\) for all \(\mathbf{d}\) — trade-off: one framework, but only for perfect nests.
Cost models: LLVM's CacheCost (loop cost by cache lines touched) vs GCC's stride model — trade-off: precision vs speed.

Strip-mining, skewing and tiling¶

Multi-level and register tiling (Polly's polly-2nd-level-tiling, polly-register-tiling) — trade-off: better reuse, more complex code.
Parametric tile sizes — trade-off: autotuning without recompiling, but non-affine bounds.
Wavefront parallelism: after skewing, the tiles along an anti-diagonal are independent.

The polyhedral model¶

Pluto's algorithm [BHRS08]: find tiling hyperplanes one by one, minimizing an upper bound of dependence distances — trade-off: good locality and coarse parallelism, integer programming at compile time.
Schedule trees (isl, used by Polly and Graphite): represent bands, filters and sequences explicitly — trade-off: easier composition of transformations.
MLIR's affine dialect: polyhedral-style analysis kept at a higher level than LLVM IR — trade-off: requires front ends to emit it.

7. In real compilers¶

Loop fusion and fission¶

LLVM

llvm/lib/Transforms/Scalar/LoopDistribute.cpp — InstPartitionContainer (partitions; mergeBeforePopulating, mergeAdjacentNonCyclic), LoopDistributeForLoop::processLoop; option enable-loop-distribute (default false) [LLVM-Distribute]. llvm/lib/Transforms/Scalar/LoopFuse.cpp — LoopFuser::fuseCandidates, dependencesAllowFusion [LLVM-Fuse] (LLVM 23.1.2).

GCC gcc/tree-loop-distribution.cc — loop_distribution::build_rdg, classify_partition, distribute_loop (GCC 15) [GCC-Distribution].

The real-world boxes for this technique are in §2.

Loop interchange and permutation¶

LLVM

llvm/lib/Transforms/Scalar/LoopInterchange.cpp — populateDependencyMatrix, isLegalToInterChangeLoops (Theorem 18.7.4 on the direction matrix), LoopInterchangeLegality::canInterchangeLoops, LoopInterchangeProfitability::isProfitable (with CacheCost) [LLVM-Interchange] (LLVM 23.1.2).

GCC gcc/gimple-loop-interchange.cc — tree_loop_interchange::valid_data_dependences, should_interchange_loops, interchange_loops (GCC 15) [GCC-Interchange].

Find where LLVM does it. Open LoopInterchange.cpp and find the function that decides, for one row of the dependence matrix at a time, whether swapping two columns keeps it lexicographically positive. Question: what is it called, and which part of each row does it skip? (Quiz interchange-find-legal.)

The real-world boxes for this technique are in §2.

Strip-mining, skewing and tiling¶

LLVM

No automatic tiling in LLVM proper. #pragma omp tile is implemented by Clang (OMPTileDirective, lowered in clang/lib/CodeGen/CGStmtOpenMP.cpp); Polly tiles in polly/lib/Transform/ScheduleOptimizer.cpp — ScheduleTreeOptimizer::isTileableBandNode, applyTileBandOpt, option polly-default-tile-size (32) [Polly-Docs] (LLVM 23.1.2).

GCC gcc/graphite-optimize-isl.cc — optimize_isl, tile size --param loop-block-tile-size (default 51) (GCC 15) [GCC-Graphite].

The real-world boxes for this technique are in §2.

The polyhedral model¶

LLVM

Polly: polly/lib/Analysis/ScopDetection.cpp and ScopInfo.cpp (SCoP extraction), polly/lib/Analysis/DependenceInfo.cpp (dependences with isl), polly/lib/Transform/ScheduleOptimizer.cpp (isl_schedule_constraints_compute_schedule, tiling), polly/lib/CodeGen/ (isl AST to LLVM IR); architecture in polly/docs/Architecture.rst [GGL12; Polly-Docs] (LLVM 23.1.2).

GCC Graphite: gcc/graphite.cc, graphite-scop-detection.cc, graphite-dependences.cc, graphite-optimize-isl.cc, graphite-isl-ast-to-gimple.cc (GCC 15) [GCC-Graphite].

The real-world boxes for this technique are in §2.

8. Comparison¶

Technique	Power / precision	Speed	Output / error quality	Implementation effort	Typical use
Loop fusion and fission	Separates cycles from vectorizable work; merges loops for reuse	SCCs linear; fusion check per dependence	More (or fewer) loops; remarks on preventing dependences	Medium	Enabling vectorization; locality of producer–consumer loops
Loop interchange and permutation	Any legal order of a perfect nest	\(O(Dd)\) per candidate	Same body, different loop order	Medium (IR surgery on headers/latches)	Unit-stride inner loops; outer-loop parallelism
Strip-mining, skewing and tiling	Cache blocking; wavefront parallelism after skewing	Cheap legality; bounds by projection	Deeper nests with `min`/`max` bounds	Medium to high	Dense linear algebra, stencils
The polyhedral model	Exact dependences and optimal affine schedules for SCoPs	Expensive (integer programming)	Arbitrary affine restructuring in one step	Very high (isl, code generation)	HPC kernels; Polly, Graphite, MLIR affine

Choose distribution to isolate a recurrence from vectorizable statements; choose fusion for producer–consumer loops over the same range. Choose interchange when the inner loop has a large stride. Choose tiling for nests that reuse data across an outer loop (matrix multiply, stencils) — after skewing if needed. Choose the polyhedral model when the kernel is a SCoP and worth compile time; otherwise the individual transformations above cover the common cases.

9. Assessment¶

Technique	Quiz ids (`solutions/quizzes/ch18.yaml`)	Drill	Flashcard tag	Exercises
Loop fusion and fission	`fission-order`, `fusion-legal`	`./course drill dependence-test` (the distances that decide fusion)	`fusion-fission`	—
Loop interchange and permutation	`interchange-legal-orders`, `interchange-direction`, `interchange-find-legal`	`./course drill interchange-legal`	`interchange`	—
Strip-mining, skewing and tiling	`tiling-skew-factor`, `tiling-permutable`	`./course drill interchange-legal` (full permutability is "every order legal")	`tiling`	—
The polyhedral model	`poly-wavefront`, `poly-dataflow`	none: the schedule LP is computed with isl in the real-world box; the drill covers its legality condition	`polyhedral`	—

Pitfall

"No dependence is carried by the inner loop, so I can interchange" is wrong: legality depends on the permuted vectors. The distance \((1, -1)\) is carried by the outer loop, yet interchange makes it \((-1, 1)\). Always permute the vectors and re-check lexicographic positivity (Theorem 18.7.4).

References¶

See the chapter references.

Lesson 18.7 — Dependence-driven transformations: fusion, fission, interchange, tiling, the polyhedral model¶

1. Problem and motivation¶

Loop fusion and fission¶

Loop interchange and permutation¶

Strip-mining, skewing and tiling¶

The polyhedral model¶

2. Definitions and algorithms¶

Loop fusion and fission¶

Loop interchange and permutation¶

Strip-mining, skewing and tiling¶

The polyhedral model¶

3. Worked example¶

Loop fusion and fission on the running example¶

Loop interchange on the running example¶

Strip-mining, skewing and tiling on the running example¶

The polyhedral model on the running example¶

4. Invariants and correctness¶

Loop fusion and fission¶

Loop interchange and permutation¶

Strip-mining, skewing and tiling¶

The polyhedral model¶

5. Complexity¶

6. Variants and refinements¶

Loop fusion and fission¶

Loop interchange and permutation¶

Strip-mining, skewing and tiling¶

The polyhedral model¶

7. In real compilers¶

Loop fusion and fission¶

Loop interchange and permutation¶

Strip-mining, skewing and tiling¶

The polyhedral model¶

8. Comparison¶

9. Assessment¶

References¶