Blocked Matrix Multiply (DPAS)

• Blocked Matrix Multiply (DPAS) is an algorithm that partitions dense matrices into b×b blocks to enhance cache locality and improve GEMM efficiency.
• The method reorders computations via various traversal schemes (k-outer, row-major, Morton) to maximize temporal reuse and minimize memory bandwidth constraints.
• Empirical results show that selecting block sizes between 16 and 64 can yield 2–4× speedups over naive multiplication, demonstrating its effectiveness across platforms.

Blocked matrix multiply (also known as DPAS-style blocked matrix–matrix multiplication) is an algorithmic and architectural strategy for computing the product of two dense $N \times N$ matrices via partitioning them into $b \times b$ submatrices (blocks) and orchestrating both storage and traversal to optimize temporal and spatial locality. The method leverages the reuse of small matrix tiles in cache to enhance computational throughput and to minimize main-memory bandwidth limitations during general matrix–matrix multiplication (GEMM) routines. Blocked multiply strategies are central to high-performance dense linear algebra libraries, and cache-oblivious heuristics are vital for generalizing these gains across architectures without relying on low-level memory layout choices (0808.1108).

1. Blocked Matrix Multiply Formalism and Index Partitioning

Let $A$, $B$, and $C$ be dense $N \times N$ matrices to be multiplied as $C = AB$. Select a block size $b$ such that $n = N / b$ is integral; otherwise, matrices are padded to the next multiple of $b$. Define submatrices as

$b \times b$0

and analogously for $b \times b$1 and $b \times b$2 for $b \times b$3. The blocked multiply is then given by

$b \times b$4

where each submatrix operation reduces to a call to GEMM on $b \times b$5 blocks.

2. Algorithmic Structure and Traversal Orders

The essence of blocked multiply is captured by iterating over all index triples $b \times b$6 in a prescribed order:

• k–outer loop:

$n = N / b$4

• i–j–k (row-major output order):

$n = N / b$5

• Morton/Peano recursion: Visit $b \times b$7 in Z-curve or space-filling curve order, recursively descending on $b \times b$8 block quadrants in $b \times b$9-space, thereby maximizing temporal and spatial locality by minimizing working set size.

The particular traversal order (the "ORDER") of $A$0 indices—rather than memory layout of blocks—is the determining factor for performance improvement when block sizes are sufficiently large to exploit cache locality (0808.1108).

3. Storage Layouts for Blocked Submatrices

Four canonical strategies exist for arranging the $A$1 blocks of matrices $A$2, $A$3, and $A$4 in physical memory (element size $A$5 bytes):

Strategy	Block Storage	In-block Order
Contiguous row-major	Linear in $A$6	row-by-row
Contiguous column-major	Linear in $A$7	column-by-column
Interleaved row-stripe	Rows of blocks contiguous (with interleaving for associativity)	arbitrary
Morton/Peano (Z order)	Blocks placed by interleaved $A$8 bits	arbitrary

• Row-major: Block $A$9 at $B$0, element $B$1 within block at $B$2.
• Column-major: Same scheme with $B$3 reversed.
• Interleaved row-stripe: Blocks $B$4 are contiguous, with adjustable stride for cache associativity.
• Morton/Z-order: Block $B$5 at $B$6; Morton index $B$7 calculated by bit-interleaving digits of $B$8 and $B$9.

Within each block, use whichever layout (row-major or column-major) is optimal for the backend GEMM implementation. In practice, results reveal that the storage layout of blocks becomes irrelevant for performance as block size increases; execution order is the controlling parameter (0808.1108).

4. Execution Order Heuristics and Cache Model

Execution-order heuristics critically impact cache usage and data reuse:

• Perfect-locality test: Repeats $C$0, keeping all blocks resident in cache.
• No-locality test: Randomizes all $C$1 tuples and block allocations, maximally destroying temporal and spatial reuse.
• k-ordered: $C$2 outermost, $C$3 innermost; maximizes temporal locality for $C$4 and $C$5 across multiple $C$6 updates.
• i–j–k (row-major): Processes all $C$7 for every $C$8 output block before moving on.
• Morton/Peano recursion: Maintains a small working set by recursive subdivision of blocks.

The cache model is characterized as follows:

• Block size $C$9 with element size $N \times N$0 implies each block is $N \times N$1 bytes.
• If $N \times N$2 ($N \times N$3 = cache size), blocks persist in cache across inner loops, ensuring temporal reuse.
• Modern prefetchers provide sufficient spatial reuse within each block when $N \times N$4 is several cache lines.
• A block should not be evicted before completion of all necessary updates, otherwise temporal locality is lost.
• Cache misses per block-multiply are predicted as
  • $N \times N$5 for cold blocks,
  • $N \times N$6 if two blocks remain cached from prior iteration.

5. Quantitative Performance and Scaling

Performance measurements on AMD Opteron 248 (2.2 GHz, 1 MiB L2) and Xeon Woodcrest (2.66 GHz, 4 MiB L2) use

$N \times N$7

where $N \times N$8 is CPU time. Speedup relative to "no locality" baseline is

$N \times N$9

$C = AB$0	no-loc	temp-only	full Peano	perfect-loc
8	1.0×	2.8×	2.9×	3.0×
16	1.0×	2.5×	2.6×	2.7×
32	1.0×	1.9×	2.0×	2.1×
64	1.0×	1.3×	1.4×	1.4×
128	1.0×	1.1×	1.1×	1.1×

On both platforms, speedup saturates for $C = AB$1, dropping to $C = AB$2; for $C = AB$3, speedup reaches $C = AB$4–$C = AB$5 over the "no locality" order (0808.1108).

6. Implementation Guidelines and Practical Recommendations

Key findings and prescriptions for DPAS-style implementation:

• Temporal locality is dominant: Arrange $C = AB$6 traversal so that all block updates for a given $C = AB$7, $C = AB$8, $C = AB$9 triplet are completed before potential cache eviction.
• Spatial locality (contiguous memory for blocks) is only consequential for prefetch efficiency within each GEMM call; block-to-block alignment is immaterial for $b$0 with $b$1 B lines.
• Morton/Peano traversal realizes nearly ideal locality, but any schedule that processes all $b$2 for each $b$3 suffices.
• Block size tuning:
  • Require $b$4.
  • Avoid $b$5 so small that GEMM vectorization degrades ($b$6 typically $b$7–$b$8).
• Blocks may be stored in any array structure (row-major is simplest).
• Loop schedule: $n = N / b$6 or use a small-recursion Z-curve that tiles the $b$9 cube in $n = N / b$0 block-octants.

By adhering to these guidelines, one achieves the $n = N / b$1–$n = N / b$2 speedups for small $n = N / b$3 as empirically demonstrated, wholly independent of whether blocks are stored in specialized Morton or other layouts (0808.1108).

7. Conclusions and Implications

Blocked matrix multiply strategies prioritize temporal locality through thoughtful traversal of block operation triples, rather than through intricate block memory layouts. The major implication is that algorithmic scheduling affords near-optimal performance across platforms without the complexity of cache-oblivious data structures or elaborate physical reorganization. The performance benefits plateau with growing block size as spatial locality between blocks becomes negligible, and as BLAS-generic GEMM tuning dominates achievable throughput. Thus, future matrix–matrix multiplication designs should focus on execution order heuristics and cache residency constraints rather than specialized block placement, consistent with the empirical and analytic findings presented in the cited study (0808.1108).