4 #include "gromacs/utility/real.h"
5 #include "../gmx_blas.h"
8 F77_FUNC(ssyr2k
,SSYR2K
)(const char *uplo
,
32 float alpha
= *alpha__
;
36 ch2
= toupper(*trans
);
43 if(n
==0 || ( ( fabs(alpha
)<GMX_FLOAT_MIN
|| k
==0 ) && fabs(beta
-1.0)<GMX_FLOAT_EPS
))
46 if(fabs(alpha
)<GMX_FLOAT_MIN
) {
48 if(fabs(beta
)<GMX_FLOAT_MIN
)
51 c
[(j
-1)*(ldc
)+(i
-1)] = 0.0;
55 c
[(j
-1)*(ldc
)+(i
-1)] *= beta
;
58 if(fabs(beta
)<GMX_FLOAT_MIN
)
61 c
[(j
-1)*(ldc
)+(i
-1)] = 0.0;
65 c
[(j
-1)*(ldc
)+(i
-1)] *= beta
;
73 if(fabs(beta
)<GMX_FLOAT_MIN
)
75 c
[(j
-1)*(ldc
)+(i
-1)] = 0.0;
76 else if(fabs(beta
-1.0)>GMX_FLOAT_EPS
)
78 c
[(j
-1)*(ldc
)+(i
-1)] *= beta
;
80 if( fabs(a
[(l
-1)*(lda
)+(j
-1)])>GMX_FLOAT_MIN
||
81 fabs(b
[(l
-1)*(ldb
)+(j
-1)])>GMX_FLOAT_MIN
) {
82 temp1
= alpha
* b
[(l
-1)*(ldb
)+(j
-1)];
83 temp2
= alpha
* a
[(l
-1)*(lda
)+(j
-1)];
85 c
[(j
-1)*(ldc
)+(i
-1)] +=
86 a
[(l
-1)*(lda
)+(i
-1)] * temp1
+
87 b
[(l
-1)*(ldb
)+(i
-1)] * temp2
;
94 if(fabs(beta
)<GMX_FLOAT_MIN
)
96 c
[(j
-1)*(ldc
)+(i
-1)] = 0.0;
97 else if(fabs(beta
-1.0)>GMX_FLOAT_EPS
)
99 c
[(j
-1)*(ldc
)+(i
-1)] *= beta
;
101 if( fabs(a
[(l
-1)*(lda
)+(j
-1)])>GMX_FLOAT_MIN
||
102 fabs(b
[(l
-1)*(ldb
)+(j
-1)])>GMX_FLOAT_MIN
) {
103 temp1
= alpha
* b
[(l
-1)*(ldb
)+(j
-1)];
104 temp2
= alpha
* a
[(l
-1)*(lda
)+(j
-1)];
106 c
[(j
-1)*(ldc
)+(i
-1)] +=
107 a
[(l
-1)*(lda
)+(i
-1)] * temp1
+
108 b
[(l
-1)*(ldb
)+(i
-1)] * temp2
;
121 temp1
+= a
[(i
-1)*(lda
)+(l
-1)] * b
[(j
-1)*(ldb
)+(l
-1)];
122 temp2
+= b
[(i
-1)*(ldb
)+(l
-1)] * a
[(j
-1)*(lda
)+(l
-1)];
124 if(fabs(beta
)<GMX_FLOAT_MIN
)
125 c
[(j
-1)*(ldc
)+(i
-1)] = alpha
* (temp1
+ temp2
);
127 c
[(j
-1)*(ldc
)+(i
-1)] = beta
* c
[(j
-1)*(ldc
)+(i
-1)] +
128 alpha
* (temp1
+ temp2
);
137 temp1
+= a
[(i
-1)*(lda
)+(l
-1)] * b
[(j
-1)*(ldb
)+(l
-1)];
138 temp2
+= b
[(i
-1)*(ldb
)+(l
-1)] * a
[(j
-1)*(lda
)+(l
-1)];
140 if(fabs(beta
)<GMX_FLOAT_MIN
)
141 c
[(j
-1)*(ldc
)+(i
-1)] = alpha
* (temp1
+ temp2
);
143 c
[(j
-1)*(ldc
)+(i
-1)] = beta
* c
[(j
-1)*(ldc
)+(i
-1)] +
144 alpha
* (temp1
+ temp2
);