how to optimize matrix multiplication (matmul) code to run fast on a single processor core

by

The state-of-the-art implementation of matrix multiplication on CPUs uses GotoBLAS algorithm. Basically the loops are organized in the following order:

Loop5 for jc = 0 to N-1 in steps of NC
Loop4   for kc = 0 to K-1 in steps of KC
          //Pack KCxNC block of B
Loop3     for ic = 0 to M-1 in steps of MC
            //Pack MCxKC block of A
//--------------------Macro Kernel------------
Loop2       for jr = 0 to NC-1 in steps of NR
Loop1         for ir = 0 to MC-1 in steps of MR
//--------------------Micro Kernel------------
Loop0           for k = 0 to KC-1 in steps of 1
                //update MRxNR block of C matrix

A key insight underlying modern high-performance implementations of matrix multiplication is to organize the computations by partitioning the operands into blocks for temporal locality (3 outer most loops), and to pack (copy) such blocks
into contiguous buffers that fit into various levels of memory for spatial locality (3 inner most loops).

GotoBLAS algorithm

The above figure (originally from this paper, directly used in this tutorial) illustrates the GotoBLAS algorithm as implemented in BLIS. Cache blocking parameters {MC, NC, KC} determine
the submatrix sizes of Bp (KC × NC) and Ai (MC × KC), such that they fit in various caches. During the computation, row panels Bp
are contiguously packed into buffer Bp to fit in the L3 cache. Blocks Ai are similarly packed into buffer Ai
to fit in the L2 cache. Register block sizes {MR, NR} relate to submatrices in registers that contribute to C. In the micro-kernel (the inner most loop), a small MR × NR micro-tile of C is updated by pair of MR × KC and KC × NR slivers of Ai and Bp.

For the Strassen’s algorithm with O(N^2.87) complexity, you might be interested in reading this paper. Other fast matrix multiplication algorithms with asymptotic complexity less than O(N^3) can be easily extended in this paper. There is a recent thesis about the practical fast matrix multiplication algorithms.

The following tutorials might be helpful if you want to learn more about how to optimize matrix multiplication on CPUs:

How to Optimize GEMM Wiki

GEMM: From Pure C to SSE Optimized Micro Kernels

BLISlab: A sandbox for optimizing GEMM for CPU and ARM

A most updated document about how to optimize GEMM on CPUs (with AVX2/FMA) step by step can be downloaded here:
https://github.com/ULAFF/LAFF-On-HPC/blob/master/LAFF-On-PfHP.pdf

A Massive Open Online Course to be offered on edX starting in June 2019 (LAFF-On Programming for High Performance):
https://github.com/ULAFF/LAFF-On-HPC
http://www.cs.utexas.edu/users/flame/laff/pfhp/LAFF-On-PfHP.html

More Related Contents:

  1. Optimization code for a debug function (c++11 with templates) [closed]
  2. What are copy elision and return value optimization?
  3. Which, if any, C++ compilers do tail-recursion optimization?
  4. What are the mechanics of short string optimization in libc++?
  5. Why is transposing a matrix of 512×512 much slower than transposing a matrix of 513×513?
  6. Is it better to use std::memcpy() or std::copy() in terms to performance?
  7. gcc optimization flag -O3 makes code slower than -O2
  8. Why does optimisation kill this function?
  9. Constants and compiler optimization in C++
  10. Fastest way to get IPv4 address from string
  11. How to get IOStream to perform better?
  12. while (1) Vs. for (;;) Is there a speed difference?
  13. Which is better option to use for dividing an integer number by 2?
  14. Why does Clang optimize away x * 1.0 but NOT x + 0.0?
  15. When do compilers inline C++ code?
  16. Why should I ever use inline code?
  17. How performing multiple matrix multiplications in CUDA?
  18. C++ error: ‘_mm_sin_ps’ was not declared in this scope
  19. How does OpenMP handle nested loops?
  20. Fast multiplication/division by 2 for floats and doubles (C/C++)
  21. Loop with a zero execution time
  22. std::vector reserve() and push_back() is faster than resize() and array index, why?
  23. Optimizations for pow() with const non-integer exponent?
  24. Optimizing member variable order in C++
  25. Does GCC inline C++ functions without the ‘inline’ keyword?
  26. Which is faster/preferred: memset or for loop to zero out an array of doubles?
  27. Different execution policies at runtime
  28. VS: unexpected optimization behavior with _BitScanReverse64 intrinsic
  29. Why is the != operator not allowed with OpenMP?
  30. List of common C++ Optimization Techniques [closed]

Leave a Comment