# Group 34 HW1 Writeup
### Members (equal contribution to code, writeup)
- Po-Yu Chou ([e-mail])
<!-- - 2-level blocking, fine loop unrolling, experiment result figures -->
- Xin Lyu ([e-mail])
<!-- - Repack & Realigning, SIMD micro-kernel, Loop reordering, report writeup -->
- Willis Wang ([e-mail])
<!-- - Alternative matmul, log parser for figures, report writeup -->

# Results

Description: Ver 7: register optimization, block (64, 192, 128), fine unrolling

Size: 31    Mflops/s: 5227.94   Percentage: 11.67
Size: 32    Mflops/s: 5298.88   Percentage: 11.83
Size: 96    Mflops/s: 10446.90  Percentage: 23.32
Size: 97    Mflops/s: 8447.04   Percentage: 18.86
Size: 127   Mflops/s: 11389.30  Percentage: 25.42
Size: 128   Mflops/s: 11667.11  Percentage: 26.04
Size: 129   Mflops/s: 9444.32   Percentage: 21.08
Size: 191   Mflops/s: 10639.69  Percentage: 23.75
Size: 192   Mflops/s: 10855.86  Percentage: 24.23
Size: 229   Mflops/s: 10096.96  Percentage: 22.54
Size: 255   Mflops/s: 11442.96  Percentage: 25.54
Size: 256   Mflops/s: 11454.71  Percentage: 25.57
Size: 257   Mflops/s: 10193.27  Percentage: 22.75
Size: 319   Mflops/s: 11051.53  Percentage: 24.67
Size: 320   Mflops/s: 11094.19  Percentage: 24.76
Size: 321   Mflops/s: 10195.49  Percentage: 22.76
Size: 417   Mflops/s: 10900.90  Percentage: 24.33
Size: 479   Mflops/s: 11457.35  Percentage: 25.57
Size: 480   Mflops/s: 11456.22  Percentage: 25.57
Size: 511   Mflops/s: 11599.75  Percentage: 25.89
Size: 512   Mflops/s: 11497.50  Percentage: 25.66
Size: 639   Mflops/s: 12083.89  Percentage: 26.97
Size: 640   Mflops/s: 12107.62  Percentage: 27.03
Size: 767   Mflops/s: 11971.71  Percentage: 26.72
Size: 768   Mflops/s: 11999.32  Percentage: 26.78
Size: 769   Mflops/s: 11490.45  Percentage: 25.65
Average percentage of Peak = 23.65

See `Final Results` below for analysis.
# Core Optimization Strategies

## Blocking

### Overview
Problematic cache condition occurs when cache lines for future reusable cache elements are swapped out. This causes a significant slowdown as L2 or L1 cache lines may need to be re-populated with the line of contention.

To minimize this issue, we first implemented blocking, and experimented with different blocking dimensions, in order to target more workload to the L2 cache. The size of L2 cache is $1MB$, which, in the ideal case, can store at most $\frac{10^{6}}{8}$ doubles. Say that we want to block matrices $A,B,C$ into blocks of size $b\times b$. For multiplying two blocks in $A,B$ and adding the results to a block in $C$, we want to load three $b\times b$ matrics into the L2 cache. Therefore, we can afford to choose $b$ such that $3b^2\le \frac{10^6}{8}$, which roughly upper bounds $b$ by $b\le 200$. Otherwise, we want to choose $b$ as large as possible so that the data movement between memory and cache is minimized. In our implementation, we set $b$ as $128$, a power of $2$, so that it can save some border checks in the microkernel computation. The experiment also showed that setting $b$ to larger values (e.g., $b = 160, 180, 200$) got fluctuated performance, and did not improve the performance in general.

Below, we document the performance of blocking, without other optimizations, when trying various blocking sizes and dimensions.

### Square blocking
We first experimented with multiple square blocking dimensions as well as transposing the `A` matrix with the starter code's `dgemm-blocked` implementation. We compared performance metrics between the different square dimensions in the below figure.]


The chosen ```BLOCK_SIZE=32``` targets the $32KB$ L1 cache (as the L1 cache can hold 4096 doubles), and ```BLOCK_SIZE=128``` targets the $1MB$ L2 cache (as the L2 cache can hold ~130k doubles). We also compare the effect of transposing the blocks in A so that the row elements are in contiguous memory.

Using solely blocking, we see an improvement in the transposed versions - as this was using the naive version, transposing the matrix will lead to contiguous memory accesses when iterating through the A matrix, making it much more cache-friendly.

Furthermore, targeting the L1 cache seemed to be more successful with only blocking. We think this may be due to inherent cache-benefits outside

### Rectangular blocking

Given our previous success in targeting the L1 cache, we investigate various rectangular blocking dimensions sized for our $32KB$ L1 cache.

The below table shows the total number of doubles in blocked A, B, and C matrices, in comparison to the theoretical maximum number of doubles able to fit into the L1 cache, 4096 doubles.

Block C: M-by-N, Block A: M-by-K, Block B: K-by-N

| (M, K, N) | (40, 24, 40) | (32, 32, 32) | (40, 32, 40) | (40, 40, 24)
| -------- | --------- | -------- | -------- | -------- |
| MK + KN + MN | 3520 doubles | 3072 doubles | 4160 doubles | 3520 doubles |

(40, 24, 40) vs (32, 32, 32): rectangle provides larger flexibility to reach the cache capacity.
(40, 24, 40) vs (40, 32, 40): the performance does not improve if MK + KN + MN exceeds capacity.
(40, 24, 40) vs (40, 40, 24): even if MK + KN + MN is the same, the dimension of block A and B matters.


<!-- [TODO: Summarize. Talk about how/why dimension actually matters (how blue is much better than red even with the same area)]. Marginal differences, cite paper. -->

<!-- As compared to the previous figure, the L1-targeted square blocking seems to outperform the rectangular blocking schemes we have tried. Even so, these are all marginal differences, with other optimizations impacting the optimality of the performance -->

### 2-Level blocking
Theoretically, by adding two-level blocking, we are able to target the L2 cache with the outer blocking, and L1 cache with inner blocking. Therefore, we wrote a 2-level blocking implementation, where we compare performance figures to 1-level L1-targeted and L2-targeted blocking below.


Unexpectedly, the performance of 2-level blocking lies that of between L1-targeting blocking and L2-targeting blocking, with an improvement in performance when compared to L2-targeted blocking, while still performing worse than L1-targeted blocking. This can be attributed to the extra work and logic needed to implement 2-level blocking, such as more instruction prefetching mispredictions with the increased amount of `for` loops alongside instruction cache contention, leading to no improvements in performance as compared to 1-level L1 blocking.

Due to the complexity involved in integrating this optimization to other strategies, alongside the insignificant changes in initial performance, we decide to stick with single-level blocking in the following sections.

### Summary
In summary, implementing **only** blocking yielded minimal increases in performance, with performance ranging from 0.5% to 1.5% of the theoretical peak. We believe that this may be due to non-aligned matrices in memory - this may cause matrix columns to use more cache lines than necessary, as well as not fully utilizing each cache line. Because of this, shrinking our block size was necessary, thereby causing more potential opportunities for instruction-level cache misses and increasing the number of instructions run. As such, we then implement realigning and repacking to solve this issue.

## Matrix Repacking + Realigning

### Overview
For this assignment, there is no guarantee for the given `A, B, C` matrices to be aligned to the 64-byte cache boundaries. For example, when loading and using a 64-byte matrix column, while it is theoretically able to be fit in one line of the cache, it might not be the case in real execution. As a consequence, cache lines will not be fully utilized when loaded in, and more cache contention will be caused, as more cache lines will be loaded in and swapped than ideally necessary when handling continuous matrix accesses. This motivates us to design and implement another optimization: realigning matrices in memory.

When we perform matrix realigning, we will also store the matrix in a special format, illustrated in the figure below. We call this step "matrix repacking".


The benefit of doing matrix repacking is that it allows us to access contiguous memory when calculating block multiplication. Note that inside a block, we store both matrix $A$ and $B$ in the column-major, and block-wise $A$ is arranged in the row-major, while $B$ is arranged in the column-major. The logic behind our design depends on the way we carry out the matrix multiplication. For the arrangement inside a block, we will elaborate on our design in the next section. For now, we explain the arrangement across blocks. Suppose $A,B,C$ are three $lda$-by-$lda$ matrices, where $lda = N \times b$. Recall in the blocking step, we partition each of $A,B,C$ into $N\times N$ blocks of size $b\times b$. The pseudocode below shows our algorithm for doing dgemm based on the blocking.

for i in [1, N] :
    for j in [1, N] :
        for k in [1, N] :
            Multiply (block) A(j,k) with (block) B(k,j);
            Add the result to (block) C(i, j);

Now, observe that fixing $i, j$ and letting $k$ run from $1$ to $N$, we are multiplying a row of blocks in $A$ with a column of blocks in $B$. Therefore, for every fixed $i, j$, by our repacking method, it is guaranteed that we are always accessing a contiguous part of $A$ and $B$ in the "$k$"-loop. This design adds some spatial locality to our algorithm, which we hope will further improve the overall performance of our algorithm.

**Work overhead:** Even though it may introduce an overhead of `O(n^2)` in realigning and repacking, this still ultimately leads to performance improvements, attributed to the fact that the workload is outshadowed by the matrix multiplications' complexity of `O(n^3)`

### Results
The below figure compares the performance between repacked+realigned and the non-repacked/realigned algorithm, as well as different blocking sizes. A naive micro-kernel (see Micro-kernel (naive)) and loop reordering were also used when running the benchmark.

Before substituting the naive micro-kernel with SIMD one, we need to have our data aligned in the memory. Meanwhile, we can also repack the data so that they follow the order of the computation of blocks.

The experiments below uses a naive micro-kernel (detailed below). Since we will copy elements of C to an aligned memory every time we are computing the block of C, this creates an overhead and hurts the performance consequently.


Unexpectedly, experiments without R&R performed better than those with R&R counterparts of the same block size. This leads us to believe that repacking and re-aligning is not be beneficial without further optimizations, such as an aligned SIMD-intrinsic micro-kernel. This can be perhaps be attributed to the overhead of realigning and repacking, alongside negligible cache benefits, given our chosen blocking dimensions.

Furthermore, in contrast to the results we obtained earlier (with only blocking implemented), there is a very noticeable benefit in L2-targeted blocking with these optimizations, as compared to more benefits of L1-targeted blocking with only blocking. This can be mainly attributed to how the addition of a microkernel, which, even with a naive implementation, can take advantage of L1 cache locality in performing memory accesses and operations.

## SIMD Microkernel

### Overview
Now let's consider multiplying two blocks. To further speed up the computation, we will exploit the instruction-level parallelism and write a SIMD microkernel.

The "atom operation" that our microkernel does is the following: it takes a value $B_{k, j}$, loads a vector of $8$ double $(A_{i,k},\dots, A_{i+7, k})$, multiplies each entry of the vector with $B_{k, j}$, and adds the result to $C_{i,k},\dots, C_{i+7, k}$. These operations can be carried out by standard AVX instructions, including `load`, `broadcast`, `fma`, and `store`. See the figure below for an illustration.


Recall that $C, A, B$ are all stored with column-major ordering inside each block. In our computation, we will enumerate $j, k$ on two outer loops, and then enumerate $i$ on the inner-most loop. Then, note that during the execution of the inner loop, we are always using the same entry of $B$. We can broadcast this entry into a vector and use that vector for the whole loop. For $A$ and $C$, we are accessing a column in each of them, which is stored in a contiguous section of memory, potentially making the data retrieving faster.

Our algorithm is formally shown in the following pseudocode.

# A, B and C are three b-by-b matrices
# Perform C = C + A * B
# Assume b is a multiple of 8
for j in [0, b-1]:
    for k in [0, b-1]:
        vecB = broadcast(B[k, j])
        for i in {0, 8, 16, ..., b - 8}:
            vecA = load(A[i:i+7, k])
            vecC = load(C[i:i+7, j])
            vecC = vecC + vecA * vecB
            store(vecC, C[i:i+7, j])

Furthermore, as we have aligned the A, B, and C matrices accesses to cache line boundaries as specified earlier, we are able to use the aligned AVX intrinsics to perform the operations detailed above. Notably, we use `_mm512_load_pd` to load a vector of C, `_mm512_fmadd_pd` to perform the fused multiply-add step, and `_mm512_store_pd` to store the resulting vector back to memory.

### Results
The previous plot showed how repacking & re-aligning did not yield benefits without a SIMD micro-kernel. However, the following plot shows how these two optimizations work together to have a large net improvement.


Overall, we see a massive improvement using SIMD intrinsic instructions (blue, green) versus without (orange, red), with a difference of ~8% theoretical peak performance with and without SIMD. This can be attributed to SIMD's inherent parallel computation of operations providing a variety of cache and CPU benefits, from contiguous memory accesses to taking advantage of the CPU's vector registers and vector ALU, overall reducing the number of CPU cycles used for memory load and arithmetic operations.

Another thing to note is that similar to our previous findings, we see a greater benefit in targeting the L2 cache for blocking versus the L1 cache.

## Optimizing the block multiplication

We further optimize the multiplication of two blocks. Recall the pseudocode so far.

def do_block(A, B, C):
    # Assume b is multiple of 8
    for j in [0, b-1]:
        for k in [0, b-1]:
            vecB = broadcast(B[i, k])
            for i in {0, 8, 16, ..., b - 8}:
                vecA = load(A[i:i+7, k])
                vecC = load(C[i:i+7, j])
                vecC = vecC + vecA * vecB
                store(vecC, C[i:i+7, j])

### Loop Unrolling

One obvious optimization is to do loop-unrolling to further utilize ILP (in particular, exploit two vector processing lanes of the architecture).

We will do an 8-way loop unrolling in our implementation. See below for an example of 2-way unrolling.

def do_block(A, B, C):
    # Assume b is multiple of 16
    for j in [0, b-1]:
        for k in [0, b-1]:
            vecB = broadcast(B[k, j])
            for i in {0, 16, ..., b - 16}:
                # load
                vecA0 = load(A[i:i+7, k])
                vecA1 = load(A[i+8:i+15, k])
                vecC0 = load(C[i:i+7, j])
                vecC1 = load(C[i+8:i+15, j])
                # fma
                vecC0 = vecC0 + vecA0 * vecB
                vecC1 = vecC1 + vecA1 * vecB
                # store
                store(vecC0, C0[i:i+7, j])
                store(vecC1, C1[i+8:i+15, j])

By conducting an 8-way unrolling, we can achieve an average performance of 15\% of the theoretical peak as shown in the Final Results graph below.

### Saving the write to memory

The bottleneck in the computation above might be the time required to write back to matrix $C$. Observe that right after every FMA operation (vecC += vecA * vecB), we will write vecC back to the matrix. Since there are $\frac{n^3}{8}$ FMA ops, it implies the writeback is performed $\frac{n^3}{8}$ times. This is really expensive, particularly in the case that $C$ is not loaded in $L1$ cache.

There is a relatively easy fix to this issue. Since a vector can handle $8$ doubles and we did an $8$-way unrolling, there is really no $i$-loop when we set $b$ as $8\times 8 = 64$. Then, after we enumerate the index "j", we can explicitly load a column of $C$ (say, $C[0:63, j]$) into registers. Then we loop over $k$ and do a bunch of fma operations on these registers. After the $k$-loop is finished, we write the registers back to the matrix. See the pseudocode below for an illustration.

def do_block(A, B, C):
    # Assume b = 64
    for j in [0, b-1]:
        vecC[0:7] = load(C[0:63, j])
        for k in [0, b-1]:
            vecA[0:7] = load(A[0:63, k])
            vecB[0:7] = broadcast(B[k, j])
            vecC = vecC + vecA * vecB
        store(C[0:63, j], vecC[0:7])

This optimization can boost the performance to around 22\% of the theoretical peak as shown in the Final Results graph below.

# Final Results

Note: LMK if you want to check out our versions on github - the repo is currently private, but we are glad to add you in.

- Green (```Ver 4.2```)
- This version loops in j-k-i order. It also uses repacking and realigning so that it is compatible with the intrinsic functions. Inside that micro-kernel, 8 elements in the column of C are calculated at the same time. In the loop that evokes the micro-kernel, manual unrolling is used so that 8 micro-kernels are being executed in one loop iteration.
- Orange (```Ver 5```)
- This version builds upon the ```Ver 4.2``` by incorporating register optimization, which involves reordering the loops to j, i, k. Besides, using different block dimensions, (M, K, N) = (64, 192, 128), also gives it a slight performance increase.

- Blue (```Ver 7```)
- This version builds upon the ```Ver 5```. While ```Ver 5``` unrolls loops if M=64, where M is the length of the column in a particular block in C, so that 8 AVX-512 instructions can be executed at once, otherwise micro-kernels are invoked one by one. In ```Ver 7```, we further unroll the loop for the case 32<=M<64, which guarantees 4 AVX-512 instructions; for 16<=M<32, it can be unrolled by 2, but this adds little performance improvement. We refer to this technique as "```fine loop unrolling```" in contrast to that used by ```Ver 5```, which is referred to as "```coarse loop unrolling```". By unrolling the loop with smaller M, we can boost the performance when dealing with smaller matrices.


# Miscellaneous Optimizations

## Loop reordering

As a preliminary treatment for incorporating the micro-kernel, we need to reorder the loops. To our surprise, we make an unexpected performance improvement for the naive GEMM, which further extends to our more advanced optimizations detailed above.


As for why simply changing the order of loops results in such tremendous performance gain, we are not certain about that. One explanation may be the compiler realizes that block A is in contiguous memory and is able to perform loop unrolling automatically.

## Micro-kernel (naive)
Prefacing our SIMD microkernel, we decided to implement a naive micro-kernel that performs 8 calculations of C entries at once to leverage the reordered loop (j, k, i). Furthermore, it serves as a placeholder for our SIMD micro-kernel.


A very interesting thing to note about the above graph is how a simple reordering of loops with the starter code, shown with the blue line, can edge out performance improvements gained from more complex optimizations, such as L1-targeted blocking and naive microkernel.

Even so, we see that adding a naive microkernel still yields positive benefits as compared to non-microkernel implementations, due to a naive microkernel by itself as a form of loop unrolling, reducing the number of loop iterations and making the code more prefetcher-friendly.

## Loop Unrolling

To further boost the performance, we unrolled the loop that calls the SIMD micro-kernel by 8.


We find it extremely surprising how simply performing loop reordering can beat the performance gains from L1-targeted blocking the microkernel. This experiment truly brings a humbling moment, where after all, computers are very complicated and often behave like an enigma.

## References
Goto, K., and van de Geijn, R. A. 2008. Anatomy of High-Performance Matrix Multiplication, ACM Transactions on Mathematical Software 34, 3, Article 12.