In 2014 I wrote a tutorial on how to write an efficient matrix multiplication that can achieve near peak performance. This is actually not achieved by some Strassen algorithm variant. It's achieved through a cache friendly algorithm (exploiting the cache hierarchy) and SIMD register instructions.

The tutorial showed how to achieve that step by step. Beginning with a naive implementation using three loops. Then using buffers to first cache matrix blocks and then performing the matrix product block wise. That actually means you make copies of blocks and then do computations with the blocks. Sounds stupid but makes sense: Assume you have square N x N blocks then cost of make making a copy is O(N^2) and for the computation it's O(N^3). If the computation with the copied/cached blocks is much faster this actually pays of.

Anyway the details are explained in the tutorial. As the university web server has changed many times after 2014 I moved it to GitHub: https://github.com/michael-lehn/gemm-tutorial

From the tutorial I created ulmBLAS (ulm makes sense, as I work at Ulm University) which is a complete BLAS implementation. Recently I started with a tutorial for the matrix multiplication with Cuda. Currently there is not much text, I just outlined a possible step by step approach to go from naive to a performant implementation (actually the goal was to be on par with NVIDIA's cuBLAS; but it's more on parish ... but I think still gives an idea). It in the cuda subdirectory of ulmBLAS: https://github.com/michael-lehn/ulmBLAS