Skip to content

Complete Mathematics of "Optimization 2: Unrolling by 4 (Column Blocking)"

This optimization introduces manual loop unrolling on the column dimension of matrix \(C\). The goal is to improve Instruction Level Parallelism (ILP) and reduce loop overhead, though the fundamental memory access patterns remain similar to Optimization 1.

1. Mathematical Formulation

The Operation

The code computes four adjacent elements of \(C\) simultaneously for a fixed row \(i\):

\[ \begin{aligned} C_{i, j} &\mathrel{+}= \sum_{p=0}^{k-1} A_{i,p} B_{p, j} \\ C_{i, j+1} &\mathrel{+}= \sum_{p=0}^{k-1} A_{i,p} B_{p, j+1} \\ C_{i, j+2} &\mathrel{+}= \sum_{p=0}^{k-1} A_{i,p} B_{p, j+2} \\ C_{i, j+3} &\mathrel{+}= \sum_{p=0}^{k-1} A_{i,p} B_{p, j+3} \end{aligned} \]

Vectorized Representation

Instead of computing a scalar dot product, we can view the inner loop as updating a vector of 4 elements \(\mathbf{c}_i = [C_{i,j}, C_{i,j+1}, C_{i,j+2}, C_{i,j+3}]^T\):

\[ \mathbf{c}_i \mathrel{+}= \sum_{p=0}^{k-1} A_{i,p} \begin{bmatrix} B_{p,j} \\ B_{p,j+1} \\ B_{p,j+2} \\ B_{p,j+3} \end{bmatrix} \]

This is mathematically equivalent to the previous version but processed in chunks of 4 columns.

2. Loop Structure & Control Flow

Original Loop (Optimization 1)

\[ \text{Iterations} = m \times n \]

Each iteration computes 1 dot product.

Unrolled Loop (Optimization 2)

\[ \text{Outer Iterations} = m \times \frac{n}{4} \]

Each iteration computes 4 dot products.

Total Dot Products:

\[ m \times \frac{n}{4} \times 4 = m \times n \]

The arithmetic workload remains exactly \(2mnk\) FLOPs.

3. Memory Access Analysis

This is the critical part for understanding the performance gain.

A. Access to Matrix A (The "Broadcast" Effect)

In the inner loop (inside AddDot), we access \(A_{i,p}\).

  • Address: \(p \cdot \text{lda} + i\)
  • Stride: \(\text{lda} = m\) (Column-major layout)

Since we call AddDot 4 times with the same &A(i,0) and lda, the CPU loads the exact same sequence of \(A\) elements 4 times.

Cache Behavior:

  1. First Call: Loads \(A_{i,0}, A_{i,1}, \dots\) from L1/L2 cache (or RAM if cold).
  2. Subsequent 3 Calls: Since \(k\) is large (2000), \(A\)'s column might not fit in L1 cache (\(2000 \times 8 \text{ bytes} = 16 \text{ KB}\)). However, the data is likely in L2.
  3. Optimization: The CPU doesn't need to fetch A from RAM 4 times, but it does need to load it into registers 4 times.

B. Access to Matrix B (Contiguous Reads)

We access \(B_{p,j}, B_{p,j+1}, B_{p,j+2}, B_{p,j+3}\).

  • Address: \(j \cdot \text{ldb} + p\)
  • Stride: 1 (Contiguous in memory for fixed \(p\), varying \(j\))

Wait! In Column Major layout (B(i,j) = b[j*ldb + i]), varying \(j\) (the outer loop here) means jumping by ldb. Varying \(p\) (the inner loop inside AddDot) means stride 1.

So, inside AddDot for \(B_{p,j}\), we read contiguous memory. For \(B_{p,j+1}\), we read a different contiguous block.

The Gain: The CPU can potentially prefetch or keep the stream of \(B\) data flowing smoothly because the access pattern inside AddDot is perfectly contiguous for each of the 4 columns.

C. Access to Matrix C (The Accumulators)

We update \(C_{i,j}, \dots, C_{i,j+3}\).

  • These are 4 distinct memory locations.
  • In the unrolled version, the compiler can keep these 4 values in registers (specifically, floating-point registers) for the duration of the inner loop \(p\).

4. Why is this "Optimization 2"?

The primary benefit is Instruction Level Parallelism (ILP).

Dependency Chain Problem (Optimization 1)

In a simple dot product:

\[ \gamma = \gamma + A_p B_p \]

The addition depends on the previous value of \(\gamma\). This creates a dependency chain. The CPU must wait for the addition to finish before starting the next one.

  • Latency of FP Addition: ~4 cycles.
  • Throughput: 1 add per cycle (ideally).
  • With dependency, effective throughput drops to 1 add per 4 cycles.

Breaking Dependencies (Optimization 2)

By unrolling by 4, we have 4 independent accumulators:

\[ \gamma_0 \mathrel{+}= A_p B_{p,j} \]
\[ \gamma_1 \mathrel{+}= A_p B_{p,j+1} \]
\[ \gamma_2 \mathrel{+}= A_p B_{p,j+2} \]
\[ \gamma_3 \mathrel{+}= A_p B_{p,j+3} \]

These 4 additions are independent. The CPU can execute them in parallel (superscalar execution).

  • Cycle 1: Issue Add 1, Add 2, Add 3, Add 4
  • Cycle 2: Issue Mult 1, Mult 2, Mult 3, Mult 4
  • ...

This allows the CPU to saturate its Floating Point Units (FPUs), significantly increasing GFLOPS.

5. Arithmetic Intensity & Roofline Model

Metric Optimization 1 Optimization 2
FLOPs \(2mnk\) \(2mnk\)
Memory Traffic (A) High (Strided) High (Strided, but reused 4x per row iteration)
Memory Traffic (B) Low (Contiguous) Low (Contiguous)
Loop Overhead High (\(mn\) iterations) Low (\(mn/4\) iterations)
Register Pressure Low Medium (4 accumulators)
ILP Low (Dependency chain) High (4 independent chains)

The "Unrolling Factor" Trade-off

Why 4?

  1. Registers: x86-64 has 16 SSE/AVX registers. Using 4 for accumulators leaves plenty for pointers and temps.
  2. Code Size: Unrolling by 8 or 16 increases code size (I-Cache pressure).
  3. Diminishing Returns: After 4-8, the dependency chain is sufficiently broken, and memory bandwidth becomes the bottleneck again.

6. Summary of Mathematical Changes

  1. Loop Bound: \(j\) goes from \(0 \to n\) with step 4.
  2. Operation: Scalar Dot Product \(\to\) Vector-Scalar Update. $$ \begin{bmatrix} C_{i,j} \ C_{i,j+1} \ C_{i,j+2} \ C_{i,j+3} \end{bmatrix} \mathrel{+}= A_{i,p} \begin{bmatrix} B_{p,j} \ B_{p,j+1} \ B_{p,j+2} \ B_{p,j+3} \end{bmatrix} $$
  3. Performance Driver: Reduced loop branch overhead and increased Instruction Level Parallelism via independent accumulators.

Limitation

This optimization does not fix the strided access to Matrix \(A\). \(A\) is still accessed with stride \(m\), which is cache-unfriendly. A true "Optimization 3" would likely involve Loop Interchange or Blocking to fix the access to \(A\).