I wonder why the cascaded triple for loop for GEMM is not just hardcoded in assembly for each architecture. And while we are at it, why is MKL > 5GB again?
Most decent BLAS implementations use some variation of the kernels in K. Goto et al. "Anatomy of High-Performance Matrix Multiplication" written in assembly.
Does this address your MKL question?[1] If so, sounds like we're two of today's lucky 10,000![0] Although for this esoterica, I'd probably reduce that to 10.