Add a "sgemm direct" mode for small matrixes

OpenBLAS has a fancy algorithm for copying the input data while laying
it out in a more CPU friendly memory layout.

This is great for large matrixes; the cost of the copy is easily
ammortized by the gains from the better memory layout.

But for small matrixes (on CPUs that can do efficient unaligned loads) this
copy can be a net loss.

This patch adds (for SKYLAKEX initially) a "sgemm direct" mode, that bypasses
the whole copy machinary for ALPHA=1/BETA=0/... standard arguments,
for small matrixes only.

What is small? For the non-threaded case this has been measured to be
in the M*N*K = 28 * 512 * 512 range, while in the threaded case it's
less, around M*N*K = 1 * 512 * 512

Author: fenrus75

common_level3.h

@@ -47,6 +47,14 @@ __global__ void cuda_dgemm_kernel(int, int, int, double *, double *, double *);
 extern "C" {
 #endif
 
+extern void sgemm_kernel_direct(BLASLONG M, BLASLONG N, BLASLONG K,
+	float * __restrict__ A, BLASLONG strideA,
+	float * __restrict__ B, BLASLONG strideB,
+	float * __restrict__ R, BLASLONG strideR);
+
+extern int sgemm_kernel_direct_performant(BLASLONG M, BLASLONG N, BLASLONG K);
+
+
 int sgemm_beta(BLASLONG, BLASLONG, BLASLONG, float,
 	       float  *, BLASLONG, float   *, BLASLONG, float  *, BLASLONG);
 int dgemm_beta(BLASLONG, BLASLONG, BLASLONG, double,

interface/gemm.c

@@ -271,6 +271,14 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANS
 
   PRINT_DEBUG_CNAME;
 
+#if !defined(COMPLEX) && !defined(DOUBLE) && defined(USE_SGEMM_KERNEL_DIRECT)
+  if (beta == 0 && alpha == 1.0 && order == CblasRowMajor && TransA == CblasNoTrans && TransB == CblasNoTrans && sgemm_kernel_direct_performant(m,n,k)) {
+	sgemm_kernel_direct(m, n, k, a, lda, b, ldb, c, ldc);
+	return;
+  }
+
+#endif
+
 #ifndef COMPLEX
   args.alpha = (void *)α
   args.beta  = (void *)β