1
2 /* -- translated by f2c (version 19940927).
3 You must link the resulting object file with the libraries:
4 -lf2c -lm (in that order)
5 */
6 #include <string.h>
7 #include "f2c.h"
8
dgemv_(char * trans,integer * m,integer * n,doublereal * alpha,doublereal * a,integer * lda,doublereal * x,integer * incx,doublereal * beta,doublereal * y,integer * incy)9 /* Subroutine */ int dgemv_(char *trans, integer *m, integer *n, doublereal *
10 alpha, doublereal *a, integer *lda, doublereal *x, integer *incx,
11 doublereal *beta, doublereal *y, integer *incy)
12 {
13
14
15 /* System generated locals */
16 integer a_dim1, a_offset, i__1, i__2;
17
18 /* Local variables */
19 static integer info;
20 static doublereal temp;
21 static integer lenx, leny, i, j;
22 static integer ix, iy, jx, jy, kx, ky;
23
24 extern int input_error(char *, int *);
25
26 /* Purpose
27 =======
28
29 DGEMV performs one of the matrix-vector operations
30
31 y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
32
33 where alpha and beta are scalars, x and y are vectors and A is an
34 m by n matrix.
35
36 Parameters
37 ==========
38
39 TRANS - CHARACTER*1.
40 On entry, TRANS specifies the operation to be performed as
41 follows:
42
43 TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
44
45 TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
46
47 TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
48
49 Unchanged on exit.
50
51 M - INTEGER.
52 On entry, M specifies the number of rows of the matrix A.
53 M must be at least zero.
54 Unchanged on exit.
55
56 N - INTEGER.
57 On entry, N specifies the number of columns of the matrix A.
58
59 N must be at least zero.
60 Unchanged on exit.
61
62 ALPHA - DOUBLE PRECISION.
63 On entry, ALPHA specifies the scalar alpha.
64 Unchanged on exit.
65
66 A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
67 Before entry, the leading m by n part of the array A must
68 contain the matrix of coefficients.
69 Unchanged on exit.
70
71 LDA - INTEGER.
72 On entry, LDA specifies the first dimension of A as declared
73
74 in the calling (sub) program. LDA must be at least
75 max( 1, m ).
76 Unchanged on exit.
77
78 X - DOUBLE PRECISION array of DIMENSION at least
79 ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
80 and at least
81 ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
82 Before entry, the incremented array X must contain the
83 vector x.
84 Unchanged on exit.
85
86 INCX - INTEGER.
87 On entry, INCX specifies the increment for the elements of
88 X. INCX must not be zero.
89 Unchanged on exit.
90
91 BETA - DOUBLE PRECISION.
92 On entry, BETA specifies the scalar beta. When BETA is
93 supplied as zero then Y need not be set on input.
94 Unchanged on exit.
95
96 Y - DOUBLE PRECISION array of DIMENSION at least
97 ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
98 and at least
99 ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
100 Before entry with BETA non-zero, the incremented array Y
101 must contain the vector y. On exit, Y is overwritten by the
102
103 updated vector y.
104
105 INCY - INTEGER.
106 On entry, INCY specifies the increment for the elements of
107 Y. INCY must not be zero.
108 Unchanged on exit.
109
110
111 Level 2 Blas routine.
112
113 -- Written on 22-October-1986.
114 Jack Dongarra, Argonne National Lab.
115 Jeremy Du Croz, Nag Central Office.
116 Sven Hammarling, Nag Central Office.
117 Richard Hanson, Sandia National Labs.
118
119 Test the input parameters.
120
121 Parameter adjustments
122 Function Body */
123 #define X(I) x[(I)-1]
124 #define Y(I) y[(I)-1]
125
126 #define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
127
128 info = 0;
129 if ( strncmp(trans, "N", 1)!=0 && strncmp(trans, "T", 1)!=0 &&
130 strncmp(trans, "C", 1)==0 ) {
131 info = 1;
132 } else if (*m < 0) {
133 info = 2;
134 } else if (*n < 0) {
135 info = 3;
136 } else if (*lda < max(1,*m)) {
137 info = 6;
138 } else if (*incx == 0) {
139 info = 8;
140 } else if (*incy == 0) {
141 info = 11;
142 }
143 if (info != 0) {
144 input_error("DGEMV ", &info);
145 return 0;
146 }
147
148 /* Quick return if possible. */
149
150 if (*m == 0 || *n == 0 || *alpha == 0. && *beta == 1.) {
151 return 0;
152 }
153
154 /* Set LENX and LENY, the lengths of the vectors x and y, and set
155
156 up the start points in X and Y. */
157
158 if (strncmp(trans, "N", 1)==0) {
159 lenx = *n;
160 leny = *m;
161 } else {
162 lenx = *m;
163 leny = *n;
164 }
165 if (*incx > 0) {
166 kx = 1;
167 } else {
168 kx = 1 - (lenx - 1) * *incx;
169 }
170 if (*incy > 0) {
171 ky = 1;
172 } else {
173 ky = 1 - (leny - 1) * *incy;
174 }
175
176 /* Start the operations. In this version the elements of A are
177 accessed sequentially with one pass through A.
178
179 First form y := beta*y. */
180
181 if (*beta != 1.) {
182 if (*incy == 1) {
183 if (*beta == 0.) {
184 i__1 = leny;
185 for (i = 1; i <= leny; ++i) {
186 Y(i) = 0.;
187 /* L10: */
188 }
189 } else {
190 i__1 = leny;
191 for (i = 1; i <= leny; ++i) {
192 Y(i) = *beta * Y(i);
193 /* L20: */
194 }
195 }
196 } else {
197 iy = ky;
198 if (*beta == 0.) {
199 i__1 = leny;
200 for (i = 1; i <= leny; ++i) {
201 Y(iy) = 0.;
202 iy += *incy;
203 /* L30: */
204 }
205 } else {
206 i__1 = leny;
207 for (i = 1; i <= leny; ++i) {
208 Y(iy) = *beta * Y(iy);
209 iy += *incy;
210 /* L40: */
211 }
212 }
213 }
214 }
215 if (*alpha == 0.) {
216 return 0;
217 }
218 if (strncmp(trans, "N", 1)==0) {
219
220 /* Form y := alpha*A*x + y. */
221
222 jx = kx;
223 if (*incy == 1) {
224 i__1 = *n;
225 for (j = 1; j <= *n; ++j) {
226 if (X(jx) != 0.) {
227 temp = *alpha * X(jx);
228 i__2 = *m;
229 for (i = 1; i <= *m; ++i) {
230 Y(i) += temp * A(i,j);
231 /* L50: */
232 }
233 }
234 jx += *incx;
235 /* L60: */
236 }
237 } else {
238 i__1 = *n;
239 for (j = 1; j <= *n; ++j) {
240 if (X(jx) != 0.) {
241 temp = *alpha * X(jx);
242 iy = ky;
243 i__2 = *m;
244 for (i = 1; i <= *m; ++i) {
245 Y(iy) += temp * A(i,j);
246 iy += *incy;
247 /* L70: */
248 }
249 }
250 jx += *incx;
251 /* L80: */
252 }
253 }
254 } else {
255
256 /* Form y := alpha*A'*x + y. */
257
258 jy = ky;
259 if (*incx == 1) {
260 i__1 = *n;
261 for (j = 1; j <= *n; ++j) {
262 temp = 0.;
263 i__2 = *m;
264 for (i = 1; i <= *m; ++i) {
265 temp += A(i,j) * X(i);
266 /* L90: */
267 }
268 Y(jy) += *alpha * temp;
269 jy += *incy;
270 /* L100: */
271 }
272 } else {
273 i__1 = *n;
274 for (j = 1; j <= *n; ++j) {
275 temp = 0.;
276 ix = kx;
277 i__2 = *m;
278 for (i = 1; i <= *m; ++i) {
279 temp += A(i,j) * X(ix);
280 ix += *incx;
281 /* L110: */
282 }
283 Y(jy) += *alpha * temp;
284 jy += *incy;
285 /* L120: */
286 }
287 }
288 }
289
290 return 0;
291
292 /* End of DGEMV . */
293
294 } /* dgemv_ */
295
296