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