1 #include <cctype>
2 #include <cmath>
3
4 #include "gromacs/utility/real.h"
5 #include "../gmx_blas.h"
6
7 void
F77_FUNC(ssymv,SSYMV)8 F77_FUNC(ssymv,SSYMV)(const char *uplo,
9 int *n__,
10 float *alpha__,
11 float *a,
12 int *lda__,
13 float *x,
14 int *incx__,
15 float *beta__,
16 float *y,
17 int *incy__)
18 {
19 const char ch=std::toupper(*uplo);
20 int kx,ky,i,j,ix,iy,jx,jy;
21 float temp1,temp2;
22
23 int n = *n__;
24 int lda = *lda__;
25 int incx = *incx__;
26 int incy = *incy__;
27 float alpha = *alpha__;
28 float beta = *beta__;
29
30 if(n<=0 || incx==0 || incy==0)
31 return;
32
33 if(incx>0)
34 kx = 1;
35 else
36 kx = 1 - (n -1)*(incx);
37
38 if(incy>0)
39 ky = 1;
40 else
41 ky = 1 - (n -1)*(incy);
42
43 if(std::abs(beta-1.0)>GMX_FLOAT_EPS) {
44 if(incy==1) {
45 if(std::abs(beta)<GMX_FLOAT_MIN)
46 for(i=1;i<=n;i++)
47 y[i-1] = 0.0;
48 else
49 for(i=1;i<=n;i++)
50 y[i-1] *= beta;
51 } else {
52 /* non-unit incr. */
53 iy = ky;
54 if(std::abs(beta)<GMX_FLOAT_MIN)
55 for(i=1;i<=n;i++) {
56 y[iy-1] = 0.0;
57 iy += incy;
58 }
59 else
60 for(i=1;i<=n;i++) {
61 y[iy-1] *= beta;
62 iy += incy;
63 }
64 }
65 }
66
67 if(std::abs(alpha)<GMX_FLOAT_MIN)
68 return;
69
70 if(ch=='U') {
71 if(incx==1 && incy==1) {
72 for(j=1;j<=n;j++) {
73 temp1 = alpha * x[j-1];
74 temp2 = 0.0;
75 for(i=1;i<j;i++) {
76 y[i-1] += temp1*a[(j-1)*(lda)+(i-1)];
77 temp2 += a[(j-1)*(lda)+(i-1)] * x[i-1];
78 }
79 y[j-1] += temp1*a[(j-1)*(lda)+(j-1)] + alpha *temp2;
80 }
81 } else {
82 /* non-unit incr. */
83 jx = kx;
84 jy = ky;
85 for(j=1;j<=n;j++) {
86 temp1 = alpha * x[jx-1];
87 temp2 = 0.0;
88 ix = kx;
89 iy = ky;
90 for(i=1;i<j;i++) {
91 y[iy-1] += temp1 * a[(j-1)*(lda)+(i-1)];
92 temp2 += a[(j-1)*(lda)+(i-1)] * x[ix-1];
93 ix += incx;
94 iy += incy;
95 }
96 y[jy-1] += temp1*a[(j-1)*(lda)+(j-1)] + alpha*temp2;
97 jx += incx;
98 jy += incy;
99 }
100 }
101 } else {
102 /* lower */
103 if(incx==1 && incy==1) {
104 for(j=1;j<=n;j++) {
105 temp1 = alpha * x[j-1];
106 temp2 = 0.0;
107 y[j-1] += temp1 * a[(j-1)*(lda)+(j-1)];
108 for(i=j+1;i<=n;i++) {
109 y[i-1] += temp1*a[(j-1)*(lda)+(i-1)];
110 temp2 += a[(j-1)*(lda)+(i-1)] * x[i-1];
111 }
112 y[j-1] += alpha *temp2;
113 }
114 } else {
115 /* non-unit incr. */
116 jx = kx;
117 jy = ky;
118 for(j=1;j<=n;j++) {
119 temp1 = alpha * x[jx-1];
120 temp2 = 0.0;
121 y[jy-1] += temp1 * a[(j-1)*(lda)+(j-1)];
122 ix = jx;
123 iy = jy;
124 for(i=j+1;i<=n;i++) {
125 ix += incx;
126 iy += incy;
127 y[iy-1] += temp1 * a[(j-1)*(lda)+(i-1)];
128 temp2 += a[(j-1)*(lda)+(i-1)] * x[ix-1];
129 }
130 y[jy-1] += alpha*temp2;
131 jx += incx;
132 jy += incy;
133 }
134 }
135 }
136 return;
137 }
138