[Paper Review] A GEMM interface and implementation on NVIDIA GPUs for multiple small matrices

Chetan Jhurani, Paul Mullowney|arXiv (Cornell University)|Apr 26, 2013

Parallel Computing and Optimization Techniques5 references4 citations

TL;DR

This paper presents a high-performance GEMM interface and CUDA kernel implementation for batched multiplication of small matrices (≤16×16) on NVIDIA GPUs, leveraging a second-level leading dimension interface to reduce memory overhead and improve data locality. The implementation achieves 30% to 600% higher performance than cuBLAS batched GEMM on Tesla K20c, with up to 216 GFlop/s for single-precision real matrices of size 16.

ABSTRACT

We present an interface and an implementation of the General Matrix Multiply (GEMM) routine for multiple small matrices processed simultaneously on NVIDIA graphics processing units (GPUs). We focus on matrix sizes under 16. The implementation can be easily extended to larger sizes. For single precision matrices, our implementation is 30% to 600% faster than the batched cuBLAS implementation distributed in the CUDA Toolkit 5.0 on NVIDIA Tesla K20c. For example, we obtain 104 GFlop/s and 216 GFlop/s when multiplying 100,000 independent matrix pairs of size 10 and 16, respectively. Similar improvement in performance is obtained for other sizes, in single and double precision for real and complex types, and when the number of matrices is smaller. Apart from our implementation, our different function interface also plays an important role in the improved performance. Applications of this software include Finite Element computation on GPUs.

Motivation & Objective

• Address the performance bottleneck in GEMM operations for small matrices (≤16×16) on GPUs, which are common in finite element methods.
• Overcome limitations of the cuBLAS batched interface, particularly the performance degradation when matrix sizes are not multiples of 16 and the overhead from pointer-to-pointer argument passing.
• Design a new interface using second-level leading dimensions to improve data locality and reduce memory transfer overhead.
• Achieve high performance across multiple data types (single/double precision, real/complex) and matrix sizes, with minimal performance drop for non-multiples of 16.
• Enable efficient, portable, and high-throughput GEMM for thousands of independent small matrix operations on modern GPUs.

Proposed method

• Design a new GEMM interface based on second-level leading dimensions instead of double-pointer arrays, reducing memory transfer and improving data coalescing.
• Implement a CUDA kernel using thread blocks optimized for small matrix sizes, with coalesced memory access patterns and efficient register usage.
• Leverage C++ templates to support multiple data types (S/D/C/Z) and matrix sizes in a single, type-safe, and efficient implementation.
• Use a uniform memory layout for batches of matrices, enabling coalesced global memory access and reducing bank conflicts in shared memory.
• Optimize kernel launch parameters and grid/block dimensions to maximize occupancy and hide memory latency.
• Integrate the kernel with a C++ wrapper function that supports various operation types (transpose, conjugate transpose) and scalar parameters (α, β).

Experimental results

Research questions

• RQ1 How can the performance of batched GEMM for small matrices be significantly improved on NVIDIA GPUs compared to existing cuBLAS batched implementations?
• RQ2 To what extent does a redesigned interface using second-level leading dimensions reduce memory transfer overhead and improve data locality compared to pointer-to-pointer interfaces?
• RQ3 What performance gains are achievable across different matrix sizes (1–16) and data types (single/double, real/complex) when using the proposed interface and kernel?
• RQ4 How does the performance of the new implementation scale with increasing batch size, and does it maintain high efficiency for non-power-of-16 sizes?
• RQ5 Can the proposed interface and kernel design be generalized to other BLAS routines and hardware platforms?

Key findings

• The proposed implementation achieves up to 600% higher performance than the cuBLAS batched GEMM for single-precision real matrices of size 10, reaching 104 GFlop/s for 100,000 matrix pairs.
• For matrix size 16, the implementation attains 216 GFlop/s (single-precision real), 173 GFlop/s (double-precision real), 609 GFlop/s (single-precision complex), and 217 GFlop/s (double-precision complex) on Tesla K20c.
• Performance degradation due to non-multiples of 16 is significantly reduced: for size 15, the new kernel achieves 150 GFlop/s (vs. 105 GFlop/s in cuBLAS), and for size 17, it achieves 150 GFlop/s (vs. 32 GFlop/s in cuBLAS).
• The interface change alone (using second-level leading dimension) contributes significantly to performance, with relative improvements up to 600% for small matrices.
• The performance gain is consistent across all data types and matrix sizes, with improvements ranging from 30% to 600% depending on size and type.
• The implementation remains efficient even for smaller batch sizes and maintains high performance for complex and double-precision types, with 504 GFlop/s for single-precision complex matrices of size 15.