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