ref: 71292834bbafcf6bc409cfd7691f9e6650df025b
dir: /libmath/gemm.c/
#include "lib9.h" #include "mathi.h" void gemm(int transa, int transb, int m, int n, int k, double alpha, double *a, int lda, double *b, int ldb, double beta, double *c, int ldc) { int i1, i2, i3, nota, notb, i, j, jb, jc, l, la; double temp; nota = transa=='N'; notb = transb=='N'; if(m == 0 || n == 0 || (alpha == 0. || k == 0) && beta == 1.){ return; } if(alpha == 0.){ if(beta == 0.){ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = 0.; } } }else{ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = beta * c[i + jc]; } } } return; } if(!a){ if(notb){ /* C := alpha*B + beta*C. */ i1 = n; for(j = 0; j < i1; ++j){ jb = j*ldb; jc = j*ldc; i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = alpha*b[i+jb] + beta*c[i+jc]; } } }else{ /* C := alpha*B' + beta*C. */ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = alpha*b[j+i*ldb] + beta*c[i+jc]; } } } return; } if(notb){ if(nota){ /* Form C := alpha*A*B + beta*C. */ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; if(beta == 0.){ i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = 0.; } }else if(beta != 1.){ i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = beta * c[i + jc]; } } i2 = k; for(l = 0; l < i2; ++l){ la = l*lda; if(b[l + j*ldb] != 0.){ temp = alpha * b[l + j*ldb]; i3 = m; for(i = 0; i < i3; ++i){ c[i + jc] += temp * a[i + la]; } } } } }else{ /* Form C := alpha*A'*B + beta*C */ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; i2 = m; for(i = 0; i < i2; ++i){ temp = 0.; i3 = k; for(l = 0; l < i3; ++l){ temp += a[l + i*lda] * b[l + j*ldb]; } if(beta == 0.){ c[i + jc] = alpha * temp; }else{ c[i + jc] = alpha * temp + beta * c[i + jc]; } } } } }else{ if(nota){ /* Form C := alpha*A*B' + beta*C */ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; if(beta == 0.){ i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = 0.; } }else if(beta != 1.){ i2 = m; for(i = 0; i < i2; ++i){ c[i + jc] = beta * c[i + jc]; } } i2 = k; for(l = 0; l < i2; ++l){ if(b[j + l*ldb] != 0.){ temp = alpha * b[j + l*ldb]; i3 = m; for(i = 0; i < i3; ++i){ c[i + jc] += temp * a[i + l*lda]; } } } } }else{ /* Form C := alpha*A'*B' + beta*C */ i1 = n; for(j = 0; j < i1; ++j){ jc = j*ldc; i2 = m; for(i = 0; i < i2; ++i){ temp = 0.; i3 = k; for(l = 0; l < i3; ++l){ temp += a[l + i*lda] * b[j + l*ldb]; } if(beta == 0.){ c[i + jc] = alpha * temp; }else{ c[i + jc] = alpha * temp + beta * c[i + jc]; } } } } } }