Pair: Alexander Cueva & Bernard Li
This report details implementing, optimizing, and analyzing the performance of matrix-matrix multiplication. We explored various optimization techniques, including algorithmic changes and compiler flags, and to evaluate their effectiveness using the Roofline performance model.
While the assignment prompt mentioned matrix-vector multiplication, we decided to go the extra mile out of interest to focus on the general case of matrix-matrix multiplication (C = A * B).
This repository has:
- Three C++ implementations for matrix multiplication: a naive baseline, a loop-ordered optimization, and a block multiplication optimization.
- A
Makefileto compile the code with different optimization flags (-O1,-O2,-O3). - A Python script (
roofline.py) to generate Roofline models from the performance data. - CSV files containing the results of performance measurements for various matrix sizes and optimization levels.
- PNG images of the generated Roofline plots.
Three different C++ programs were used to measure the performance of matrix-matrix multiplication. All implementations use a single, contiguous block of memory for each matrix (e.g., new long double[rows * cols]), which is an efficient strategy that promotes good spatial locality. We found that earlier attempts with non-contiguous memory would make our program too slow and inefficient in
creating measurements, especially with rather large dimensioned matrices.
This file contains the baseline, "naive" implementation using a standard ijk triple-loop structure. This approach is simple to write but is known to have poor cache performance because it accesses the second matrix (B) in a column-wise fashion, leading to non-sequential memory access and frequent cache misses.
// Naive ijk loop order
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
for (int k = 0; k < shared; k++) {
result[i * cols + j] += matrixA[i * shared + k] * matrixB[k * cols + j];
}
}
}
This version reorders the loops to an ikj structure. This is an optimization for row-major matrices. In this ordering, the innermost loop iterates over j, resulting in sequential access patterns for both the result matrix C and the second input matrix B. The element from matrix A is fetched once and reused across the entire inner loop, significantly improving cache performance.
// Optimized ikj loop order
for (int i = 0; i < rows; i++) {
for (int k = 0; k < shared; k++) {
long double a_ik = matrixA[i * shared + k];
for (int j = 0; j < cols; j++) {
result[i * cols + j] += a_ik * matrixB[k * cols + j];
}
}
}
This optimization uses a block matrix multiplication (aka tiling). The matrices are divided into smaller, cache-sized blocks. The computation proceeds block by block, which dramatically increases data reuse due to temporal locality. Within each block, the ikj loop ordering is used to maintain efficient sequential memory access due to spatial locality.
// Outer loops iterate over blocks
for (int ii = 0; ii < rows; ii += blockSize) {
for (int kk = 0; kk < shared; kk += blockSize) {
for (int jj = 0; jj < cols; jj += blockSize) {
// Inner loops perform ikj multiplication on the block
for (int i = ii; i < iMax; i++) {
for (int k = kk; k < kMax; k++) {
// ...
}
}
}
}
}
The results stored in the results_*.csv files show that compiler optimizations have a massive impact on performance.
- No Flags: Performance is uniformly poor across all implementations.
- -O1: Provides a significant boost over the unoptimized version.
- -O2: Offers further improvements. This level enables a large number of optimizations, including instruction scheduling and automatic vectorization, that are highly effective for numerical code.
- -O3: Provides a marginal improvement over
-O2and in some cases can be slightly slower. For this reason,-O2is often the recommended "go-to" optimization level.
The block multiplication and ikj versions consistently outperform the naive ijk version at every optimization level, implying that while compiler optimizations are powerful, a better algorithm ultimately provides a stronger basis.
The Roofline model helps visualize the performance of an application in the context of the hardware's limitations—specifically, its peak memory bandwidth and peak computational rate. For this assignment, the assumed hardware specifications were:
- Peak Memory Bandwidth: 160 GB/s
- Peak FLOP Rate: 1 TFLOP/s (1000 GFLOP/s)
The ridge point of the model occurs at an arithmetic intensity (AI) of 1000 / 160 = 6.25 FLOPs/byte.
- If the AI of an application is less than 6.25, its performance is limited by memory bandwidth (it is memory-bound).
- If the AI is greater than 6.25, its performance is limited by the CPU's peak computational rate (it is compute-bound).
Figure 1: Roofline model for the block multiplication implementation with -O3 optimization.
The plot shows that for all tested matrix sizes, the AI is greater than the ridge point. This means the matrix-matrix multiplication is compute-bound. Our goal is therefore to move the performance points vertically towards the "compute roof" of 1000 GFLOP/s.
The results show that even our best-optimized code does not reach the theoretical peak. However, the plots clearly show that the optimized versions (ikj and block) achieve higher performance (are higher on the chart) than the naive version. We should be able to reach the theoretical peak (or become memory-bound) with parallelization, which on CPUs can be via OpenMP.
Overall, we see that a simple change in loop ordering (ijk to ikj) provided a substantial performance boost by aligning memory access patterns with the row-major storage format. However, tiling provided the best performance by maximizing both spatial and temporal cache locality.
Compiler flags like -O2 and -O3 provide huge performance gains, but they cannot fix a fundamentally inefficient algorithm. The best results are achieved when an efficient algorithm is paired with aggressive compiler optimization. One more compiler optimization we did not implement in our code is the use of loop unrolling (#pragma unroll). This is especially effective on smaller matrices which can be easily unrolled and bounds entirely fit into local cache memory.
The Roofline model showed that our application is compute-bound, meaning further significant improvements would require more advanced techniques like parallelization via OpenMP techniques.
For all of these, we can see that without the optimizations, the performance is quite poor. For better performance, we need a more efficient, more parallelizable algorithm that will allow us to cut down on the computation cost (or at the very least parallelize it).
