This work presents a comprehensive performance analysis and progressive optimization of a Generic Matrix Multiplication Algorithm (GEMM) for large matrices running natively on a 64-bit RISC-V architecture utilizing the RVV 1.0 vector extension and advanced cache locality techniques (Tiling), ILP, Out Of Order execution and multithreading via OpenMP. All improvements are cumulative:
This work develops the algorithm through various stages:
GEMM_O2.c: Pure scalar code.GEMM_O3.c: Code vectorized by the compiler using the-O3flag.GEMM_vaddmul.c: First vector approach usingvle32.v,vfmul.vf,vfadd.vvwith LMUL of 8.GEMM_vmacc.c: Latency reduction by replacing the Mul/Add pair with the combined instructionvfmacc.vv. This version uses a fixed LMUL of 8.GEMM_tiled.c: Implementation of Loop Nest Blocking (Tiling) by dividing the spatial problem into submatrices that fit into the L1 cache, keeping the data in L1 and eliminating stores to RAM. This version also uses LMUL of 8.GEMM_unroll.c: Enabling ILP (Instruction Level Parallelism) with the implementation of loop unrolling with depth 4 and interleaving to enable OOO (Out-of-Order) execution, utilizing the VPU pipeline and mitigating RAW risks. To prevent register spilling, the LMUL is reduced to 4.GEMM_omp.c: Multithreading application using OpenMP, to make the most of the SBC’s 8 physical cores.
The experiment was carried out on a Banana BPI-F3 SBC with X60 SoC system SpaceMit K1 octacore, measuring the pure computation time required to calculate the resulting matrix
The benchmark matrices have been set to dimensions of
perf stat -e cycles,instructions,branches,branch-misses,L1-dcache-loads,L1-dcache-load-misses,L1-dcache-stores,L1-dcache-store-misses ./P3GEMM_version 2000 3000 100| Metric | GEMM_O2 (Base) |
GEMM_O3 |
GEMM_vaddmul |
GEMM_vmacc |
GEMM_tiled |
GEMM_unroll |
GEMM_omp |
|---|---|---|---|---|---|---|---|
| Computing Time | 2.881 s | 0.768 s | 0.781 s | 0.710 s | 0.320 s | 0.192 s | 0.043 s |
| Speedup | 1.0x | 3.75x | 3.69x | 4.06x | 8.99x | 15.02x | 66.58x |
| LMUL | - | 1 | 8 | 8 | 8 | 4 | 4 |
CPU Cycles (cycles) |
5,484 M | 2,177 M | 2,138 M | 2,010 M | 1,389 M | 1,234 M | 1,436 M |
| Instructions | 5,412 M | 1,562 M | 793 M | 782 M | 675 M | 674 M | 741 M |
IPC (insn per cycle) |
0.99 | 0.72 | 0.37 | 0.39 | 0.49 | 0.55 | 0.52 |
L1 Loads (loads) |
1,949 M | 466 M | 574 M | 574 M | 361 M | 223 M | 231 M |
| L1 Load Misses | 1.9 M | 1.8 M | 2.0 M | 2.0 M | 2.0 M | 2.4 M | 2.3 M |
L1 Stores (stores) |
689 M | 244 M | 280 M | 280 M | 91 M | 91 M | 91 M |
| L1 Store Misses | 568 K | 553 K | 550 K | 550 K | 546 K | 543 K | 549 K |
Key takeaway: Maximum optimization is not achieved simply by inserting vector instructions. Leaving aside the unfair comparison with the multithreaded version, the
GEMM_unrollversion, which uses only one core, achieves performance$4x$ times that of the automatically vectorized code, simply by applying reordering techniques and optimizing how data moves between main memory, the caches and the processor’s vector registers.
- Instruction reduction: Manual vector optimization (_vaddmul and _vmacc) and block-based optimization (tiled) significantly reduce the total number of instructions compared to automatic vectorization _O3. This is because the compiler is very conservative in its application, and uses an LMUL of 1 to avoid excessive pressure on registers. This can be examined in the assembly code generated by the compiler:
90 0082 D7F7060D vsetvli a5,a3,e32,m1,ta,ma # LMUL=1
91 .loc 1 20 8 is_stmt 1
92 .loc 1 20 22 is_stmt 0
93 0086 87600502 vle32.v v1,0(a0)
94 008a 93952700 slli a1,a5,2
95 .loc 1 20 40
96 008e 07610602 vle32.v v2,0(a2)
19:P3GEMM.c **** for (j = 0; j < m; j++)
97 .loc 1 19 24 discriminator 1
98 0092 9D8E sub a3,a3,a5
99 0094 2E95 add a0,a0,a1
100 0096 2E96 add a2,a2,a1
101 .loc 1 20 26
102 0098 D79021B2 vfmacc.vv v1,v3,v2 # Uses fused vmacc
103 .loc 1 20 16
104 009c A7600702 vse32.v v1,0(a4)- The misleading IPC of _O2: Although _O2 has the highest IPC (0.99), it is inefficient; the processor runs fast but only executes scalar instructions. Vector versions, on the other hand, have a much lower IPC, but execute much more complex instructions on a much larger amount of data.
-
Efficient operation fusion: The
GEMM_vmaccversion reduces the number of instructions compared to _vaddmul (from 793M to 782M), demonstrating that the multiply and accumulate operations are fused into a single CPU cycle. -
Numerical drift with separate addition and multiplication operations: in the case of
GEMM_vaddmul.c, which requires a significant number of multiplication and accumulation operations, a numerical drift of the order of$10^{-5}$ has been observed. Whilst this may be of little significance for certain types of applications, for scientific computing applications it can have a huge impact on how a numerical method converges. -
Loop Nest Blocking (Tiling): the tiled version achieves an
$\approx 9x$ speedup (down from 2.881 s to 0.320 s) compared to the base _O2 version. This performance improvement can be achieved by implementing an efficient load-and-store strategy. - Hardware-Saturating Instruction Interleaving: Implementing a row-wise loop unrolling factor of 4 effectively hides FMA (Fused Multiply-Add) latency, maximizing pipeline throughput and driving instruction per cycle (IPC) efficiency up to 0.55.
-
Store traffic:
GEMM_tiledreduces L1 cache writes from 689 million to 91 million (a reduction of 86.7%), confirming that partial sums are retained in the registers before being written to memory. - Data bus independence: _tiled reduces L1 reads (loads) to one-fifth of the base version, preventing the CPU from suffering from data starvation and raising its actual IPC to 0.49.
-
Consistency of store failures: Write failures (store-misses) remain constant (
$\approx$ 550 K) across all versions because they correspond to the initialization of the data structures.
A strong scaling experiment was conducted on the Banana Pi board to evaluate the parallel efficiency of the 2D Register-Blocked GEMM kernel. The matrix sizes for the benchmark were set to
The table below summarizes the execution time, speedup factor, and parallel efficiency achieved as the thread count scales from 1 to 8 cores, with OpenMP thread affinity fully enabled (OMP_PLACES=cores, OMP_PROC_BIND=close).
| Threads ( |
Compute Time (s) | Measured Speedup ( |
Ideal Speedup | Parallel Efficiency ( |
|---|---|---|---|---|
| 1 | 0.199959 | 1.00x | 1.00x | 100.0% |
| 2 | 0.109174 | 1.83x | 2.00x | 91.5% |
| 3 | 0.081090 | 2.47x | 3.00x | 82.3% |
| 4 | 0.065696 | 3.04x | 4.00x | 76.0% |
| 5 | 0.057543 | 3.47x | 5.00x | 69.4% |
| 6 | 0.051422 | 3.89x | 6.00x | 64.8% |
| 7 | 0.048268 | 4.14x | 7.00x | 59.1% |
| 8 | 0.043392 | 4.61x | 8.00x | 57.6% |
-
High Initial Efficiency (1 to 4 Cores): The application scales well up to 4 threads, maintaining a parallel efficiency of 76%. This confirms that the OpenMP work-sharing directive
#pragma omp parallel for collapse(2)successfully achieves a good load balance when distributing the matrix tiles across the available processors. -
Memory Bandwidth Bottleneck (5 to 8 Cores): Beyond 4 threads, the speedup curve experiences a standard sub-linear saturation, reaching a final acceleration of 4.61x at 8 threads. This behavior is a textbook example of a memory-bound limitation in Dense Linear Algebra. As more processing cores concurrently issue vectorized vector-register loads (
__riscv_vle32), the shared L2 cache system and the main RAM bus bandwidth become saturated, introducing slight memory stalls. -
Thread Affinity Impact: Forcing tight hardware constraints prevents the scheduler from migrating threads across different physical cores. This guarantees strict cache-locality for the
$64 \times 64$ tiles, eliminating cache thrashing and providing highly reproducible compute timings. -
Numerical Consistency: The output matrix
$C$ yields identical floating-point precision data across all thread configurations, verifying the mathematical correctness of the loop tail cleaning logic and ensuring that the implementation remains entirely race-condition free.
The Banana Pi features a 32 KB L1 data cache and operates with 32-bit floating-point precision (4 bytes per element). An empirical study was carried out by varying the tile size to find the hardware’s optimal saturation point. Here are the results for the tiled version without unrolling :
| Tile Size | Computation Time | L1 Loads | L1 Misses | % Misses | L1 Cache Status |
|---|---|---|---|---|---|
| 32 x 32 | 0.712 s | 588 M | 2.2 M | 0.39% | Underutilized |
| 50 x 50 | 0.310 s | 368 M | 2.7 M | 0.75% | Optimal net balance |
| 64 x 64 | 0.320 s | 361 M | 2.4 M | 0.68% | Alignment optimum |
| 128 x 128 | 0.452 s | 359 M | 8.3 M | 2.33% | Saturation and Overflow |
The algorithm’s behaviour is as expected, depending on the tile size:
- For small sizes that fit comfortably into L1, the number of load operations increases.
- For sizes close to the optimum, utilisation is maximised: the number of loads decreases whilst the cache miss rate remains relatively stable.
- For large tile sizes, which exceed the capacity of L1, the number of cache misses skyrockets.
To compile natively in the RISC-V environment using vector support (requires gcc with support for RVV 1.0):
# Compile all versions
make
# Compile one file
gcc -march=rv64gcv -Wall -O3 <GEMM_version.c> -o <GEMM_version>The objective of this work is to optimize the calculation of the sine of a 32-bit floating-point data array using the polynomial approximation of the Taylor series. The core of the algorithm requires iteratively calculating odd powers and accumulated factorials using the mathematical relationship:
-
Test Platform: Banana Pi BPI-F3 with Spacemit K1 processor, 64-bit RISC-V architecture with 256 bits Vector Extension 1.0.
-
Dataset Configuration: N=200,000 elements, 12 terms of the sine Taylor series.
By analyzing the compiler's behavior with the flags -fopt-info-vec-optimized and -fopt-info-vec-missed, we can see how the compiler behaves when attempting to vectorize the loop:
gcc -O3 -march=rv64gcv -fopt-info-vec-optimized -fopt-info-vec-missed P3sin.c -o P3sin
P3sin.c:26:14: missed: couldn't vectorize loop
P3sin.c:26:14: missed: not vectorized: unsupported control flow in loop.
P3sin.c:33:16: missed: couldn't vectorize loop
P3sin.c:35:26: missed: not vectorized: unsupported use in stmt.
P3sin.c:92:14: optimized: loop vectorized using variable length vectors
The only loop the compiler has been able to optimize is the data initialization loop, but it has not been able to vectorize any of the loops that perform the Taylor series calculation. Starting from this baseline, we will conduct a study of manual vectorization, varying the length of the vectorization to examine performance.
The use of vector intrinsics improves performance by about 10 times compared to code optimized with the -O3 flag. However, the value of LMUL does not affect performance; this is mainly due to two factors:
- The RAW dependencies between multiplication operations.
- The vector processing unit (VPU) dedicated to multiplication can process, in a single cycle, only an amount of data equivalent to
$\frac{VLEN}{SEW}$ , that is,$8$ floats.
To mitigate the RAW latency bottleneck, a loop unrolling and instruction interleaving technique was implemented. Instead of processing a single vector block sequentially, the algorithm was restructured to handle four independent blocks simultaneously using vectors with a base register multiplier (LMUL=1).
Load Phase: Four consecutive loads are performed from memory __riscv_vle32, offset by the size of the physical vector (vl_max), initializing four accumulators and four independent terms in separate registers (logical registers v0 through v15).
Computational Interleaving (Latency Hiding): In Taylor's inner loop, multiplication instructions are cross-interleaved. When the instruction from Block 0 is issued and held in the pipeline waiting for its operands to become ready, the out-of-order processor takes advantage of the hardware’s free lanes to immediately issue the instructions from Blocks 1, 2, and 3.
Residue Robustness: To avoid numerical errors or out-of-range violations at the end of the array, the code was divided into two phases: a fast, optimized loop that processes only exact multiples of the combined size of the 4 blocks (4×vlmax), and a sequential vector cleanup loop for the remaining elements (0 to 31 elements). The graph shows how the calculated series compares with the original sine function.
The following table presents the data collected using the Linux performance profiler (perf stat), comparing the base scalar version, flat vector configurations with different LMUL factors, and the optimized version with interleaving
| Metric | Scalar | Vector LMUL=1 |
Vector LMUL=8 |
Vector UNROLLED (4 Blocks) |
|---|---|---|---|---|
| Compute Time (s) | 0.038544 | 0.003781 | 0.003716 | 0.001793 |
| Speedup | Base (1×) | 10.19× | 10.37× | 21.50× |
| LMUL | - | 1 | 8 | 1 |
| Cycles | 618 M | 562 M | 563 M | 561 M |
| Instructions Executed | 682 M | 651 M | 651 M | 650 M |
| IPC | 1.10 | 1.16 | 1.16 | 1.16 |
| L1-DCache Misses | 0.11% | 0.12% | 0.13% | 0.13% |
| Total Elapsed Time (s) | 0.464794 | 0.419387 | 0.431193 | 0.429963 |
Analysis of the hardware counters reveals critical insights into the microarchitecture of the Spacemit K1:
The LMUL Glass Ceiling: In flat vector tests, increasing the LMUL parameter from 1 to 8 did not result in any improvement in computation time (0.003781s vs 0.003716s). This demonstrates that the algorithm’s bottleneck was not caused by the size of the vector register or bandwidth, but rather by the internal latency of sequential execution.
The Success of Interleaving: The unrolled version reduced the computation time to 0.001793 seconds, achieving a speedup of about 2× compared to standard vectorization and 21.50× compared to the original scalar code.
The Global Counters Paradox: When looking at global perf metrics, the process’s total cycles and instructions remain virtually flat. This is because perf stat measures the entire lifecycle of the executable, including the heavy overhead of operating system calls and memory allocation. Pure computation time is cut in half, but the remaining time is dominated almost entirely by the sequential initialization of data in RAM and the rest of the program’s computation.
In conclusion, the instruction interleaving technique exploiting independent intrinsic variables has proven to be an effective strategy for saturating vector functional units, improving instruction-level parallelism (ILP) and enabling OoO in vectorial RISC-V architectures.