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