gemmt.c

/*********************************************************************/
/* Copyright 2022, 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. */
/* */
/* THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF TEXAS AT */
/* AUSTIN ``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 UNIVERSITY OF TEXAS AT */
/* AUSTIN 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<stdio.h>
#include<stdlib.h>
#include"common.h"
#ifdefFUNCTION_PROFILE
#include"functable.h"
#endif

#ifndefCOMPLEX
#defineSMP_THRESHOLD_MIN 65536.0
#ifdefXDOUBLE
#defineERROR_NAME "QGEMT "
#elif defined(DOUBLE)
#defineERROR_NAME "DGEMT "
#elif defined(BFLOAT16)
#defineERROR_NAME "SBGEMT "
#else
#defineERROR_NAME "SGEMT "
#endif
#else
#defineSMP_THRESHOLD_MIN 8192.0
#ifdefXDOUBLE
#defineERROR_NAME "XGEMT "
#elif defined(DOUBLE)
#defineERROR_NAME "ZGEMT "
#else
#defineERROR_NAME "CGEMT "
#endif
#endif

#ifndefGEMM_MULTITHREAD_THRESHOLD
#defineGEMM_MULTITHREAD_THRESHOLD 4
#endif

#ifndefCBLAS

voidNAME(char*UPLO, char*TRANSA, char*TRANSB,
blasint*M, blasint*N, blasint*K,
FLOAT*Alpha,
IFLOAT*a, blasint*ldA,
IFLOAT*b, blasint*ldB, FLOAT*Beta, FLOAT*c, blasint*ldC)
{

blasintm, n, k;
blasintlda, ldb, ldc;
inttransa, transb, uplo;
blasintinfo;

chartransA, transB, Uplo;
IFLOAT*buffer;
IFLOAT*aa, *bb;
FLOAT*cc;
#if defined(COMPLEX)
FLOATalpha_r, alpha_i, beta_r, beta_i;
#else
FLOATalpha, beta;
#endif

PRINT_DEBUG_NAME;

m=*M;
n=*N;
k=*K;

#if defined(COMPLEX)
FLOAT*alpha=Alpha;
alpha_r=*(Alpha+0);
alpha_i=*(Alpha+1);

beta_r=*(Beta+0);
beta_i=*(Beta+1);
#else
alpha=*Alpha;
beta=*Beta;
#endif

lda=*ldA;
ldb=*ldB;
ldc=*ldC;

transA=*TRANSA;
transB=*TRANSB;
Uplo=*UPLO;
TOUPPER(transA);
TOUPPER(transB);
TOUPPER(Uplo);

transa=-1;
transb=-1;
uplo=-1;

if (transA=='N')
transa=0;
if (transA=='T')
transa=1;
#ifndefCOMPLEX
if (transA=='R')
transa=0;
if (transA=='C')
transa=1;
#else
if (transA=='R')
transa=2;
if (transA=='C')
transa=3;
#endif

if (transB=='N')
transb=0;
if (transB=='T')
transb=1;
#ifndefCOMPLEX
if (transB=='R')
transb=0;
if (transB=='C')
transb=1;
#else
if (transB=='R')
transb=2;
if (transB=='C')
transb=3;
#endif

if (Uplo=='U')
uplo=0;
if (Uplo=='L')
uplo=1;

info=0;

if (uplo<0)
info=14;
if (ldc<m)
info=13;
if (k<0)
info=5;
if (n<0)
info=4;
if (m<0)
info=3;
if (transb<0)
info=2;
if (transa<0)
info=1;

if (info) {
BLASFUNC(xerbla) (ERROR_NAME, &info, sizeof(ERROR_NAME));
return;
 }
#else

voidCNAME(enumCBLAS_ORDERorder, enumCBLAS_UPLOUplo,
enumCBLAS_TRANSPOSETransA, enumCBLAS_TRANSPOSETransB, blasintM,
blasintN, blasintk,
#ifndefCOMPLEX
FLOATalpha,
IFLOAT*A, blasintLDA,
IFLOAT*B, blasintLDB, FLOATbeta, FLOAT*c, blasintldc)
{
#else
void*valpha,
void*va, blasintLDA,
void*vb, blasintLDB, void*vbeta, void*vc, blasintldc)
{
FLOAT*alpha= (FLOAT*) valpha;
FLOAT*beta= (FLOAT*) vbeta;
FLOAT*A= (FLOAT*) va;
FLOAT*B= (FLOAT*) vb;
FLOAT*c= (FLOAT*) vc;
#endif
FLOAT*aa, *bb, *cc;

inttransa, transb, uplo;
blasintinfo;
blasintm, n, lda, ldb;
FLOAT*a, *b;
XFLOAT*buffer;

PRINT_DEBUG_CNAME;

transa=-1;
transb=-1;
info=0;

if (order==CblasColMajor) {

if (TransA==CblasNoTrans)
transa=0;
if (TransA==CblasTrans)
transa=1;
#ifndefCOMPLEX
if (TransA==CblasConjNoTrans)
transa=0;
if (TransA==CblasConjTrans)
transa=1;
#else
if (TransA==CblasConjNoTrans)
transa=2;
if (TransA==CblasConjTrans)
transa=3;
#endif
if (TransB==CblasNoTrans)
transb=0;
if (TransB==CblasTrans)
transb=1;
#ifndefCOMPLEX
if (TransB==CblasConjNoTrans)
transb=0;
if (TransB==CblasConjTrans)
transb=1;
#else
if (TransB==CblasConjNoTrans)
transb=2;
if (TransB==CblasConjTrans)
transb=3;
#endif

m=M;
n=N;

a= (void*)A;
b= (void*)B;
lda=LDA;
ldb=LDB;

info=-1;

if (ldc<m)
info=13;
if (k<0)
info=5;
if (n<0)
info=4;
if (m<0)
info=3;
if (transb<0)
info=2;
if (transa<0)
info=1;
 }

if (order==CblasRowMajor) {
m=N;
n=M;

a= (void*)B;
b= (void*)A;

lda=LDB;
ldb=LDA;

if (TransB==CblasNoTrans)
transa=0;
if (TransB==CblasTrans)
transa=1;
#ifndefCOMPLEX
if (TransB==CblasConjNoTrans)
transa=0;
if (TransB==CblasConjTrans)
transa=1;
#else
if (TransB==CblasConjNoTrans)
transa=2;
if (TransB==CblasConjTrans)
transa=3;
#endif
if (TransA==CblasNoTrans)
transb=0;
if (TransA==CblasTrans)
transb=1;
#ifndefCOMPLEX
if (TransA==CblasConjNoTrans)
transb=0;
if (TransA==CblasConjTrans)
transb=1;
#else
if (TransA==CblasConjNoTrans)
transb=2;
if (TransA==CblasConjTrans)
transb=3;
#endif

info=-1;

if (ldc<m)
info=13;
if (k<0)
info=5;
if (n<0)
info=4;
if (m<0)
info=3;
if (transb<0)
info=2;
if (transa<0)
info=1;

 }

uplo=-1;
if (Uplo==CblasUpper)
uplo=0;
if (Uplo==CblasLower)
uplo=1;
if (uplo<0)
info=14;

if (info >= 0) {
BLASFUNC(xerbla) (ERROR_NAME, &info, sizeof(ERROR_NAME));
return;
 }
#if defined(COMPLEX)
FLOATalpha_r=*(alpha+0);
FLOATalpha_i=*(alpha+1);

FLOATbeta_r=*(beta+0);
FLOATbeta_i=*(beta+1);
#endif

#endif
intbuffer_size;
blasintl;
blasinti, j;

#ifdefSMP
intnthreads;
#endif

#if defined(COMPLEX)

#ifdefSMP
staticint (*gemv_thread[]) (BLASLONG, BLASLONG, FLOAT*, FLOAT*,
BLASLONG, FLOAT*, BLASLONG, FLOAT*,
BLASLONG, FLOAT*, int) = {
#ifdefXDOUBLE
xgemv_thread_n, xgemv_thread_t, xgemv_thread_r, xgemv_thread_c,
xgemv_thread_o, xgemv_thread_u, xgemv_thread_s,
xgemv_thread_d,
#elif defined DOUBLE
zgemv_thread_n, zgemv_thread_t, zgemv_thread_r, zgemv_thread_c,
zgemv_thread_o, zgemv_thread_u, zgemv_thread_s,
zgemv_thread_d,
#else
cgemv_thread_n, cgemv_thread_t, cgemv_thread_r, cgemv_thread_c,
cgemv_thread_o, cgemv_thread_u, cgemv_thread_s,
cgemv_thread_d,
#endif
 };
#endif

int (*gemv[]) (BLASLONG, BLASLONG, BLASLONG, FLOAT, FLOAT, FLOAT*,
BLASLONG, FLOAT*, BLASLONG, FLOAT*, BLASLONG,
FLOAT*) = {
GEMV_N, GEMV_T, GEMV_R, GEMV_C, GEMV_O, GEMV_U, GEMV_S, GEMV_D,};

#else

#ifdefSMP
staticint (*gemv_thread[]) (BLASLONG, BLASLONG, FLOAT, FLOAT*,
BLASLONG, FLOAT*, BLASLONG, FLOAT*,
BLASLONG, FLOAT*, int) = {
#ifdefXDOUBLE
qgemv_thread_n, qgemv_thread_t,
#elif defined DOUBLE
dgemv_thread_n, dgemv_thread_t,
#else
sgemv_thread_n, sgemv_thread_t,
#endif
 };
#endif
int (*gemv[]) (BLASLONG, BLASLONG, BLASLONG, FLOAT, FLOAT*, BLASLONG,
FLOAT*, BLASLONG, FLOAT*, BLASLONG, FLOAT*) = {
GEMV_N, GEMV_T,};

#endif

if ((m==0) || (n==0))
return;

IDEBUG_START;

FUNCTION_PROFILE_START();

constblasintincb= (transb==0) ? 1 : ldb;

if (uplo==1) {
for (i=0; i<n; i++) {
j=n-i;

l=j;
#if defined(COMPLEX)
aa=a+i*2;
bb=b+i*ldb*2;
if (transa) {
l=k;
aa=a+lda*i*2;
bb=b+i*2;
 }
cc=c+i*2*ldc+i*2;
#else
aa=a+i;
bb=b+i*ldb;
if (transa) {
l=k;
aa=a+lda*i;
bb=b+i;
 }
cc=c+i*ldc+i;
#endif

#if defined(COMPLEX)
if (beta_r!=ONE||beta_i!=ZERO)
SCAL_K(l, 0, 0, beta_r, beta_i, cc, 1, NULL, 0,
NULL, 0);

if (alpha_r==ZERO&&alpha_i==ZERO)
return;
#else
if (beta!=ONE)
SCAL_K(l, 0, 0, beta, cc, 1, NULL, 0, NULL, 0);

if (alpha==ZERO)
continue;
#endif

IDEBUG_START;

FUNCTION_PROFILE_START();

buffer_size=j+k+128 / sizeof(FLOAT);
#ifdefWINDOWS_ABI
buffer_size+=160 / sizeof(FLOAT);
#endif
// for alignment
buffer_size= (buffer_size+3) & ~3;
STACK_ALLOC(buffer_size, FLOAT, buffer);

#ifdefSMP

if (1L*j*k<2304L*GEMM_MULTITHREAD_THRESHOLD)
nthreads=1;
else
nthreads=num_cpu_avail(2);

if (nthreads==1) {
#endif

#if defined(COMPLEX)
 (gemv[(int)transa]) (j, k, 0, alpha_r, alpha_i,
aa, lda, bb, incb, cc, 1,
buffer);
#else
 (gemv[(int)transa]) (j, k, 0, alpha, aa, lda,
bb, incb, cc, 1, buffer);
#endif
#ifdefSMP
 } else {

 (gemv_thread[(int)transa]) (j, k, alpha, aa,
lda, bb, incb, cc,
1, buffer,
nthreads);

 }
#endif

STACK_FREE(buffer);
 }
 } else {

for (i=0; i<n; i++) {
j=i+1;

l=j;
#if defined COMPLEX
bb=b+i*ldb*2;
if (transa) {
l=k;
bb=b+i*2;
 }
cc=c+i*2*ldc;
#else
bb=b+i*ldb;
if (transa) {
l=k;
bb=b+i;
 }
cc=c+i*ldc;
#endif

#if defined(COMPLEX)
if (beta_r!=ONE||beta_i!=ZERO)
SCAL_K(l, 0, 0, beta_r, beta_i, cc, 1, NULL, 0,
NULL, 0);

if (alpha_r==ZERO&&alpha_i==ZERO)
return;
#else
if (beta!=ONE)
SCAL_K(l, 0, 0, beta, cc, 1, NULL, 0, NULL, 0);

if (alpha==ZERO)
continue;
#endif
IDEBUG_START;

FUNCTION_PROFILE_START();

buffer_size=j+k+128 / sizeof(FLOAT);
#ifdefWINDOWS_ABI
buffer_size+=160 / sizeof(FLOAT);
#endif
// for alignment
buffer_size= (buffer_size+3) & ~3;
STACK_ALLOC(buffer_size, FLOAT, buffer);

#ifdefSMP

if (1L*j*k<2304L*GEMM_MULTITHREAD_THRESHOLD)
nthreads=1;
else
nthreads=num_cpu_avail(2);

if (nthreads==1) {
#endif

#if defined(COMPLEX)
 (gemv[(int)transa]) (j, k, 0, alpha_r, alpha_i,
a, lda, bb, incb, cc, 1,
buffer);
#else
 (gemv[(int)transa]) (j, k, 0, alpha, a, lda, bb,
incb, cc, 1, buffer);
#endif

#ifdefSMP
 } else {

 (gemv_thread[(int)transa]) (j, k, alpha, a, lda,
bb, incb, cc, 1,
buffer, nthreads);

 }
#endif

STACK_FREE(buffer);
 }
 }
FUNCTION_PROFILE_END(COMPSIZE*COMPSIZE,
args.m*args.k+args.k*args.n+
args.m*args.n, 2*args.m*args.n*args.k);

IDEBUG_END;

return;
}