Add small marix optimization kernel interface.
authorXianyi Zhang <xianyi@perfxlab.com>
Tue, 28 Apr 2020 11:01:36 +0000 (19:01 +0800)
committerWangyang Guo <wangyang.guo@intel.com>
Mon, 2 Aug 2021 07:01:47 +0000 (07:01 +0000)
make SMALL_MATRIX_OPT=1

Makefile.system
common_d.h
common_level3.h
common_macro.h
common_s.h
interface/gemm.c
kernel/Makefile.L3
kernel/generic/gemm_small_matrix_kernel_nn.c [new file with mode: 0644]
kernel/generic/gemm_small_matrix_kernel_nt.c [new file with mode: 0644]
kernel/generic/gemm_small_matrix_kernel_tn.c [new file with mode: 0644]
kernel/generic/gemm_small_matrix_kernel_tt.c [new file with mode: 0644]

index 13c946b..20d8d2f 100644 (file)
@@ -244,6 +244,11 @@ else
 ONLY_CBLAS = 0
 endif
 
+#For small matrix optimization
+ifeq ($(SMALL_MATRIX_OPT), 1)
+CCOMMON_OPT += -DSMALL_MATRIX_OPT
+endif
+
 # This operation is expensive, so execution should be once.
 ifndef GOTOBLAS_MAKEFILE
 export GOTOBLAS_MAKEFILE = 1
index 94dc3ee..dad304a 100644 (file)
 #define DIMATCOPY_K_RT      dimatcopy_k_rt
 #define DGEADD_K                dgeadd_k 
 
+
+#define DGEMM_SMALL_KERNEL_NN   dgemm_small_kernel_nn
+#define DGEMM_SMALL_KERNEL_NT   dgemm_small_kernel_nt
+#define DGEMM_SMALL_KERNEL_TN   dgemm_small_kernel_tn
+#define DGEMM_SMALL_KERNEL_TT   dgemm_small_kernel_tt
+
 #else
 
 #define        DAMAX_K                 gotoblas -> damax_k
index c4f9435..751592b 100644 (file)
@@ -515,6 +515,18 @@ int qgemm_kernel(BLASLONG, BLASLONG, BLASLONG, xidouble *, xidouble *, xidouble
 int qgemm_kernel(BLASLONG, BLASLONG, BLASLONG, xdouble, xdouble *, xdouble *, xdouble *, BLASLONG);
 #endif
 
+#ifdef SMALL_MATRIX_OPT
+int sgemm_small_kernel_nn(BLASLONG m, BLASLONG n, BLASLONG k, float * A, BLASLONG lda, float alpha, float * B, BLASLONG ldb, float beta, float * C, BLASLONG ldc);
+int sgemm_small_kernel_nt(BLASLONG m, BLASLONG n, BLASLONG k, float * A, BLASLONG lda, float alpha, float * B, BLASLONG ldb, float beta, float * C, BLASLONG ldc);
+int sgemm_small_kernel_tn(BLASLONG m, BLASLONG n, BLASLONG k, float * A, BLASLONG lda, float alpha, float * B, BLASLONG ldb, float beta, float * C, BLASLONG ldc);
+int sgemm_small_kernel_tt(BLASLONG m, BLASLONG n, BLASLONG k, float * A, BLASLONG lda, float alpha, float * B, BLASLONG ldb, float beta, float * C, BLASLONG ldc);
+
+int dgemm_small_kernel_nn(BLASLONG m, BLASLONG n, BLASLONG k, double * A, BLASLONG lda, double alpha, double * B, BLASLONG ldb, double beta, double * C, BLASLONG ldc);
+int dgemm_small_kernel_nt(BLASLONG m, BLASLONG n, BLASLONG k, double * A, BLASLONG lda, double alpha, double * B, BLASLONG ldb, double beta, double * C, BLASLONG ldc);
+int dgemm_small_kernel_tn(BLASLONG m, BLASLONG n, BLASLONG k, double * A, BLASLONG lda, double alpha, double * B, BLASLONG ldb, double beta, double * C, BLASLONG ldc);
+int dgemm_small_kernel_tt(BLASLONG m, BLASLONG n, BLASLONG k, double * A, BLASLONG lda, double alpha, double * B, BLASLONG ldb, double beta, double * C, BLASLONG ldc);
+#endif
+
 int cgemm_kernel_n(BLASLONG, BLASLONG, BLASLONG, float,  float,  float  *, float  *, float  *, BLASLONG);
 int cgemm_kernel_l(BLASLONG, BLASLONG, BLASLONG, float,  float,  float  *, float  *, float  *, BLASLONG);
 int cgemm_kernel_r(BLASLONG, BLASLONG, BLASLONG, float,  float,  float  *, float  *, float  *, BLASLONG);
index 0136f18..eb2abcd 100644 (file)
 
 #define GEADD_K                 DGEADD_K
 
+#define GEMM_SMALL_KERNEL_NN    DGEMM_SMALL_KERNEL_NN
+#define GEMM_SMALL_KERNEL_NT    DGEMM_SMALL_KERNEL_NT
+#define GEMM_SMALL_KERNEL_TN    DGEMM_SMALL_KERNEL_TN
+#define GEMM_SMALL_KERNEL_TT    DGEMM_SMALL_KERNEL_TT
+
 #elif defined(BFLOAT16)
 
 #define D_TO_BF16_K     SBDTOBF16_K
 
 #define GEADD_K                SGEADD_K
 
+#define GEMM_SMALL_KERNEL_NN    SGEMM_SMALL_KERNEL_NN
+#define GEMM_SMALL_KERNEL_NT    SGEMM_SMALL_KERNEL_NT
+#define GEMM_SMALL_KERNEL_TN    SGEMM_SMALL_KERNEL_TN
+#define GEMM_SMALL_KERNEL_TT    SGEMM_SMALL_KERNEL_TT
+
 #endif
 
 #else
 #define IMATCOPY_K_RT          SIMATCOPY_K_RT
 
 #define GEADD_K                SGEADD_K
+
+#define GEMM_SMALL_KERNEL_NN    SGEMM_SMALL_KERNEL_NN
+#define GEMM_SMALL_KERNEL_NT    SGEMM_SMALL_KERNEL_NT
+#define GEMM_SMALL_KERNEL_TN    SGEMM_SMALL_KERNEL_TN
+#define GEMM_SMALL_KERNEL_TT    SGEMM_SMALL_KERNEL_TT
+
 #endif
 #else
 #ifdef XDOUBLE
index 34903ec..6ad98ba 100644 (file)
 
 #define SGEADD_K                sgeadd_k 
 
+#define SGEMM_SMALL_KERNEL_NN   sgemm_small_kernel_nn
+#define SGEMM_SMALL_KERNEL_NT   sgemm_small_kernel_nt
+#define SGEMM_SMALL_KERNEL_TN   sgemm_small_kernel_tn
+#define SGEMM_SMALL_KERNEL_TT   sgemm_small_kernel_tt
+
 #else
 
 #define        SAMAX_K                 gotoblas -> samax_k
index 10426fd..d2fb42f 100644 (file)
@@ -105,6 +105,18 @@ static int (*gemm[])(blas_arg_t *, BLASLONG *, BLASLONG *, IFLOAT *, IFLOAT *, B
 #endif
 };
 
+#ifdef SMALL_MATRIX_OPT
+//Only support s/dgemm small matrix optimiztion so far.
+static int (*gemm_small_kernel[])(BLASLONG, BLASLONG, BLASLONG, FLOAT *, BLASLONG, FLOAT ,FLOAT *, BLASLONG, FLOAT, FLOAT *, BLASLONG) = {
+#ifndef GEMM3M
+#ifndef COMPLEX
+       GEMM_SMALL_KERNEL_NN, GEMM_SMALL_KERNEL_TN, NULL, NULL,
+       GEMM_SMALL_KERNEL_NT, GEMM_SMALL_KERNEL_TT, NULL, NULL,
+#endif
+#endif
+};
+#endif
+
 #ifndef CBLAS
 
 void NAME(char *TRANSA, char *TRANSB,
@@ -417,6 +429,20 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANS
 
   FUNCTION_PROFILE_START();
 
+  MNK = (double) args.m * (double) args.n * (double) args.k;
+
+#ifdef SMALL_MATRIX_OPT
+#if !defined(COMPLEX)
+  //need to tune small matrices cases.
+  if(MNK <= 100.0*100.0*100.0){
+         (gemm_small_kernel[(transb << 2) | transa])(args.m, args.n, args.k, args.a, args.lda, *(FLOAT *)(args.alpha), args.b,
+                                                     args.ldb, *(FLOAT *)(args.beta), args.c, args.ldc);
+         return;
+  }
+#endif
+#endif
+  
+
   buffer = (XFLOAT *)blas_memory_alloc(0);
 
   sa = (XFLOAT *)((BLASLONG)buffer +GEMM_OFFSET_A);
@@ -428,7 +454,7 @@ void CNAME(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE TransA, enum CBLAS_TRANS
   mode |= (transb << BLAS_TRANSB_SHIFT);
 #endif
 
-  MNK = (double) args.m * (double) args.n * (double) args.k;
+
   if ( MNK <= (SMP_THRESHOLD_MIN  * (double) GEMM_MULTITHREAD_THRESHOLD)  )
        args.nthreads = 1;
   else
index 2d9e3ec..88e5eb2 100644 (file)
@@ -447,6 +447,19 @@ XBLASOBJS += \
 
 endif
 
+######  BLAS small matrix optimization #####
+ifeq ($(SMALL_MATRIX_OPT), 1)
+
+SBLASOBJS += \
+       sgemm_small_kernel_nn$(TSUFFIX).$(SUFFIX) sgemm_small_kernel_nt$(TSUFFIX).$(SUFFIX) \
+       sgemm_small_kernel_tn$(TSUFFIX).$(SUFFIX) sgemm_small_kernel_tt$(TSUFFIX).$(SUFFIX)
+
+DBLASOBJS += \
+       dgemm_small_kernel_nn$(TSUFFIX).$(SUFFIX) dgemm_small_kernel_nt$(TSUFFIX).$(SUFFIX) \
+       dgemm_small_kernel_tn$(TSUFFIX).$(SUFFIX) dgemm_small_kernel_tt$(TSUFFIX).$(SUFFIX)
+
+endif
+
 ######  BLAS extensions #####
 
 ifeq ($(BUILD_SINGLE),1)
@@ -4237,3 +4250,63 @@ endif
 $(KDIR)zgeadd_k$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(ZGEADD_K)
        $(CC) $(CFLAGS) -c -DDOUBLE -DCOMPLEX -UROWM $< -o $@
 
+
+
+######  BLAS small matrix optimization #####
+
+ifndef DGEMM_SAMLL_K_NN
+DGEMM_SAMLL_K_NN = ../generic/gemm_small_matrix_kernel_nn.c
+endif
+
+ifndef DGEMM_SAMLL_K_NT
+DGEMM_SAMLL_K_NT = ../generic/gemm_small_matrix_kernel_nt.c
+endif
+
+ifndef DGEMM_SAMLL_K_TN
+DGEMM_SAMLL_K_TN = ../generic/gemm_small_matrix_kernel_tn.c
+endif
+
+ifndef DGEMM_SAMLL_K_TT
+DGEMM_SAMLL_K_TT = ../generic/gemm_small_matrix_kernel_tt.c
+endif
+
+$(KDIR)dgemm_small_kernel_nn$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(DGEMM_SAMLL_K_NN)
+       $(CC) $(CFLAGS) -c -DDOUBLE -UCOMPLEX $< -o $@
+
+$(KDIR)dgemm_small_kernel_nt$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(DGEMM_SAMLL_K_NT)
+       $(CC) $(CFLAGS) -c -DDOUBLE -UCOMPLEX $< -o $@
+
+$(KDIR)dgemm_small_kernel_tn$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(DGEMM_SAMLL_K_TN)
+       $(CC) $(CFLAGS) -c -DDOUBLE -UCOMPLEX $< -o $@
+
+$(KDIR)dgemm_small_kernel_tt$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(DGEMM_SAMLL_K_TT)
+       $(CC) $(CFLAGS) -c -DDOUBLE -UCOMPLEX $< -o $@
+
+
+ifndef SGEMM_SAMLL_K_NN
+SGEMM_SAMLL_K_NN = ../generic/gemm_small_matrix_kernel_nn.c
+endif
+
+ifndef SGEMM_SAMLL_K_NT
+SGEMM_SAMLL_K_NT = ../generic/gemm_small_matrix_kernel_nt.c
+endif
+
+ifndef SGEMM_SAMLL_K_TN
+SGEMM_SAMLL_K_TN = ../generic/gemm_small_matrix_kernel_tn.c
+endif
+
+ifndef SGEMM_SAMLL_K_TT
+SGEMM_SAMLL_K_TT = ../generic/gemm_small_matrix_kernel_tt.c
+endif
+
+$(KDIR)sgemm_small_kernel_nn$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(SGEMM_SAMLL_K_NN)
+       $(CC) $(CFLAGS) -c -UDOUBLE -UCOMPLEX $< -o $@
+
+$(KDIR)sgemm_small_kernel_nt$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(SGEMM_SAMLL_K_NT)
+       $(CC) $(CFLAGS) -c -UDOUBLE -UCOMPLEX $< -o $@
+
+$(KDIR)sgemm_small_kernel_tn$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(SGEMM_SAMLL_K_TN)
+       $(CC) $(CFLAGS) -c -UDOUBLE -UCOMPLEX $< -o $@
+
+$(KDIR)sgemm_small_kernel_tt$(TSUFFIX).$(SUFFIX) : $(KERNELDIR)/$(SGEMM_SAMLL_K_TT)
+       $(CC) $(CFLAGS) -c -UDOUBLE -UCOMPLEX $< -o $@
diff --git a/kernel/generic/gemm_small_matrix_kernel_nn.c b/kernel/generic/gemm_small_matrix_kernel_nn.c
new file mode 100644 (file)
index 0000000..efcc27c
--- /dev/null
@@ -0,0 +1,49 @@
+/***************************************************************************
+Copyright (c) 2020, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+#include "common.h"
+
+int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, FLOAT * A, BLASLONG lda, FLOAT alpha, FLOAT * B, BLASLONG ldb, FLOAT beta, FLOAT * C, BLASLONG ldc)
+{
+       //naive implemtation
+       //Column major
+       
+       BLASLONG i,j,k;
+       FLOAT result=0.0;
+
+       for(i=0; i<M; i++){
+               for(j=0; j<N; j++){
+                       result=0.0;
+                       for(k=0; k<K; k++){
+                               result += A[i+k*lda] * B[k+j*ldb];
+                       }
+                       C[i+j*ldc]=C[i+j*ldc] * beta + alpha * result;
+               }
+       }
+       
+       return 0;
+}
diff --git a/kernel/generic/gemm_small_matrix_kernel_nt.c b/kernel/generic/gemm_small_matrix_kernel_nt.c
new file mode 100644 (file)
index 0000000..e8d9a6c
--- /dev/null
@@ -0,0 +1,49 @@
+/***************************************************************************
+Copyright (c) 2020, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+#include "common.h"
+
+int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, FLOAT * A, BLASLONG lda, FLOAT alpha, FLOAT * B, BLASLONG ldb, FLOAT beta, FLOAT * C, BLASLONG ldc)
+{
+       //naive implemtation
+       //Column major
+
+       BLASLONG i,j,k;
+       FLOAT result=0.0;
+
+       for(i=0; i<M; i++){
+               for(j=0; j<N; j++){
+                       result=0.0;
+                       for(k=0; k<K; k++){
+                               result += A[i+k*lda] * B[k*ldb+j];
+                       }
+                       C[i+j*ldc]=C[i+j*ldc] * beta + alpha * result;
+               }
+       }
+       
+       return 0;
+}
diff --git a/kernel/generic/gemm_small_matrix_kernel_tn.c b/kernel/generic/gemm_small_matrix_kernel_tn.c
new file mode 100644 (file)
index 0000000..f7b7f2f
--- /dev/null
@@ -0,0 +1,49 @@
+/***************************************************************************
+Copyright (c) 2020, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+#include "common.h"
+
+int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, FLOAT * A, BLASLONG lda, FLOAT alpha, FLOAT * B, BLASLONG ldb, FLOAT beta, FLOAT * C, BLASLONG ldc)
+{
+       //naive implemtation
+       //Column major
+
+       BLASLONG i,j,k;
+       FLOAT result=0.0;
+
+       for(i=0; i<M; i++){
+               for(j=0; j<N; j++){
+                       result=0.0;
+                       for(k=0; k<K; k++){
+                               result += A[i*lda+k] * B[k+j*ldb];
+                       }
+                       C[i+j*ldc]=C[i+j*ldc] * beta + alpha * result;
+               }
+       }
+
+       return 0;
+}
diff --git a/kernel/generic/gemm_small_matrix_kernel_tt.c b/kernel/generic/gemm_small_matrix_kernel_tt.c
new file mode 100644 (file)
index 0000000..40a5b1b
--- /dev/null
@@ -0,0 +1,49 @@
+/***************************************************************************
+Copyright (c) 2020, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+#include "common.h"
+
+int CNAME(BLASLONG M, BLASLONG N, BLASLONG K, FLOAT * A, BLASLONG lda, FLOAT alpha, FLOAT * B, BLASLONG ldb, FLOAT beta, FLOAT * C, BLASLONG ldc)
+{
+       //naive implemtation
+       //Column major
+
+       BLASLONG i,j,k;
+       FLOAT result=0.0;
+
+       for(i=0; i<M; i++){
+               for(j=0; j<N; j++){
+                       result=0.0;
+                       for(k=0; k<K; k++){
+                               result += A[i*lda+k] * B[k*ldb+j];
+                       }
+                       C[i+j*ldc]=C[i+j*ldc] * beta + alpha * result;
+               }
+       }
+
+       return 0;
+}