diff options
Diffstat (limited to 'libmath/gemm.c')
| -rw-r--r-- | libmath/gemm.c | 167 |
1 files changed, 167 insertions, 0 deletions
diff --git a/libmath/gemm.c b/libmath/gemm.c new file mode 100644 index 00000000..908a0998 --- /dev/null +++ b/libmath/gemm.c @@ -0,0 +1,167 @@ +#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]; + } + } + } + } + } +} |