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