1 /* dsbmv.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
dsbmv_(char * uplo,integer * n,integer * k,doublereal * alpha,doublereal * a,integer * lda,doublereal * x,integer * incx,doublereal * beta,doublereal * y,integer * incy,ftnlen uplo_len)15 /* Subroutine */ int dsbmv_(char *uplo, integer *n, integer *k, doublereal *
16 alpha, doublereal *a, integer *lda, doublereal *x, integer *incx,
17 doublereal *beta, doublereal *y, integer *incy, ftnlen uplo_len)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
21
22 /* Local variables */
23 integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
24 doublereal temp1, temp2;
25 extern logical lsame_(char *, char *, ftnlen, ftnlen);
26 integer kplus1;
27 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
28
29 /* .. Scalar Arguments .. */
30 /* .. */
31 /* .. Array Arguments .. */
32 /* .. */
33
34 /* Purpose */
35 /* ======= */
36
37 /* DSBMV performs the matrix-vector operation */
38
39 /* y := alpha*A*x + beta*y, */
40
41 /* where alpha and beta are scalars, x and y are n element vectors and */
42 /* A is an n by n symmetric band matrix, with k super-diagonals. */
43
44 /* Arguments */
45 /* ========== */
46
47 /* UPLO - CHARACTER*1. */
48 /* On entry, UPLO specifies whether the upper or lower */
49 /* triangular part of the band matrix A is being supplied as */
50 /* follows: */
51
52 /* UPLO = 'U' or 'u' The upper triangular part of A is */
53 /* being supplied. */
54
55 /* UPLO = 'L' or 'l' The lower triangular part of A is */
56 /* being supplied. */
57
58 /* Unchanged on exit. */
59
60 /* N - INTEGER. */
61 /* On entry, N specifies the order of the matrix A. */
62 /* N must be at least zero. */
63 /* Unchanged on exit. */
64
65 /* K - INTEGER. */
66 /* On entry, K specifies the number of super-diagonals of the */
67 /* matrix A. K must satisfy 0 .le. K. */
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 with UPLO = 'U' or 'u', the leading ( k + 1 ) */
76 /* by n part of the array A must contain the upper triangular */
77 /* band part of the symmetric matrix, supplied column by */
78 /* column, with the leading diagonal of the matrix in row */
79 /* ( k + 1 ) of the array, the first super-diagonal starting at */
80 /* position 2 in row k, and so on. The top left k by k triangle */
81 /* of the array A is not referenced. */
82 /* The following program segment will transfer the upper */
83 /* triangular part of a symmetric band matrix from conventional */
84 /* full matrix storage to band storage: */
85
86 /* DO 20, J = 1, N */
87 /* M = K + 1 - J */
88 /* DO 10, I = MAX( 1, J - K ), J */
89 /* A( M + I, J ) = matrix( I, J ) */
90 /* 10 CONTINUE */
91 /* 20 CONTINUE */
92
93 /* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
94 /* by n part of the array A must contain the lower triangular */
95 /* band part of the symmetric matrix, supplied column by */
96 /* column, with the leading diagonal of the matrix in row 1 of */
97 /* the array, the first sub-diagonal starting at position 1 in */
98 /* row 2, and so on. The bottom right k by k triangle of the */
99 /* array A is not referenced. */
100 /* The following program segment will transfer the lower */
101 /* triangular part of a symmetric band matrix from conventional */
102 /* full matrix storage to band storage: */
103
104 /* DO 20, J = 1, N */
105 /* M = 1 - J */
106 /* DO 10, I = J, MIN( N, J + K ) */
107 /* A( M + I, J ) = matrix( I, J ) */
108 /* 10 CONTINUE */
109 /* 20 CONTINUE */
110
111 /* Unchanged on exit. */
112
113 /* LDA - INTEGER. */
114 /* On entry, LDA specifies the first dimension of A as declared */
115 /* in the calling (sub) program. LDA must be at least */
116 /* ( k + 1 ). */
117 /* Unchanged on exit. */
118
119 /* X - DOUBLE PRECISION array of DIMENSION at least */
120 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
121 /* Before entry, the incremented array X must contain the */
122 /* vector x. */
123 /* Unchanged on exit. */
124
125 /* INCX - INTEGER. */
126 /* On entry, INCX specifies the increment for the elements of */
127 /* X. INCX must not be zero. */
128 /* Unchanged on exit. */
129
130 /* BETA - DOUBLE PRECISION. */
131 /* On entry, BETA specifies the scalar beta. */
132 /* Unchanged on exit. */
133
134 /* Y - DOUBLE PRECISION array of DIMENSION at least */
135 /* ( 1 + ( n - 1 )*abs( INCY ) ). */
136 /* Before entry, the incremented array Y must contain the */
137 /* vector y. On exit, Y is overwritten by the updated vector y. */
138
139 /* INCY - INTEGER. */
140 /* On entry, INCY specifies the increment for the elements of */
141 /* Y. INCY must not be zero. */
142 /* Unchanged on exit. */
143
144
145 /* Level 2 Blas routine. */
146
147 /* -- Written on 22-October-1986. */
148 /* Jack Dongarra, Argonne National Lab. */
149 /* Jeremy Du Croz, Nag Central Office. */
150 /* Sven Hammarling, Nag Central Office. */
151 /* Richard Hanson, Sandia National Labs. */
152
153 /* ===================================================================== */
154
155 /* .. Parameters .. */
156 /* .. */
157 /* .. Local Scalars .. */
158 /* .. */
159 /* .. External Functions .. */
160 /* .. */
161 /* .. External Subroutines .. */
162 /* .. */
163 /* .. Intrinsic Functions .. */
164 /* .. */
165
166 /* Test the input parameters. */
167
168 /* Parameter adjustments */
169 a_dim1 = *lda;
170 a_offset = 1 + a_dim1;
171 a -= a_offset;
172 --x;
173 --y;
174
175 /* Function Body */
176 info = 0;
177 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
178 ftnlen)1, (ftnlen)1)) {
179 info = 1;
180 } else if (*n < 0) {
181 info = 2;
182 } else if (*k < 0) {
183 info = 3;
184 } else if (*lda < *k + 1) {
185 info = 6;
186 } else if (*incx == 0) {
187 info = 8;
188 } else if (*incy == 0) {
189 info = 11;
190 }
191 if (info != 0) {
192 xerbla_("DSBMV ", &info, (ftnlen)6);
193 return 0;
194 }
195
196 /* Quick return if possible. */
197
198 if (*n == 0 || (*alpha == 0. && *beta == 1.)) {
199 return 0;
200 }
201
202 /* Set up the start points in X and Y. */
203
204 if (*incx > 0) {
205 kx = 1;
206 } else {
207 kx = 1 - (*n - 1) * *incx;
208 }
209 if (*incy > 0) {
210 ky = 1;
211 } else {
212 ky = 1 - (*n - 1) * *incy;
213 }
214
215 /* Start the operations. In this version the elements of the array A */
216 /* are accessed sequentially with one pass through A. */
217
218 /* First form y := beta*y. */
219
220 if (*beta != 1.) {
221 if (*incy == 1) {
222 if (*beta == 0.) {
223 i__1 = *n;
224 for (i__ = 1; i__ <= i__1; ++i__) {
225 y[i__] = 0.;
226 /* L10: */
227 }
228 } else {
229 i__1 = *n;
230 for (i__ = 1; i__ <= i__1; ++i__) {
231 y[i__] = *beta * y[i__];
232 /* L20: */
233 }
234 }
235 } else {
236 iy = ky;
237 if (*beta == 0.) {
238 i__1 = *n;
239 for (i__ = 1; i__ <= i__1; ++i__) {
240 y[iy] = 0.;
241 iy += *incy;
242 /* L30: */
243 }
244 } else {
245 i__1 = *n;
246 for (i__ = 1; i__ <= i__1; ++i__) {
247 y[iy] = *beta * y[iy];
248 iy += *incy;
249 /* L40: */
250 }
251 }
252 }
253 }
254 if (*alpha == 0.) {
255 return 0;
256 }
257 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
258
259 /* Form y when upper triangle of A is stored. */
260
261 kplus1 = *k + 1;
262 if (*incx == 1 && *incy == 1) {
263 i__1 = *n;
264 for (j = 1; j <= i__1; ++j) {
265 temp1 = *alpha * x[j];
266 temp2 = 0.;
267 l = kplus1 - j;
268 /* Computing MAX */
269 i__2 = 1, i__3 = j - *k;
270 i__4 = j - 1;
271 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
272 y[i__] += temp1 * a[l + i__ + j * a_dim1];
273 temp2 += a[l + i__ + j * a_dim1] * x[i__];
274 /* L50: */
275 }
276 y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
277 /* L60: */
278 }
279 } else {
280 jx = kx;
281 jy = ky;
282 i__1 = *n;
283 for (j = 1; j <= i__1; ++j) {
284 temp1 = *alpha * x[jx];
285 temp2 = 0.;
286 ix = kx;
287 iy = ky;
288 l = kplus1 - j;
289 /* Computing MAX */
290 i__4 = 1, i__2 = j - *k;
291 i__3 = j - 1;
292 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
293 y[iy] += temp1 * a[l + i__ + j * a_dim1];
294 temp2 += a[l + i__ + j * a_dim1] * x[ix];
295 ix += *incx;
296 iy += *incy;
297 /* L70: */
298 }
299 y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha *
300 temp2;
301 jx += *incx;
302 jy += *incy;
303 if (j > *k) {
304 kx += *incx;
305 ky += *incy;
306 }
307 /* L80: */
308 }
309 }
310 } else {
311
312 /* Form y when lower triangle of A is stored. */
313
314 if (*incx == 1 && *incy == 1) {
315 i__1 = *n;
316 for (j = 1; j <= i__1; ++j) {
317 temp1 = *alpha * x[j];
318 temp2 = 0.;
319 y[j] += temp1 * a[j * a_dim1 + 1];
320 l = 1 - j;
321 /* Computing MIN */
322 i__4 = *n, i__2 = j + *k;
323 i__3 = min(i__4,i__2);
324 for (i__ = j + 1; i__ <= i__3; ++i__) {
325 y[i__] += temp1 * a[l + i__ + j * a_dim1];
326 temp2 += a[l + i__ + j * a_dim1] * x[i__];
327 /* L90: */
328 }
329 y[j] += *alpha * temp2;
330 /* L100: */
331 }
332 } else {
333 jx = kx;
334 jy = ky;
335 i__1 = *n;
336 for (j = 1; j <= i__1; ++j) {
337 temp1 = *alpha * x[jx];
338 temp2 = 0.;
339 y[jy] += temp1 * a[j * a_dim1 + 1];
340 l = 1 - j;
341 ix = jx;
342 iy = jy;
343 /* Computing MIN */
344 i__4 = *n, i__2 = j + *k;
345 i__3 = min(i__4,i__2);
346 for (i__ = j + 1; i__ <= i__3; ++i__) {
347 ix += *incx;
348 iy += *incy;
349 y[iy] += temp1 * a[l + i__ + j * a_dim1];
350 temp2 += a[l + i__ + j * a_dim1] * x[ix];
351 /* L110: */
352 }
353 y[jy] += *alpha * temp2;
354 jx += *incx;
355 jy += *incy;
356 /* L120: */
357 }
358 }
359 }
360
361 return 0;
362
363 /* End of DSBMV . */
364
365 } /* dsbmv_ */
366
367