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
sgemv_(char * trans,integer * m,integer * n,real * alpha,real * a,integer * lda,real * x,integer * incx,real * beta,real * y,integer * incy)9 /* Subroutine */ int sgemv_(char *trans, integer *m, integer *n, real *alpha,
10 real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
11 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 real 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 SGEMV 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 - REAL .
63 On entry, ALPHA specifies the scalar alpha.
64 Unchanged on exit.
65
66 A - REAL 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 - REAL 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 - REAL .
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 - REAL 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
120
121 Test the input parameters.
122
123
124 Parameter adjustments
125 Function Body */
126 #define X(I) x[(I)-1]
127 #define Y(I) y[(I)-1]
128
129 #define A(I,J) a[(I)-1 + ((J)-1)* ( *lda)]
130
131 info = 0;
132 if ( strncmp(trans, "N", 1)!=0 && strncmp(trans, "T", 1)!=0 &&
133 strncmp(trans, "C", 1)!=0 ) {
134 info = 1;
135 } else if (*m < 0) {
136 info = 2;
137 } else if (*n < 0) {
138 info = 3;
139 } else if (*lda < max(1,*m)) {
140 info = 6;
141 } else if (*incx == 0) {
142 info = 8;
143 } else if (*incy == 0) {
144 info = 11;
145 }
146 if (info != 0) {
147 input_error("SGEMV ", &info);
148 return 0;
149 }
150
151 /* Quick return if possible. */
152
153 if (*m == 0 || *n == 0 || *alpha == 0.f && *beta == 1.f) {
154 return 0;
155 }
156
157 /* Set LENX and LENY, the lengths of the vectors x and y, and set
158
159 up the start points in X and Y. */
160
161 if (strncmp(trans, "N", 1)==0) {
162 lenx = *n;
163 leny = *m;
164 } else {
165 lenx = *m;
166 leny = *n;
167 }
168 if (*incx > 0) {
169 kx = 1;
170 } else {
171 kx = 1 - (lenx - 1) * *incx;
172 }
173 if (*incy > 0) {
174 ky = 1;
175 } else {
176 ky = 1 - (leny - 1) * *incy;
177 }
178
179 /* Start the operations. In this version the elements of A are
180 accessed sequentially with one pass through A.
181
182 First form y := beta*y. */
183
184 if (*beta != 1.f) {
185 if (*incy == 1) {
186 if (*beta == 0.f) {
187 i__1 = leny;
188 for (i = 1; i <= leny; ++i) {
189 Y(i) = 0.f;
190 /* L10: */
191 }
192 } else {
193 i__1 = leny;
194 for (i = 1; i <= leny; ++i) {
195 Y(i) = *beta * Y(i);
196 /* L20: */
197 }
198 }
199 } else {
200 iy = ky;
201 if (*beta == 0.f) {
202 i__1 = leny;
203 for (i = 1; i <= leny; ++i) {
204 Y(iy) = 0.f;
205 iy += *incy;
206 /* L30: */
207 }
208 } else {
209 i__1 = leny;
210 for (i = 1; i <= leny; ++i) {
211 Y(iy) = *beta * Y(iy);
212 iy += *incy;
213 /* L40: */
214 }
215 }
216 }
217 }
218 if (*alpha == 0.f) {
219 return 0;
220 }
221 if (strncmp(trans, "N", 1)==0) {
222
223 /* Form y := alpha*A*x + y. */
224
225 jx = kx;
226 if (*incy == 1) {
227 i__1 = *n;
228 for (j = 1; j <= *n; ++j) {
229 if (X(jx) != 0.f) {
230 temp = *alpha * X(jx);
231 i__2 = *m;
232 for (i = 1; i <= *m; ++i) {
233 Y(i) += temp * A(i,j);
234 /* L50: */
235 }
236 }
237 jx += *incx;
238 /* L60: */
239 }
240 } else {
241 i__1 = *n;
242 for (j = 1; j <= *n; ++j) {
243 if (X(jx) != 0.f) {
244 temp = *alpha * X(jx);
245 iy = ky;
246 i__2 = *m;
247 for (i = 1; i <= *m; ++i) {
248 Y(iy) += temp * A(i,j);
249 iy += *incy;
250 /* L70: */
251 }
252 }
253 jx += *incx;
254 /* L80: */
255 }
256 }
257 } else {
258
259 /* Form y := alpha*A'*x + y. */
260
261 jy = ky;
262 if (*incx == 1) {
263 i__1 = *n;
264 for (j = 1; j <= *n; ++j) {
265 temp = 0.f;
266 i__2 = *m;
267 for (i = 1; i <= *m; ++i) {
268 temp += A(i,j) * X(i);
269 /* L90: */
270 }
271 Y(jy) += *alpha * temp;
272 jy += *incy;
273 /* L100: */
274 }
275 } else {
276 i__1 = *n;
277 for (j = 1; j <= *n; ++j) {
278 temp = 0.f;
279 ix = kx;
280 i__2 = *m;
281 for (i = 1; i <= *m; ++i) {
282 temp += A(i,j) * X(ix);
283 ix += *incx;
284 /* L110: */
285 }
286 Y(jy) += *alpha * temp;
287 jy += *incy;
288 /* L120: */
289 }
290 }
291 }
292
293 return 0;
294
295 /* End of SGEMV . */
296
297 } /* sgemv_ */
298
299