1 /* ssbmv.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
ssbmv_(char * uplo,integer * n,integer * k,real * alpha,real * a,integer * lda,real * x,integer * incx,real * beta,real * y,integer * incy,ftnlen uplo_len)15 /* Subroutine */ int ssbmv_(char *uplo, integer *n, integer *k, real *alpha,
16 real *a, integer *lda, real *x, integer *incx, real *beta, real *y,
17 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 real 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 /* SSBMV 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 - REAL . */
71 /* On entry, ALPHA specifies the scalar alpha. */
72 /* Unchanged on exit. */
73
74 /* A - REAL 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 - REAL 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 - REAL . */
131 /* On entry, BETA specifies the scalar beta. */
132 /* Unchanged on exit. */
133
134 /* Y - REAL 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 /* Further Details */
145 /* =============== */
146
147 /* Level 2 Blas routine. */
148
149 /* -- Written on 22-October-1986. */
150 /* Jack Dongarra, Argonne National Lab. */
151 /* Jeremy Du Croz, Nag Central Office. */
152 /* Sven Hammarling, Nag Central Office. */
153 /* Richard Hanson, Sandia National Labs. */
154
155 /* ===================================================================== */
156
157 /* .. Parameters .. */
158 /* .. */
159 /* .. Local Scalars .. */
160 /* .. */
161 /* .. External Functions .. */
162 /* .. */
163 /* .. External Subroutines .. */
164 /* .. */
165 /* .. Intrinsic Functions .. */
166 /* .. */
167
168 /* Test the input parameters. */
169
170 /* Parameter adjustments */
171 a_dim1 = *lda;
172 a_offset = 1 + a_dim1;
173 a -= a_offset;
174 --x;
175 --y;
176
177 /* Function Body */
178 info = 0;
179 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
180 ftnlen)1, (ftnlen)1)) {
181 info = 1;
182 } else if (*n < 0) {
183 info = 2;
184 } else if (*k < 0) {
185 info = 3;
186 } else if (*lda < *k + 1) {
187 info = 6;
188 } else if (*incx == 0) {
189 info = 8;
190 } else if (*incy == 0) {
191 info = 11;
192 }
193 if (info != 0) {
194 xerbla_("SSBMV ", &info, (ftnlen)6);
195 return 0;
196 }
197
198 /* Quick return if possible. */
199
200 if (*n == 0 || (*alpha == 0.f && *beta == 1.f)) {
201 return 0;
202 }
203
204 /* Set up the start points in X and Y. */
205
206 if (*incx > 0) {
207 kx = 1;
208 } else {
209 kx = 1 - (*n - 1) * *incx;
210 }
211 if (*incy > 0) {
212 ky = 1;
213 } else {
214 ky = 1 - (*n - 1) * *incy;
215 }
216
217 /* Start the operations. In this version the elements of the array A */
218 /* are accessed sequentially with one pass through A. */
219
220 /* First form y := beta*y. */
221
222 if (*beta != 1.f) {
223 if (*incy == 1) {
224 if (*beta == 0.f) {
225 i__1 = *n;
226 for (i__ = 1; i__ <= i__1; ++i__) {
227 y[i__] = 0.f;
228 /* L10: */
229 }
230 } else {
231 i__1 = *n;
232 for (i__ = 1; i__ <= i__1; ++i__) {
233 y[i__] = *beta * y[i__];
234 /* L20: */
235 }
236 }
237 } else {
238 iy = ky;
239 if (*beta == 0.f) {
240 i__1 = *n;
241 for (i__ = 1; i__ <= i__1; ++i__) {
242 y[iy] = 0.f;
243 iy += *incy;
244 /* L30: */
245 }
246 } else {
247 i__1 = *n;
248 for (i__ = 1; i__ <= i__1; ++i__) {
249 y[iy] = *beta * y[iy];
250 iy += *incy;
251 /* L40: */
252 }
253 }
254 }
255 }
256 if (*alpha == 0.f) {
257 return 0;
258 }
259 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
260
261 /* Form y when upper triangle of A is stored. */
262
263 kplus1 = *k + 1;
264 if (*incx == 1 && *incy == 1) {
265 i__1 = *n;
266 for (j = 1; j <= i__1; ++j) {
267 temp1 = *alpha * x[j];
268 temp2 = 0.f;
269 l = kplus1 - j;
270 /* Computing MAX */
271 i__2 = 1, i__3 = j - *k;
272 i__4 = j - 1;
273 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
274 y[i__] += temp1 * a[l + i__ + j * a_dim1];
275 temp2 += a[l + i__ + j * a_dim1] * x[i__];
276 /* L50: */
277 }
278 y[j] = y[j] + temp1 * a[kplus1 + j * a_dim1] + *alpha * temp2;
279 /* L60: */
280 }
281 } else {
282 jx = kx;
283 jy = ky;
284 i__1 = *n;
285 for (j = 1; j <= i__1; ++j) {
286 temp1 = *alpha * x[jx];
287 temp2 = 0.f;
288 ix = kx;
289 iy = ky;
290 l = kplus1 - j;
291 /* Computing MAX */
292 i__4 = 1, i__2 = j - *k;
293 i__3 = j - 1;
294 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
295 y[iy] += temp1 * a[l + i__ + j * a_dim1];
296 temp2 += a[l + i__ + j * a_dim1] * x[ix];
297 ix += *incx;
298 iy += *incy;
299 /* L70: */
300 }
301 y[jy] = y[jy] + temp1 * a[kplus1 + j * a_dim1] + *alpha *
302 temp2;
303 jx += *incx;
304 jy += *incy;
305 if (j > *k) {
306 kx += *incx;
307 ky += *incy;
308 }
309 /* L80: */
310 }
311 }
312 } else {
313
314 /* Form y when lower triangle of A is stored. */
315
316 if (*incx == 1 && *incy == 1) {
317 i__1 = *n;
318 for (j = 1; j <= i__1; ++j) {
319 temp1 = *alpha * x[j];
320 temp2 = 0.f;
321 y[j] += temp1 * a[j * a_dim1 + 1];
322 l = 1 - j;
323 /* Computing MIN */
324 i__4 = *n, i__2 = j + *k;
325 i__3 = min(i__4,i__2);
326 for (i__ = j + 1; i__ <= i__3; ++i__) {
327 y[i__] += temp1 * a[l + i__ + j * a_dim1];
328 temp2 += a[l + i__ + j * a_dim1] * x[i__];
329 /* L90: */
330 }
331 y[j] += *alpha * temp2;
332 /* L100: */
333 }
334 } else {
335 jx = kx;
336 jy = ky;
337 i__1 = *n;
338 for (j = 1; j <= i__1; ++j) {
339 temp1 = *alpha * x[jx];
340 temp2 = 0.f;
341 y[jy] += temp1 * a[j * a_dim1 + 1];
342 l = 1 - j;
343 ix = jx;
344 iy = jy;
345 /* Computing MIN */
346 i__4 = *n, i__2 = j + *k;
347 i__3 = min(i__4,i__2);
348 for (i__ = j + 1; i__ <= i__3; ++i__) {
349 ix += *incx;
350 iy += *incy;
351 y[iy] += temp1 * a[l + i__ + j * a_dim1];
352 temp2 += a[l + i__ + j * a_dim1] * x[ix];
353 /* L110: */
354 }
355 y[jy] += *alpha * temp2;
356 jx += *incx;
357 jy += *incy;
358 /* L120: */
359 }
360 }
361 }
362
363 return 0;
364
365 /* End of SSBMV . */
366
367 } /* ssbmv_ */
368
369