1 #include <cctype>
2 #include <cmath>
3
4 #include "gromacs/utility/real.h"
5 #include "../gmx_blas.h"
6
7 void
F77_FUNC(sgemm,SGEMM)8 F77_FUNC(sgemm,SGEMM)(const char *transa,
9 const char *transb,
10 int *m__,
11 int *n__,
12 int *k__,
13 float *alpha__,
14 float *a,
15 int *lda__,
16 float *b,
17 int *ldb__,
18 float *beta__,
19 float *c,
20 int *ldc__)
21 {
22 const char tra=std::toupper(*transa);
23 const char trb=std::toupper(*transb);
24 float temp;
25 int i,j,l;
26
27 int m = *m__;
28 int n = *n__;
29 int k = *k__;
30 int lda = *lda__;
31 int ldb = *ldb__;
32 int ldc = *ldc__;
33
34 float alpha = *alpha__;
35 float beta = *beta__;
36
37 if(m==0 || n==0 || (( std::abs(alpha)<GMX_FLOAT_MIN || k==0) && std::abs(beta-1.0)<GMX_FLOAT_EPS))
38 return;
39
40 if(std::abs(alpha)<GMX_FLOAT_MIN) {
41 if(std::abs(beta)<GMX_FLOAT_MIN) {
42 for(j=0;j<n;j++)
43 for(i=0;i<m;i++)
44 c[j*(ldc)+i] = 0.0;
45 } else {
46 /* nonzero beta */
47 for(j=0;j<n;j++)
48 for(i=0;i<m;i++)
49 c[j*(ldc)+i] *= beta;
50 }
51 return;
52 }
53
54 if(trb=='N') {
55 if(tra=='N') {
56
57 for(j=0;j<n;j++) {
58 if(std::abs(beta)<GMX_FLOAT_MIN) {
59 for(i=0;i<m;i++)
60 c[j*(ldc)+i] = 0.0;
61 } else if(std::abs(beta-1.0)>GMX_FLOAT_EPS) {
62 for(i=0;i<m;i++)
63 c[j*(ldc)+i] *= beta;
64 }
65 for(l=0;l<k;l++) {
66 if( std::abs(b[ j*(ldb) + l ])>GMX_FLOAT_MIN) {
67 temp = alpha * b[ j*(ldb) + l ];
68 for(i=0;i<m;i++)
69 c[j*(ldc)+i] += temp * a[l*(lda)+i];
70 }
71 }
72 }
73 } else {
74 /* transpose A, but not B */
75 for(j=0;j<n;j++) {
76 for(i=0;i<m;i++) {
77 temp = 0.0;
78 for(l=0;l<k;l++)
79 temp += a[i*(lda)+l] * b[j*(ldb)+l];
80 if(std::abs(beta)<GMX_FLOAT_MIN)
81 c[j*(ldc)+i] = alpha * temp;
82 else
83 c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
84 }
85 }
86 }
87 } else {
88 /* transpose B */
89 if(tra=='N') {
90
91 /* transpose B, but not A */
92
93 for(j=0;j<n;j++) {
94 if(std::abs(beta)<GMX_FLOAT_MIN) {
95 for(i=0;i<m;i++)
96 c[j*(ldc)+i] = 0.0;
97 } else if(std::abs(beta-1.0)>GMX_FLOAT_EPS) {
98 for(i=0;i<m;i++)
99 c[j*(ldc)+i] *= beta;
100 }
101 for(l=0;l<k;l++) {
102 if( std::abs(b[ l*(ldb) + j ])>GMX_FLOAT_MIN) {
103 temp = alpha * b[ l*(ldb) + j ];
104 for(i=0;i<m;i++)
105 c[j*(ldc)+i] += temp * a[l*(lda)+i];
106 }
107 }
108 }
109
110 } else {
111 /* Transpose both A and B */
112 for(j=0;j<n;j++) {
113 for(i=0;i<m;i++) {
114 temp = 0.0;
115 for(l=0;l<k;l++)
116 temp += a[i*(lda)+l] * b[l*(ldb)+j];
117 if(std::abs(beta)<GMX_FLOAT_MIN)
118 c[j*(ldc)+i] = alpha * temp;
119 else
120 c[j*(ldc)+i] = alpha * temp + beta * c[j*(ldc)+i];
121 }
122 }
123 }
124 }
125 }
126