Cache-Aware Blocking / Tiling

Why Cache-Aware Blocking / Tiling Matters

Two programs can compute exactly the same matrix product, perform the same number of multiplications, and move different amounts of data through memory. The difference is often locality: whether a value is still nearby when it is needed again. A cache fetches a whole line of neighboring words, so using a small region repeatedly can be much cheaper than returning to it after its lines have been evicted.

Cache-aware blocking, also called tiling, divides a loop iteration space into small pieces chosen with the memory hierarchy in mind. Finish useful work on one piece while its data is resident, then move on. The transformation changes traversal and reuse; it does not automatically change the mathematical task or reduce its arithmetic complexity.

This lesson uses dense matrix multiplication as a concrete example. The Python, C++, and Java programs compare an untiled row-major traversal with a tiled traversal through the same fully associative LRU cache simulator. Both calculate every product term, preserve partial sums, handle rectangular edge tiles, and count line fills plus dirty writebacks.

You will trace every tile of a small product, relate tile dimensions to cache capacity, and distinguish word references from block transfers. The same reasoning applies to transpose, image blocks, stencil neighborhoods, and suitable dynamic-programming loops, provided their dependencies permit the reordering. Page-oriented locality from the earlier external-memory lessons becomes cache-line locality here.

Core Concepts and Terminology

Let $A$ have shape $r\times k$, $B$ shape $k\times c$, and $C$ shape $r\times c$. The task is $$ C_{ij}=\sum_{x=0}^{k-1}A_{ix}B_{xj}, \qquad 0\le i\lt r,\quad 0\le j\lt c. $$ Matrices are flat row-major arrays: $A_{ix}$ is at offset $ik+x$, $B_{xj}$ at offset $xc+j$, and $C_{ij}$ at offset $ic+j$. The innermost $j$ loop walks contiguous row segments of $B$ and $C$. The scalar $A_{ix}$ is loaded once and held across that segment.

A tile is a bounded group of iterations, not a unit of storage transfer. Tile width $T$ creates row starts ii, column starts jj, and inner-dimension starts kk. One full cubic tile product updates a $T\times T$ output block using two $T\times T$ input blocks. With scalar word width $w=8$ logical bytes, its useful data is about $3T^2w$ bytes.

The cache model has $\lambda$ words per line and $Q$ line frames, so its data capacity is $M=Q\lambda$ words or $8Q\lambda$ bytes. Address $a$ maps to line $\lfloor a/\lambda\rfloor$ and slot $a\bmod\lambda$. The model is fully associative with LRU replacement: any line can occupy any frame, and the least recently referenced resident line is evicted when space is needed.

Write policy matters. A missing load or store fetches one full line: write allocate applies to stores. A store updates the cached line and marks it dirty. A clean eviction has no writeback; a dirty eviction writes a full line to backing memory. A final flush writes all remaining dirty lines, then marks them clean. Thus $$ \text{transfers}=\text{line reads}+\text{dirty eviction writes}+\text{final flush writes}. $$ Word references and cache hits are different counts. Every load/store is one reference, but only a miss causes a line read. The kernel reads $C_{ij}$ immediately before storing its update, so its stores hit; the general simulator still implements store-miss allocation for other access patterns.

Address layout and measurement boundary. $A$ begins at word address zero. $B$ starts at the next line boundary after $A$, and $C$ at the next line boundary after $B$. Padding is zero and is never treated as a matrix entry. Input loading, initial zero-filled $C$, and final inspection of flushed backing memory are outside the cache-event count. Both compared runs start cold with identical input images and include their own final dirty flush.

The simple sizing rule $3T^2\lesssim M$ leaves room for the useful entries, but actual line footprint depends on alignment and row stride. A two-word row segment can span two lines if misaligned. Leave headroom and inspect the distinct lines touched; fitting a count of words alone is not a universal performance guarantee. The worked $T=2$, $\lambda=2$, $Q=7$ case has fourteen cache words and twelve useful tile words.

How Cache-Aware Blocking / Tiling Works, Step by Step

1. Fix layout and a baseline. Both programs compute $C=AB$ from a zero output using an inner i,k,j kernel. This already makes $j$ contiguous in row-major $B$ and $C$. The comparison isolates tiling of that kernel; it does not compare against a deliberately strided i,j,k baseline or a tuned BLAS implementation.

2. Choose a tile size. Use $T$ as a candidate, estimate the three-block working set, and account for line fragments and other live data. Cache-aware code chooses/tunes $T$ using a target cache or measured behavior. The simulator accepts a tile even if it exceeds cache capacity: that remains mathematically valid and is useful for demonstrating poor choices.

3. Partition all three dimensions. Loop ii over output row blocks, jj over output column blocks, then kk over the reduction/inner blocks. Ends are clamped independently with min(start + T, dimension). Nonmultiple and rectangular shapes require this in every dimension.

4. Reuse within a tile. For each $i$ in the row block and each inner index $x$ in the current reduction block, load $A_{ix}$ once. For each $j$ in the column block, load $B_{xj}$, load the current $C_{ij}$, and store $C_{ij}+A_{ix}B_{xj}$. Cached words remain available until replacement evicts their lines. There is no explicit copy into packed tile arrays in these programs.

5. Preserve the output between inner tiles. $C$ is zeroed once during setup. The next kk block adds another partial dot product to the same output region. Resetting $C$ inside the kk loop would discard earlier terms. The next figure follows all row-pair updates of the first tile product.

6. Move to the next output tile. Only after processing all reduction blocks for a given ii,jj pair does the outer traversal select another output tile. The loops still compute every $(i,x,j)$ product exactly once.

7. Materialize dirty results fairly. After the kernel, flush every remaining dirty cache line and include those writes in the transfer count. Compare complete output arrays only after this flush. Output inspection is an uninstrumented observation of backing memory, not a counted second cache pass. Empty dimensions perform no kernel accesses; their output zeros already exist in the initialized image.

flowchart TD
    A["Validate shapes and cache parameters"] --> B["Align A, B, C regions; zero C; start cold cache"]
    B --> C{"All dimensions positive?"}
    C -->|No| P["Flush remaining dirty lines and return result"]
    C -->|Yes| D["Choose next ii,jj output tile"]
    D --> E["Choose next kk reduction tile and clamp all ends"]
    E --> F["Choose i,x; load A scalar"]
    F --> G["Choose j; load B and C; store accumulated C"]
    G --> H{"More j in tile?"}
    H -->|Yes| G
    H -->|No| I{"More i,x in tile?"}
    I -->|Yes| F
    I -->|No| J{"More kk blocks?"}
    J -->|Yes| E
    J -->|No| K{"More ii,jj blocks?"}
    K -->|Yes| D
    K -->|No| P
    G -.-> L["Each reference updates LRU; a miss reads one line"]
    L -.-> M["A dirty eviction writes one line"]

Worked Example with a Full Trace

Set $r=k=c=4$, $T=2$, $\lambda=2$ words per line, and $Q=7$ lines. Each line is $16$ bytes; cache capacity is $112$ bytes. Use $$ A=\begin{bmatrix} 1&2&3&4\\5&6&7&8\\9&10&11&12\\13&14&15&16 \end{bmatrix}, \qquad B=\begin{bmatrix} 1&2&3&4\\0&1&0&1\\2&0&1&0\\1&1&1&1 \end{bmatrix}. $$ $A$ occupies word addresses $0$ through $15$ and lines $0$ through $7$. $B$ occupies addresses $16$ through $31$, lines $8$ through $15$. Zero-initialized $C$ occupies addresses $32$ through $47$, lines $16$ through $23$. Each row segment of two aligned words is exactly one line.

First tile's complete cache-reference trace. For ii=jj=kk=0, a reference sequence uses L for load and S for store, followed by the global line ID. M means miss/read, H means hit. The four row-pair updates each make seven references:

1. $i=0,x=0$: L0 M, L8 M, L16 M, S16 H, L8 H, L16 H, S16 H. Output row pair becomes $[1,2]$.
2. $i=0,x=1$: L0 H, L10 M, L16 H, S16 H, L10 H, L16 H, S16 H. Pair becomes $[1,4]$.
3. $i=1,x=0$: L2 M, L8 H, L18 M, S18 H, L8 H, L18 H, S18 H. Pair becomes $[5,10]$.
4. $i=1,x=1$: L2 H, L10 H, L18 H, S18 H, L10 H, L18 H, S18 H. Pair becomes $[5,16]$.

This tile has twenty-eight references, six line reads, twenty-two hits, and zero writes so far. Lines $16$ and $18$ are dirty; the seven-line cache has only six residents at this point. LRU order is $[0,16,8,2,10,18]$, oldest first.

All eight tile products. The following list records every row-pair arithmetic update, the complete affected C tile afterward, and all cumulative cache state. Each row-pair result gives both column values after its two consecutive $j$ updates. Together they account for all sixty-four multiply-adds. Unmentioned output tiles remain unchanged. Dirty lists contain global line IDs; counters are before the final flush.

1. Tile start (0,0,0). i=0, x=0 → $[1,2]$; i=0, x=1 → $[1,4]$; i=1, x=0 → $[5,10]$; i=1, x=1 → $[5,16]$. C tile = $\begin{bmatrix}1&4\\5&16\end{bmatrix}$. Cumulative references $28$, hits $22$, line reads $6$, eviction writes $0$. LRU = [0, 16, 8, 2, 10, 18]; dirty = [16, 18].

1. Tile start (0,0,2). i=0, x=2 → $[7,4]$; i=0, x=3 → $[11,8]$; i=1, x=2 → $[19,16]$; i=1, x=3 → $[27,24]$. C tile = $\begin{bmatrix}11&8\\27&24\end{bmatrix}$. Cumulative references $56$, hits $46$, line reads $10$, eviction writes $0$. LRU = [10, 1, 16, 12, 3, 14, 18]; dirty = [16, 18].

1. Tile start (0,2,0). i=0, x=0 → $[3,4]$; i=0, x=1 → $[3,6]$; i=1, x=0 → $[15,20]$; i=1, x=1 → $[15,26]$. C tile = $\begin{bmatrix}3&6\\15&26\end{bmatrix}$. Cumulative references $84$, hits $68$, line reads $16$, eviction writes $1$. LRU = [18, 0, 17, 9, 2, 11, 19]; dirty = [18, 17, 19].

1. Tile start (0,2,2). i=0, x=2 → $[6,6]$; i=0, x=3 → $[10,10]$; i=1, x=2 → $[22,26]$; i=1, x=3 → $[30,34]$. C tile = $\begin{bmatrix}10&10\\30&34\end{bmatrix}$. Cumulative references $112$, hits $92$, line reads $20$, eviction writes $2$. LRU = [11, 1, 17, 13, 3, 15, 19]; dirty = [17, 19].

1. Tile start (2,0,0). i=2, x=0 → $[9,18]$; i=2, x=1 → $[9,28]$; i=3, x=0 → $[13,26]$; i=3, x=1 → $[13,40]$. C tile = $\begin{bmatrix}9&28\\13&40\end{bmatrix}$. Cumulative references $140$, hits $114$, line reads $26$, eviction writes $3$. LRU = [19, 4, 20, 8, 6, 10, 22]; dirty = [19, 20, 22].

1. Tile start (2,0,2). i=2, x=2 → $[31,28]$; i=2, x=3 → $[43,40]$; i=3, x=2 → $[43,40]$; i=3, x=3 → $[59,56]$. C tile = $\begin{bmatrix}43&40\\59&56\end{bmatrix}$. Cumulative references $168$, hits $138$, line reads $30$, eviction writes $4$. LRU = [10, 5, 20, 12, 7, 14, 22]; dirty = [20, 22].

1. Tile start (2,2,0). i=2, x=0 → $[27,36]$; i=2, x=1 → $[27,46]$; i=3, x=0 → $[39,52]$; i=3, x=1 → $[39,66]$. C tile = $\begin{bmatrix}27&46\\39&66\end{bmatrix}$. Cumulative references $196$, hits $160$, line reads $36$, eviction writes $5$. LRU = [22, 4, 21, 9, 6, 11, 23]; dirty = [22, 21, 23].

1. Tile start (2,2,2). i=2, x=2 → $[38,46]$; i=2, x=3 → $[50,58]$; i=3, x=2 → $[54,66]$; i=3, x=3 → $[70,82]$. C tile = $\begin{bmatrix}50&58\\70&82\end{bmatrix}$. Cumulative references $224$, hits $184$, line reads $40$, eviction writes $6$. LRU = [11, 5, 21, 13, 7, 15, 23]; dirty = [21, 23].

The final matrix is $$ C=\begin{bmatrix} 11&8&10&10\\27&24&30&34\\43&40&50&58\\59&56&70&82 \end{bmatrix}. $$ Before flush, tiled execution has forty line reads and six dirty eviction writes. Resident dirty lines $21$ and $23$ contain the last output-column segments; flushing them adds two writes. Total tiled transfers are $40+6+2=48$ lines, or $768$ bytes. Its $224$ references equal $184$ hits plus $40$ misses.

The untiled i,k,j traversal produces the identical matrix with $208$ references, $160$ hits, $48$ line reads, six eviction writes, and two flush writes: $56$ transfers, or $896$ bytes. Tiling therefore saves eight line transfers here while making sixteen additional A references. It has not reduced the sixty-four multiplications or the sixty-four additions, and the simulator does not claim an elapsed-time speedup.

Implementation in Python, C++, and Java

Save Python as main.py and run python3 main.py. Save C++ as main.cpp and compile with g++ -std=c++17 -O2 main.cpp -o tiled-multiply. Save Java as Main.java and run javac Main.java followed by java Main. Python 3.7 or newer supplies dataclasses. C++ uses long long and Java uses long for integer words; all examples use exact integer arithmetic within the declared bounds.

The API is multiply(a, b, rows, inner, cols, tile, line_words, frames, blocked), with corresponding camel-case parameter names in Java. Pass flat row-major Python lists, C++ vector<Word>, or Java long[]. Null/nonsequence inputs are outside the contract. blocked=false runs the untiled i,k,j kernel; true uses ii,jj,kk,i,k,j traversal. Python defaults to blocked mode; C++ also supplies that default, while Java calls pass the boolean explicitly.

Dimensions must be integers from zero through $64$. Tile size and line width are each $1$ through $64$, and cache line capacity is $1$ through $128$. Array lengths must equal rows*inner and inner*cols, and entries must be integers in $[-1000,1000]$. With at most sixty-four terms, the magnitude of any complete or partial output sum is at most $64\cdot10^6$, so these test-domain products fit the 64-bit accumulators without overflow.

Out-of-range dimensions raise DIMENSION, invalid tile/cache limits raise PARAMETER, wrong array lengths raise SHAPE, and invalid entries raise VALUE, in that order. Python first checks exact parameter types and an exact boolean mode, rejecting booleans/floats used as integer dimensions, parameters, or entries. Both inputs remain unchanged. A zero dimension performs no kernel accesses; if the inner dimension is zero, output is the zero matrix already present in the initialized backing image.

Cache mechanics. The cache stores actual copied line values and dirty flags, with at most frames resident lines. Every modeled word load/store goes through fetch. Eviction writes dirty values back before discarding a line; final flush writes all remaining dirty lines and clears their dirty flags. A clean line is never written merely because it is evicted. C++ uses an LRU list plus hash map; Python uses OrderedDict; Java uses access-ordered LinkedHashMap.

Results contain the flat output values and references, hits, line reads, eviction writebacks, final-flush writebacks, evictions, and peak resident lines. Transfers are the sum of reads and both write categories. The final array is observed directly from flushed backing memory without creating more cache events. Setup/validation and output inspection are outside measurement; the matrix kernel and final dirty flush are inside it.

The simulator itself uses host arrays for the entire backing memory and output observations. Only resident line payloads obey the modeled $Q\lambda$-word capacity; metadata, scalar registers, inputs, backing arrays, and returned output use additional host memory. Python's objects are not physically eight-byte integers here. The printed byte counts describe the logical model, not measured machine traffic. Compilers, actual caches, prefetchers, SIMD, and allocation are not being benchmarked by these instrumented programs.

from dataclasses import dataclass
from collections import OrderedDict


@dataclass
class Line:
    values: list
    dirty: bool = False


@dataclass
class Result:
    values: list
    accesses: int
    hits: int
    reads: int
    eviction_writes: int
    flush_writes: int
    evictions: int
    peak_lines: int


class Cache:
    def __init__(self, memory, line_words, frames):
        self.memory, self.line_words, self.frames = memory, line_words, frames
        self.lines = OrderedDict()
        self.accesses = self.hits = self.reads = 0
        self.eviction_writes = self.flush_writes = self.evictions = self.peak_lines = 0

    def writeback(self, line_id, line):
        start = line_id * self.line_words
        self.memory[start:start + self.line_words] = line.values
        line.dirty = False

    def fetch(self, address):
        assert 0 <= address < len(self.memory)
        self.accesses += 1
        line_id = address // self.line_words
        if line_id in self.lines:
            self.hits += 1
            self.lines.move_to_end(line_id)
            return self.lines[line_id]
        self.reads += 1
        if len(self.lines) == self.frames:
            old_id, old = self.lines.popitem(last=False)
            self.evictions += 1
            if old.dirty:
                self.writeback(old_id, old)
                self.eviction_writes += 1
        start = line_id * self.line_words
        line = Line(self.memory[start:start + self.line_words])
        self.lines[line_id] = line
        self.peak_lines = max(self.peak_lines, len(self.lines))
        return line

    def load(self, address):
        return self.fetch(address).values[address % self.line_words]

    def store(self, address, value):
        line = self.fetch(address)
        line.values[address % self.line_words] = value
        line.dirty = True

    def flush(self):
        for line_id, line in self.lines.items():
            if line.dirty:
                self.writeback(line_id, line)
                self.flush_writes += 1


def multiply(a, b, rows, inner, cols, tile, line_words, frames, blocked=True):
    params = (rows, inner, cols, tile, line_words, frames)
    if any(type(x) is not int for x in params) or type(blocked) is not bool:
        raise ValueError("PARAMETER")
    if not all(0 <= x <= 64 for x in (rows, inner, cols)):
        raise ValueError("DIMENSION")
    if not (1 <= tile <= 64 and 1 <= line_words <= 64 and 1 <= frames <= 128):
        raise ValueError("PARAMETER")
    if len(a) != rows * inner or len(b) != inner * cols:
        raise ValueError("SHAPE")
    if any(type(x) is not int or not -1000 <= x <= 1000 for x in a + b):
        raise ValueError("VALUE")
    def padded(n):
        return (n + line_words - 1) // line_words * line_words
    b_base = padded(len(a))
    c_base = b_base + padded(len(b))
    memory = [0] * (c_base + padded(rows * cols))
    memory[:len(a)] = a
    memory[b_base:b_base + len(b)] = b
    cache = Cache(memory, line_words, frames)

    def kernel(i0, i1, k0, k1, j0, j1):
        for i in range(i0, i1):
            for k in range(k0, k1):
                value_a = cache.load(i * inner + k)
                for j in range(j0, j1):
                    value_b = cache.load(b_base + k * cols + j)
                    address_c = c_base + i * cols + j
                    old_c = cache.load(address_c)
                    cache.store(address_c, old_c + value_a * value_b)
    if rows and inner and cols:
        if blocked:
            for ii in range(0, rows, tile):
                for jj in range(0, cols, tile):
                    for kk in range(0, inner, tile):
                        kernel(ii, min(ii + tile, rows), kk, min(kk + tile, inner),
                               jj, min(jj + tile, cols))
        else:
            kernel(0, rows, 0, inner, 0, cols)
    cache.flush()
    return Result(memory[c_base:c_base + rows * cols], cache.accesses, cache.hits,
                  cache.reads, cache.eviction_writes, cache.flush_writes,
                  cache.evictions, cache.peak_lines)


def show(label, result):
    writes = result.eviction_writes + result.flush_writes
    print(f"{label} accesses={result.accesses} hits={result.hits} reads={result.reads} "
          f"eviction_writes={result.eviction_writes} flush_writes={result.flush_writes} "
          f"transfers={result.reads + writes} evictions={result.evictions} "
          f"peak_lines={result.peak_lines}")


def main():
    a = list(range(1, 17))
    b = [1,2,3,4, 0,1,0,1, 2,0,1,0, 1,1,1,1]
    plain = multiply(a, b, 4, 4, 4, 2, 2, 7, False)
    blocked = multiply(a, b, 4, 4, 4, 2, 2, 7, True)
    assert plain.values == blocked.values
    print("shape=4x4 * 4x4 tile=2 line_words=2 cache_lines=7")
    print("line_bytes=16 cache_bytes=112 equal=true")
    print("C=[" + ",".join(map(str, blocked.values)) + "]")
    show("untiled", plain)
    show("tiled", blocked)


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cassert>
#include <iostream>
#include <list>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>
using namespace std;
using Word = long long;

struct Result {
    vector<Word> values;
    int accesses, hits, reads, eviction_writes, flush_writes, evictions, peak_lines;
};
struct Line {
    vector<Word> values;
    bool dirty;
    list<int>::iterator position;
};
struct Cache {
    vector<Word>& memory;
    int line_words, frames;
    list<int> lru;
    unordered_map<int,Line> lines;
    int accesses = 0, hits = 0, reads = 0;
    int eviction_writes = 0, flush_writes = 0, evictions = 0, peak_lines = 0;
    Cache(vector<Word>& m, int l, int q) : memory(m), line_words(l), frames(q) {}
    Cache(const Cache&) = delete;
    Cache& operator=(const Cache&) = delete;
    void writeback(int id, Line& line) {
        copy(line.values.begin(), line.values.end(), memory.begin() + id * line_words);
        line.dirty = false;
    }
    Line& fetch(int address) {
        assert(address >= 0 && address < static_cast<int>(memory.size()));
        ++accesses;
        int id = address / line_words;
        auto found = lines.find(id);
        if (found != lines.end()) {
            ++hits; lru.splice(lru.end(), lru, found->second.position); return found->second;
        }
        ++reads;
        if (static_cast<int>(lines.size()) == frames) {
            int old_id = lru.front(); auto old = lines.find(old_id);
            ++evictions;
            if (old->second.dirty) { writeback(old_id, old->second); ++eviction_writes; }
            lines.erase(old); lru.pop_front();
        }
        int start = id * line_words;
        vector<Word> values(memory.begin() + start, memory.begin() + start + line_words);
        lru.push_back(id);
        auto inserted = lines.emplace(id, Line{move(values), false, prev(lru.end())});
        peak_lines = max(peak_lines, static_cast<int>(lines.size()));
        return inserted.first->second;
    }
    Word load(int address) { return fetch(address).values[address % line_words]; }
    void store(int address, Word value) {
        Line& line = fetch(address); line.values[address % line_words] = value; line.dirty = true;
    }
    void flush() {
        for (auto& entry : lines)
            if (entry.second.dirty) { writeback(entry.first, entry.second); ++flush_writes; }
    }
};
Result multiply(const vector<Word>& a, const vector<Word>& b, int rows, int inner, int cols,
                int tile, int line_words, int frames, bool blocked = true) {
    if (rows < 0 || rows > 64 || inner < 0 || inner > 64 || cols < 0 || cols > 64)
        throw invalid_argument("DIMENSION");
    if (tile < 1 || tile > 64 || line_words < 1 || line_words > 64 || frames < 1 || frames > 128)
        throw invalid_argument("PARAMETER");
    if (a.size() != static_cast<size_t>(rows * inner) || b.size() != static_cast<size_t>(inner * cols))
        throw invalid_argument("SHAPE");
    for (Word x : a) if (x < -1000 || x > 1000) throw invalid_argument("VALUE");
    for (Word x : b) if (x < -1000 || x > 1000) throw invalid_argument("VALUE");
    auto padded = [line_words](int n) { return (n + line_words - 1) / line_words * line_words; };
    int b_base = padded(static_cast<int>(a.size()));
    int c_base = b_base + padded(static_cast<int>(b.size()));
    vector<Word> memory(c_base + padded(rows * cols), 0);
    copy(a.begin(), a.end(), memory.begin()); copy(b.begin(), b.end(), memory.begin() + b_base);
    Cache cache(memory, line_words, frames);
    auto kernel = [&](int i0, int i1, int k0, int k1, int j0, int j1) {
        for (int i = i0; i < i1; ++i) for (int k = k0; k < k1; ++k) {
            Word value_a = cache.load(i * inner + k);
            for (int j = j0; j < j1; ++j) {
                Word value_b = cache.load(b_base + k * cols + j);
                int address_c = c_base + i * cols + j;
                Word old_c = cache.load(address_c);
                cache.store(address_c, old_c + value_a * value_b);
            }
        }
    };
    if (rows && inner && cols) {
        if (blocked) {
            for (int ii = 0; ii < rows; ii += tile)
                for (int jj = 0; jj < cols; jj += tile)
                    for (int kk = 0; kk < inner; kk += tile)
                        kernel(ii, min(ii + tile, rows), kk, min(kk + tile, inner),
                               jj, min(jj + tile, cols));
        } else kernel(0, rows, 0, inner, 0, cols);
    }
    cache.flush();
    vector<Word> values(memory.begin() + c_base, memory.begin() + c_base + rows * cols);
    return Result{move(values), cache.accesses, cache.hits, cache.reads, cache.eviction_writes,
                  cache.flush_writes, cache.evictions, cache.peak_lines};
}
void show(const string& label, const Result& r) {
    cout << label << " accesses=" << r.accesses << " hits=" << r.hits << " reads=" << r.reads
         << " eviction_writes=" << r.eviction_writes << " flush_writes=" << r.flush_writes
         << " transfers=" << r.reads + r.eviction_writes + r.flush_writes
         << " evictions=" << r.evictions << " peak_lines=" << r.peak_lines << '\n';
}
int main() {
    vector<Word> a;
    for (int i = 1; i <= 16; ++i) a.push_back(i);
    vector<Word> b = {1,2,3,4,0,1,0,1,2,0,1,0,1,1,1,1};
    Result plain = multiply(a,b,4,4,4,2,2,7,false);
    Result blocked = multiply(a,b,4,4,4,2,2,7,true);
    assert(plain.values == blocked.values);
    cout << "shape=4x4 * 4x4 tile=2 line_words=2 cache_lines=7\n"
         << "line_bytes=16 cache_bytes=112 equal=true\nC=[";
    for (size_t i = 0; i < blocked.values.size(); ++i) { if (i) cout << ','; cout << blocked.values[i]; }
    cout << "]\n"; show("untiled",plain); show("tiled",blocked); return 0;
}

import java.util.Arrays;
import java.util.LinkedHashMap;
import java.util.Map;

public class Main {
    static final class Result {
        final long[] values;
        final int accesses, hits, reads, evictionWrites, flushWrites, evictions, peakLines;
        Result(long[] v, Cache c) {
            values = v; accesses = c.accesses; hits = c.hits; reads = c.reads;
            evictionWrites = c.evictionWrites; flushWrites = c.flushWrites;
            evictions = c.evictions; peakLines = c.peakLines;
        }
    }
    static final class Line {
        final long[] values;
        boolean dirty;
        Line(long[] values) { this.values = values; }
    }
    static final class Cache {
        final long[] memory;
        final int lineWords, frames;
        final LinkedHashMap<Integer,Line> lines = new LinkedHashMap<>(16,0.75f,true);
        int accesses, hits, reads, evictionWrites, flushWrites, evictions, peakLines;
        Cache(long[] memory, int lineWords, int frames) {
            this.memory = memory; this.lineWords = lineWords; this.frames = frames;
        }
        void writeback(int id, Line line) {
            System.arraycopy(line.values,0,memory,id*lineWords,lineWords); line.dirty = false;
        }
        Line fetch(int address) {
            if (address < 0 || address >= memory.length) throw new AssertionError("address");
            ++accesses;
            int id = address / lineWords; Line line = lines.get(id);
            if (line != null) { ++hits; return line; }
            ++reads;
            if (lines.size() == frames) {
                Map.Entry<Integer,Line> old = lines.entrySet().iterator().next();
                ++evictions;
                if (old.getValue().dirty) { writeback(old.getKey(),old.getValue()); ++evictionWrites; }
                lines.remove(old.getKey());
            }
            int start = id*lineWords;
            line = new Line(Arrays.copyOfRange(memory,start,start+lineWords));
            lines.put(id,line); peakLines = Math.max(peakLines,lines.size()); return line;
        }
        long load(int address) { return fetch(address).values[address%lineWords]; }
        void store(int address, long value) {
            Line line = fetch(address); line.values[address%lineWords] = value; line.dirty = true;
        }
        void flush() {
            for (Map.Entry<Integer,Line> e : lines.entrySet())
                if (e.getValue().dirty) { writeback(e.getKey(),e.getValue()); ++flushWrites; }
        }
    }
    static int padded(int n, int lineWords) { return (n+lineWords-1)/lineWords*lineWords; }
    static void kernel(Cache cache, int inner, int cols, int bBase, int cBase,
                       int i0, int i1, int k0, int k1, int j0, int j1) {
        for (int i = i0; i < i1; ++i) for (int k = k0; k < k1; ++k) {
            long valueA = cache.load(i*inner+k);
            for (int j = j0; j < j1; ++j) {
                long valueB = cache.load(bBase+k*cols+j);
                int addressC = cBase+i*cols+j;
                long oldC = cache.load(addressC);
                cache.store(addressC,oldC+valueA*valueB);
            }
        }
    }
    static Result multiply(long[] a, long[] b, int rows, int inner, int cols,
                           int tile, int lineWords, int frames, boolean blocked) {
        if (rows < 0 || rows > 64 || inner < 0 || inner > 64 || cols < 0 || cols > 64)
            throw new IllegalArgumentException("DIMENSION");
        if (tile < 1 || tile > 64 || lineWords < 1 || lineWords > 64 || frames < 1 || frames > 128)
            throw new IllegalArgumentException("PARAMETER");
        if (a.length != rows*inner || b.length != inner*cols) throw new IllegalArgumentException("SHAPE");
        for (long x : a) if (x < -1000 || x > 1000) throw new IllegalArgumentException("VALUE");
        for (long x : b) if (x < -1000 || x > 1000) throw new IllegalArgumentException("VALUE");
        int bBase = padded(a.length,lineWords), cBase = bBase+padded(b.length,lineWords);
        long[] memory = new long[cBase+padded(rows*cols,lineWords)];
        System.arraycopy(a,0,memory,0,a.length); System.arraycopy(b,0,memory,bBase,b.length);
        Cache cache = new Cache(memory,lineWords,frames);
        if (rows != 0 && inner != 0 && cols != 0) {
            if (blocked) {
                for (int ii = 0; ii < rows; ii += tile)
                    for (int jj = 0; jj < cols; jj += tile)
                        for (int kk = 0; kk < inner; kk += tile)
                            kernel(cache,inner,cols,bBase,cBase,ii,Math.min(ii+tile,rows),
                                   kk,Math.min(kk+tile,inner),jj,Math.min(jj+tile,cols));
            } else kernel(cache,inner,cols,bBase,cBase,0,rows,0,inner,0,cols);
        }
        cache.flush();
        return new Result(Arrays.copyOfRange(memory,cBase,cBase+rows*cols),cache);
    }
    static void show(String label, Result r) {
        System.out.println(label+" accesses="+r.accesses+" hits="+r.hits+" reads="+r.reads
            +" eviction_writes="+r.evictionWrites+" flush_writes="+r.flushWrites
            +" transfers="+(r.reads+r.evictionWrites+r.flushWrites)
            +" evictions="+r.evictions+" peak_lines="+r.peakLines);
    }
    public static void main(String[]args) {
        long[] a = new long[16]; for (int i = 0; i < a.length; ++i) a[i] = i+1;
        long[] b = {1,2,3,4,0,1,0,1,2,0,1,0,1,1,1,1};
        Result plain = multiply(a,b,4,4,4,2,2,7,false);
        Result blocked = multiply(a,b,4,4,4,2,2,7,true);
        if (!Arrays.equals(plain.values,blocked.values)) throw new AssertionError("product");
        System.out.println("shape=4x4 * 4x4 tile=2 line_words=2 cache_lines=7");
        System.out.println("line_bytes=16 cache_bytes=112 equal=true");
        StringBuilder s = new StringBuilder("C=[");
        for (int i = 0; i < blocked.values.length; ++i) { if (i != 0) s.append(','); s.append(blocked.values[i]); }
        System.out.println(s.append(']')); show("untiled",plain); show("tiled",blocked);
    }
}

Expected standard output, identical in all three languages:

• shape=4x4 * 4x4 tile=2 line_words=2 cache_lines=7
• line_bytes=16 cache_bytes=112 equal=true
• C=[11,8,10,10,27,24,30,34,43,40,50,58,59,56,70,82]
• untiled accesses=208 hits=160 reads=48 eviction_writes=6 flush_writes=2 transfers=56 evictions=41 peak_lines=7
• tiled accesses=224 hits=184 reads=40 eviction_writes=6 flush_writes=2 transfers=48 evictions=33 peak_lines=7

Complexity, Correctness, and Invariants

Arithmetic and references. For positive dimensions, both kernels perform exactly $rkc$ multiplications and $rkc$ additions, so ordinary arithmetic time is $\Theta(rkc)$. The tiled traversal adds block-loop work and may reload an A scalar once per output-column tile. Its exact modeled word-reference count is $$ U_{\mathrm{untiled}}=rk+3rkc, \qquad U_{\mathrm{tiled}}=rk\left\lceil\frac{c}{T}\right\rceil+3rkc. $$ The factor three counts one B load, one C load, and one C store per product. These formulas assume every dimension is positive; otherwise the guarded kernels make zero references. For the sample, the counts are $16+192=208$ and $32+192=224$. A reference is not necessarily a transfer.

Exact transfer accounting. If $R$ is line reads, $W_e$ dirty eviction writes, and $W_f$ final flush writes, measured transfers are $R+W_e+W_f$, with $(R+W_e+W_f)8\lambda$ logical bytes. Also $U=\text{hits}+R$. These counts depend on dimensions, tile order, alignment, line width, and cache capacity. Since there is no data-dependent sparsity shortcut, they are identical for best/average/worst numeric values at fixed geometry and parameters. There is no general exact closed form for LRU misses over every shape; replay the address trace as the simulator does.

For positive dimensions and independently line-aligned matrix regions, let $$ D_A=\left\lceil\frac{rk}{\lambda}\right\rceil, \quad D_B=\left\lceil\frac{kc}{\lambda}\right\rceil, \quad D_C=\left\lceil\frac{rc}{\lambda}\right\rceil. $$ Every input/output line is touched, so cold line reads are at least $D_A+D_B+D_C$, and materializing all updated output requires at least $D_C$ line writes in this write-back model. The compulsory lower bound is therefore $D_A+D_B+2D_C$ transfers. If the cache can retain the entire touched working set, both traversals achieve it. Otherwise repeated loads or dirty evictions can increase traffic. The sample has eight lines per matrix, so its compulsory bound is $32$ transfers; neither $56$ nor $48$ reaches that full-residency bound.

Choosing a tile and stating assumptions. Three useful square tiles need about $3T^2$ words. A practical candidate uses $3T^2\le\alpha M$ with a safety factor $0\lt\alpha\lt1$, but actual resident lines, other live data, and set mapping matter. An unaligned row segment of length $T$ can require up to $\lceil T/\lambda\rceil+1$ lines; three tiles can therefore touch up to roughly $3T(\lceil T/\lambda\rceil+1)$ lines even before accounting for neighbors retained from earlier work. Packing and alignment can reduce this waste. Cornell's array-blocking workshop motivates keeping the input/output working blocks together in cache.

Under an idealized suitable layout/cache with enough capacity and headroom for those blocks, a full tile does $\Theta(T^3)$ arithmetic work for $O(T^2/\lambda)$ line movement when $T$ is at least on the scale of $\lambda$. For large square $n\times n$ inputs with compatible block geometry, this gives the classical blocked scale $$ O\left(\frac{n^3}{\lambda T}+\frac{n^2}{\lambda}\right). $$ Choosing $T=\Theta(\sqrt M)$ yields $O(n^3/(\lambda\sqrt M)+n^2/\lambda)$ under a tall-cache/layout regime where line fragments do not dominate. This is an idealized asymptotic analysis, not a guarantee for every $T,Q,\lambda$ accepted by the finite simulator or for a real set-associative processor cache. Poor tile choices can erase or reverse the benefit, as the capacity experiments show.

Space and simulator cost. The numerical kernel requires the matrices/output and constant scalar state; this version does not allocate packed tile copies. Cache payload capacity is $Q\lambda$ words with $O(Q)$ LRU metadata. Backing/returned arrays use $O(rk+kc+rc+\lambda)$ host words, plus resident lines. Validation, initial zeroing, and result copying are linear in stored array sizes. With expected constant-time hash lookup, the instrumented simulator additionally copies $\lambda$ words per line read/write, so its software runtime includes $O(U+\lambda(R+W_e+W_f))$ work beyond setup. That simulation overhead is not the uninstrumented kernel's hardware execution time.

Correctness by coverage. Half-open, clamped row/column/inner tiles partition each dimension without gaps or overlap. Every triple $(i,x,j)$ belongs to exactly one combination of outer blocks and one inner iteration. Thus every required product $A_{ix}B_{xj}$ is added exactly once to the correct row-major output address.

Partial-sum invariant. For a fixed output tile, before reduction block kk, each entry equals the sum over inner indices below kk. The kernel adds exactly the indices from kk to the clamped block end; it never resets the previous sum. After the last reduction block, the entry equals the full dot product. Within an output entry, inner indices are processed in increasing order in both traversals. Integer bounds make every partial sum exact in these programs.

Cache coherence invariant. For every address, the newest value is in its dirty resident line if one exists, otherwise in backing memory. Reads use resident data or copy the current backing line. Dirty eviction writes the newest values before removal. Final flush therefore leaves all C values materialized in backing memory for observation. Locality affects transfers, while these invariants preserve the answer.

Variants, Pitfalls, and Optimizations

Keep reduction state across tiles. Initialize output once, outside all reduction blocks. Reinitializing a C tile each time kk changes keeps only the last partial product. The diagram below contrasts the first output block's correct sum with that failure. Test with multiple reduction blocks; one-block examples cannot expose it.

Clamp all dimensions. A square example divisible by $T$ hides tail bugs. Use independent row, column, and inner end bounds, and compute row-major offsets with each array's actual leading dimension. $A$ has row stride $k$; $B$ and $C$ have row stride $c$. Confusing them can appear correct for square inputs and fail for rectangles.

Measure write traffic as well as misses. Fewer line reads need not mean fewer total transfers. A schedule that evicts dirty output more often can cancel its read savings. Flush resident dirty lines at the same boundary for both variants; omitting the tiled run's final flush would unfairly count a result still waiting in cache as already materialized.

Tile size is a tuning parameter. Larger tiles reduce block-loop overhead but increase the live footprint; very small tiles can waste lines or repeat scalar/panel loads. Capacity-fit estimates need headroom and line geometry. This fully associative model excludes set conflicts; real row strides and power-of-two dimensions can map many lines to the same set. Padding, alignment, associativity, and packing can change outcomes.

More levels need more than one block size. A large cache panel and a small register tile serve different purposes. Hierarchical blocking can target several cache levels, while a microkernel keeps part of C in registers and vectorizes contiguous work. Non-square tile dimensions can reflect matrix shape and panel reuse. The finite code uses one cubic tile size to keep the invariant and trace visible.

Loop interchange and tiling are separate choices. Changing i,j,k to i,k,j improves row-major stride before tiling is considered. The supplied comparison gives both variants the same inner order. Intel's memory-access optimization guide discusses both transformations and measurement. Do not attribute every improvement from a changed loop nest solely to the block boundaries.

Dependencies and aliasing. Independent dense output entries permit this reordering because each receives the same ordered reduction. In-place updates, triangular solves, stencils, and dynamic programs may have dependencies crossing tiles. Use legal dependency order, halos, or wavefronts as needed. GPU shared-memory tiling also requires correct ownership, synchronization, and boundary handling; a CPU cache does not have the same explicit programming contract.

Numerical behavior. The sample uses bounded integers so exact output comparison is appropriate. Floating-point reductions can change rounding when a transformation changes accumulation order, vector reduction structure, or fused operations. These loops preserve per-entry inner-index order, but a different register/parallel microkernel may not. Compare with a suitable tolerance and numerical error model when moving beyond the integer laboratory.

Avoid reading too much into the simulator. The modeled line width, fully associative LRU policy, write allocation, and zero-initialized output are explicit choices. Real prefetchers, multiple cache levels, TLB misses, coherence, non-temporal stores, store buffers, compiler transformations, and memory bandwidth are absent. Values held in CPU temporaries and metadata are outside the cache-data model. Production multiplication with C initialized in registers or specialized stores may have a different C read/write schedule.

Validate outputs before timing. Use an independent dot-product oracle, rectangular and tail dimensions, negative/zero values, and configurations that force eviction. Then benchmark uninstrumented kernels or established numerical libraries with repeated controlled measurements. The simulator teaches data movement; its map/list calls and line copies should not be mistaken for realistic GEMM timings.

Practice Problems and Real-World Use

Compare both traversals with an independent reference that computes each output as a direct dot product. Test non-square shapes, every dimension smaller than the tile, dimensions just above multiples, zero dimensions, all-zero/all-one data, negative values, and tiles larger than the whole problem. A tile covering every dimension should reproduce the untiled address order and all counters.

For cache verification, generate the reference address trace separately and replay recency using timestamps rather than the implementation's LRU container. Check all reads, dirty writes, hits, evictions, peak residents, and the final matrix. Exercise generic store misses and dirty evictions independently too: multiplication's immediate C load/store pair does not by itself cover every cache-policy branch.

All three implementations agree on $3851$ shared inputs: $3834$ valid cases and $17$ rejected cases. The tests compare every output entry with independent dot products and every cache counter with a timestamp-based replay. They cover exhaustive small shape/parameter combinations, randomized rectangular products, all three dimension limits including $64\times64\times64$, ragged tiles, tiny/large caches, partial lines, invalid limits/shapes, and bounded extreme values. Python also passed $36$ exact-type rejection checks, $3200$ generic cache load/store checks with independent backing-memory verification and idempotent flush, and $200$ value-independence comparisons. C++ ran with address and undefined-behavior sanitizers; all sample output matches byte for byte.

To explore a machine, first identify the actual scalar representation and memory layout, then measure candidate rectangular tiles and cache levels. Record output correctness, elapsed time, and relevant hardware counters. Keep source initialization, kernel, flush/observation boundaries, thread placement, and compiler choices consistent so an apparent gain has a clear cause.