1 /* dtbmv.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
dtbmv_(char * uplo,char * trans,char * diag,integer * n,integer * k,doublereal * a,integer * lda,doublereal * x,integer * incx,ftnlen uplo_len,ftnlen trans_len,ftnlen diag_len)15 /* Subroutine */ int dtbmv_(char *uplo, char *trans, char *diag, integer *n,
16 integer *k, doublereal *a, integer *lda, doublereal *x, integer *incx,
17 ftnlen uplo_len, ftnlen trans_len, ftnlen diag_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, jx, kx, info;
24 doublereal temp;
25 extern logical lsame_(char *, char *, ftnlen, ftnlen);
26 integer kplus1;
27 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
28 logical nounit;
29
30 /* .. Scalar Arguments .. */
31 /* .. */
32 /* .. Array Arguments .. */
33 /* .. */
34
35 /* Purpose */
36 /* ======= */
37
38 /* DTBMV performs one of the matrix-vector operations */
39
40 /* x := A*x, or x := A'*x, */
41
42 /* where x is an n element vector and A is an n by n unit, or non-unit, */
43 /* upper or lower triangular band matrix, with ( k + 1 ) diagonals. */
44
45 /* Arguments */
46 /* ========== */
47
48 /* UPLO - CHARACTER*1. */
49 /* On entry, UPLO specifies whether the matrix is an upper or */
50 /* lower triangular matrix as follows: */
51
52 /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
53
54 /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
55
56 /* Unchanged on exit. */
57
58 /* TRANS - CHARACTER*1. */
59 /* On entry, TRANS specifies the operation to be performed as */
60 /* follows: */
61
62 /* TRANS = 'N' or 'n' x := A*x. */
63
64 /* TRANS = 'T' or 't' x := A'*x. */
65
66 /* TRANS = 'C' or 'c' x := A'*x. */
67
68 /* Unchanged on exit. */
69
70 /* DIAG - CHARACTER*1. */
71 /* On entry, DIAG specifies whether or not A is unit */
72 /* triangular as follows: */
73
74 /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
75
76 /* DIAG = 'N' or 'n' A is not assumed to be unit */
77 /* triangular. */
78
79 /* Unchanged on exit. */
80
81 /* N - INTEGER. */
82 /* On entry, N specifies the order of the matrix A. */
83 /* N must be at least zero. */
84 /* Unchanged on exit. */
85
86 /* K - INTEGER. */
87 /* On entry with UPLO = 'U' or 'u', K specifies the number of */
88 /* super-diagonals of the matrix A. */
89 /* On entry with UPLO = 'L' or 'l', K specifies the number of */
90 /* sub-diagonals of the matrix A. */
91 /* K must satisfy 0 .le. K. */
92 /* Unchanged on exit. */
93
94 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
95 /* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
96 /* by n part of the array A must contain the upper triangular */
97 /* band part of the matrix of coefficients, supplied column by */
98 /* column, with the leading diagonal of the matrix in row */
99 /* ( k + 1 ) of the array, the first super-diagonal starting at */
100 /* position 2 in row k, and so on. The top left k by k triangle */
101 /* of the array A is not referenced. */
102 /* The following program segment will transfer an upper */
103 /* triangular band matrix from conventional full matrix storage */
104 /* to band storage: */
105
106 /* DO 20, J = 1, N */
107 /* M = K + 1 - J */
108 /* DO 10, I = MAX( 1, J - K ), J */
109 /* A( M + I, J ) = matrix( I, J ) */
110 /* 10 CONTINUE */
111 /* 20 CONTINUE */
112
113 /* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
114 /* by n part of the array A must contain the lower triangular */
115 /* band part of the matrix of coefficients, supplied column by */
116 /* column, with the leading diagonal of the matrix in row 1 of */
117 /* the array, the first sub-diagonal starting at position 1 in */
118 /* row 2, and so on. The bottom right k by k triangle of the */
119 /* array A is not referenced. */
120 /* The following program segment will transfer a lower */
121 /* triangular band matrix from conventional full matrix storage */
122 /* to band storage: */
123
124 /* DO 20, J = 1, N */
125 /* M = 1 - J */
126 /* DO 10, I = J, MIN( N, J + K ) */
127 /* A( M + I, J ) = matrix( I, J ) */
128 /* 10 CONTINUE */
129 /* 20 CONTINUE */
130
131 /* Note that when DIAG = 'U' or 'u' the elements of the array A */
132 /* corresponding to the diagonal elements of the matrix are not */
133 /* referenced, but are assumed to be unity. */
134 /* Unchanged on exit. */
135
136 /* LDA - INTEGER. */
137 /* On entry, LDA specifies the first dimension of A as declared */
138 /* in the calling (sub) program. LDA must be at least */
139 /* ( k + 1 ). */
140 /* Unchanged on exit. */
141
142 /* X - DOUBLE PRECISION array of dimension at least */
143 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
144 /* Before entry, the incremented array X must contain the n */
145 /* element vector x. On exit, X is overwritten with the */
146 /* tranformed vector x. */
147
148 /* INCX - INTEGER. */
149 /* On entry, INCX specifies the increment for the elements of */
150 /* X. INCX must not be zero. */
151 /* Unchanged on exit. */
152
153 /* Further Details */
154 /* =============== */
155
156 /* Level 2 Blas routine. */
157
158 /* -- Written on 22-October-1986. */
159 /* Jack Dongarra, Argonne National Lab. */
160 /* Jeremy Du Croz, Nag Central Office. */
161 /* Sven Hammarling, Nag Central Office. */
162 /* Richard Hanson, Sandia National Labs. */
163
164 /* ===================================================================== */
165
166 /* .. Parameters .. */
167 /* .. */
168 /* .. Local Scalars .. */
169 /* .. */
170 /* .. External Functions .. */
171 /* .. */
172 /* .. External Subroutines .. */
173 /* .. */
174 /* .. Intrinsic Functions .. */
175 /* .. */
176
177 /* Test the input parameters. */
178
179 /* Parameter adjustments */
180 a_dim1 = *lda;
181 a_offset = 1 + a_dim1;
182 a -= a_offset;
183 --x;
184
185 /* Function Body */
186 info = 0;
187 if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
188 ftnlen)1, (ftnlen)1)) {
189 info = 1;
190 } else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans,
191 "T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
192 ftnlen)1)) {
193 info = 2;
194 } else if (! lsame_(diag, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(diag,
195 "N", (ftnlen)1, (ftnlen)1)) {
196 info = 3;
197 } else if (*n < 0) {
198 info = 4;
199 } else if (*k < 0) {
200 info = 5;
201 } else if (*lda < *k + 1) {
202 info = 7;
203 } else if (*incx == 0) {
204 info = 9;
205 }
206 if (info != 0) {
207 xerbla_("DTBMV ", &info, (ftnlen)6);
208 return 0;
209 }
210
211 /* Quick return if possible. */
212
213 if (*n == 0) {
214 return 0;
215 }
216
217 nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
218
219 /* Set up the start point in X if the increment is not unity. This */
220 /* will be ( N - 1 )*INCX too small for descending loops. */
221
222 if (*incx <= 0) {
223 kx = 1 - (*n - 1) * *incx;
224 } else if (*incx != 1) {
225 kx = 1;
226 }
227
228 /* Start the operations. In this version the elements of A are */
229 /* accessed sequentially with one pass through A. */
230
231 if (lsame_(trans, "N", (ftnlen)1, (ftnlen)1)) {
232
233 /* Form x := A*x. */
234
235 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
236 kplus1 = *k + 1;
237 if (*incx == 1) {
238 i__1 = *n;
239 for (j = 1; j <= i__1; ++j) {
240 if (x[j] != 0.) {
241 temp = x[j];
242 l = kplus1 - j;
243 /* Computing MAX */
244 i__2 = 1, i__3 = j - *k;
245 i__4 = j - 1;
246 for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
247 x[i__] += temp * a[l + i__ + j * a_dim1];
248 /* L10: */
249 }
250 if (nounit) {
251 x[j] *= a[kplus1 + j * a_dim1];
252 }
253 }
254 /* L20: */
255 }
256 } else {
257 jx = kx;
258 i__1 = *n;
259 for (j = 1; j <= i__1; ++j) {
260 if (x[jx] != 0.) {
261 temp = x[jx];
262 ix = kx;
263 l = kplus1 - j;
264 /* Computing MAX */
265 i__4 = 1, i__2 = j - *k;
266 i__3 = j - 1;
267 for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
268 x[ix] += temp * a[l + i__ + j * a_dim1];
269 ix += *incx;
270 /* L30: */
271 }
272 if (nounit) {
273 x[jx] *= a[kplus1 + j * a_dim1];
274 }
275 }
276 jx += *incx;
277 if (j > *k) {
278 kx += *incx;
279 }
280 /* L40: */
281 }
282 }
283 } else {
284 if (*incx == 1) {
285 for (j = *n; j >= 1; --j) {
286 if (x[j] != 0.) {
287 temp = x[j];
288 l = 1 - j;
289 /* Computing MIN */
290 i__1 = *n, i__3 = j + *k;
291 i__4 = j + 1;
292 for (i__ = min(i__1,i__3); i__ >= i__4; --i__) {
293 x[i__] += temp * a[l + i__ + j * a_dim1];
294 /* L50: */
295 }
296 if (nounit) {
297 x[j] *= a[j * a_dim1 + 1];
298 }
299 }
300 /* L60: */
301 }
302 } else {
303 kx += (*n - 1) * *incx;
304 jx = kx;
305 for (j = *n; j >= 1; --j) {
306 if (x[jx] != 0.) {
307 temp = x[jx];
308 ix = kx;
309 l = 1 - j;
310 /* Computing MIN */
311 i__4 = *n, i__1 = j + *k;
312 i__3 = j + 1;
313 for (i__ = min(i__4,i__1); i__ >= i__3; --i__) {
314 x[ix] += temp * a[l + i__ + j * a_dim1];
315 ix -= *incx;
316 /* L70: */
317 }
318 if (nounit) {
319 x[jx] *= a[j * a_dim1 + 1];
320 }
321 }
322 jx -= *incx;
323 if (*n - j >= *k) {
324 kx -= *incx;
325 }
326 /* L80: */
327 }
328 }
329 }
330 } else {
331
332 /* Form x := A'*x. */
333
334 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
335 kplus1 = *k + 1;
336 if (*incx == 1) {
337 for (j = *n; j >= 1; --j) {
338 temp = x[j];
339 l = kplus1 - j;
340 if (nounit) {
341 temp *= a[kplus1 + j * a_dim1];
342 }
343 /* Computing MAX */
344 i__4 = 1, i__1 = j - *k;
345 i__3 = max(i__4,i__1);
346 for (i__ = j - 1; i__ >= i__3; --i__) {
347 temp += a[l + i__ + j * a_dim1] * x[i__];
348 /* L90: */
349 }
350 x[j] = temp;
351 /* L100: */
352 }
353 } else {
354 kx += (*n - 1) * *incx;
355 jx = kx;
356 for (j = *n; j >= 1; --j) {
357 temp = x[jx];
358 kx -= *incx;
359 ix = kx;
360 l = kplus1 - j;
361 if (nounit) {
362 temp *= a[kplus1 + j * a_dim1];
363 }
364 /* Computing MAX */
365 i__4 = 1, i__1 = j - *k;
366 i__3 = max(i__4,i__1);
367 for (i__ = j - 1; i__ >= i__3; --i__) {
368 temp += a[l + i__ + j * a_dim1] * x[ix];
369 ix -= *incx;
370 /* L110: */
371 }
372 x[jx] = temp;
373 jx -= *incx;
374 /* L120: */
375 }
376 }
377 } else {
378 if (*incx == 1) {
379 i__3 = *n;
380 for (j = 1; j <= i__3; ++j) {
381 temp = x[j];
382 l = 1 - j;
383 if (nounit) {
384 temp *= a[j * a_dim1 + 1];
385 }
386 /* Computing MIN */
387 i__1 = *n, i__2 = j + *k;
388 i__4 = min(i__1,i__2);
389 for (i__ = j + 1; i__ <= i__4; ++i__) {
390 temp += a[l + i__ + j * a_dim1] * x[i__];
391 /* L130: */
392 }
393 x[j] = temp;
394 /* L140: */
395 }
396 } else {
397 jx = kx;
398 i__3 = *n;
399 for (j = 1; j <= i__3; ++j) {
400 temp = x[jx];
401 kx += *incx;
402 ix = kx;
403 l = 1 - j;
404 if (nounit) {
405 temp *= a[j * a_dim1 + 1];
406 }
407 /* Computing MIN */
408 i__1 = *n, i__2 = j + *k;
409 i__4 = min(i__1,i__2);
410 for (i__ = j + 1; i__ <= i__4; ++i__) {
411 temp += a[l + i__ + j * a_dim1] * x[ix];
412 ix += *incx;
413 /* L150: */
414 }
415 x[jx] = temp;
416 jx += *incx;
417 /* L160: */
418 }
419 }
420 }
421 }
422
423 return 0;
424
425 /* End of DTBMV . */
426
427 } /* dtbmv_ */
428
429