dtbmv.c
Go to the documentation of this file.
00001 /* dtbmv.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Subroutine */ int dtbmv_(char *uplo, char *trans, char *diag, integer *n, 
00017         integer *k, doublereal *a, integer *lda, doublereal *x, integer *incx)
00018 {
00019     /* System generated locals */
00020     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
00021 
00022     /* Local variables */
00023     integer i__, j, l, ix, jx, kx, info;
00024     doublereal temp;
00025     extern logical lsame_(char *, char *);
00026     integer kplus1;
00027     extern /* Subroutine */ int xerbla_(char *, integer *);
00028     logical nounit;
00029 
00030 /*     .. Scalar Arguments .. */
00031 /*     .. */
00032 /*     .. Array Arguments .. */
00033 /*     .. */
00034 
00035 /*  Purpose */
00036 /*  ======= */
00037 
00038 /*  DTBMV  performs one of the matrix-vector operations */
00039 
00040 /*     x := A*x,   or   x := A'*x, */
00041 
00042 /*  where x is an n element vector and  A is an n by n unit, or non-unit, */
00043 /*  upper or lower triangular band matrix, with ( k + 1 ) diagonals. */
00044 
00045 /*  Arguments */
00046 /*  ========== */
00047 
00048 /*  UPLO   - CHARACTER*1. */
00049 /*           On entry, UPLO specifies whether the matrix is an upper or */
00050 /*           lower triangular matrix as follows: */
00051 
00052 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
00053 
00054 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
00055 
00056 /*           Unchanged on exit. */
00057 
00058 /*  TRANS  - CHARACTER*1. */
00059 /*           On entry, TRANS specifies the operation to be performed as */
00060 /*           follows: */
00061 
00062 /*              TRANS = 'N' or 'n'   x := A*x. */
00063 
00064 /*              TRANS = 'T' or 't'   x := A'*x. */
00065 
00066 /*              TRANS = 'C' or 'c'   x := A'*x. */
00067 
00068 /*           Unchanged on exit. */
00069 
00070 /*  DIAG   - CHARACTER*1. */
00071 /*           On entry, DIAG specifies whether or not A is unit */
00072 /*           triangular as follows: */
00073 
00074 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
00075 
00076 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
00077 /*                                  triangular. */
00078 
00079 /*           Unchanged on exit. */
00080 
00081 /*  N      - INTEGER. */
00082 /*           On entry, N specifies the order of the matrix A. */
00083 /*           N must be at least zero. */
00084 /*           Unchanged on exit. */
00085 
00086 /*  K      - INTEGER. */
00087 /*           On entry with UPLO = 'U' or 'u', K specifies the number of */
00088 /*           super-diagonals of the matrix A. */
00089 /*           On entry with UPLO = 'L' or 'l', K specifies the number of */
00090 /*           sub-diagonals of the matrix A. */
00091 /*           K must satisfy  0 .le. K. */
00092 /*           Unchanged on exit. */
00093 
00094 /*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ). */
00095 /*           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
00096 /*           by n part of the array A must contain the upper triangular */
00097 /*           band part of the matrix of coefficients, supplied column by */
00098 /*           column, with the leading diagonal of the matrix in row */
00099 /*           ( k + 1 ) of the array, the first super-diagonal starting at */
00100 /*           position 2 in row k, and so on. The top left k by k triangle */
00101 /*           of the array A is not referenced. */
00102 /*           The following program segment will transfer an upper */
00103 /*           triangular band matrix from conventional full matrix storage */
00104 /*           to band storage: */
00105 
00106 /*                 DO 20, J = 1, N */
00107 /*                    M = K + 1 - J */
00108 /*                    DO 10, I = MAX( 1, J - K ), J */
00109 /*                       A( M + I, J ) = matrix( I, J ) */
00110 /*              10    CONTINUE */
00111 /*              20 CONTINUE */
00112 
00113 /*           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
00114 /*           by n part of the array A must contain the lower triangular */
00115 /*           band part of the matrix of coefficients, supplied column by */
00116 /*           column, with the leading diagonal of the matrix in row 1 of */
00117 /*           the array, the first sub-diagonal starting at position 1 in */
00118 /*           row 2, and so on. The bottom right k by k triangle of the */
00119 /*           array A is not referenced. */
00120 /*           The following program segment will transfer a lower */
00121 /*           triangular band matrix from conventional full matrix storage */
00122 /*           to band storage: */
00123 
00124 /*                 DO 20, J = 1, N */
00125 /*                    M = 1 - J */
00126 /*                    DO 10, I = J, MIN( N, J + K ) */
00127 /*                       A( M + I, J ) = matrix( I, J ) */
00128 /*              10    CONTINUE */
00129 /*              20 CONTINUE */
00130 
00131 /*           Note that when DIAG = 'U' or 'u' the elements of the array A */
00132 /*           corresponding to the diagonal elements of the matrix are not */
00133 /*           referenced, but are assumed to be unity. */
00134 /*           Unchanged on exit. */
00135 
00136 /*  LDA    - INTEGER. */
00137 /*           On entry, LDA specifies the first dimension of A as declared */
00138 /*           in the calling (sub) program. LDA must be at least */
00139 /*           ( k + 1 ). */
00140 /*           Unchanged on exit. */
00141 
00142 /*  X      - DOUBLE PRECISION array of dimension at least */
00143 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
00144 /*           Before entry, the incremented array X must contain the n */
00145 /*           element vector x. On exit, X is overwritten with the */
00146 /*           tranformed vector x. */
00147 
00148 /*  INCX   - INTEGER. */
00149 /*           On entry, INCX specifies the increment for the elements of */
00150 /*           X. INCX must not be zero. */
00151 /*           Unchanged on exit. */
00152 
00153 
00154 /*  Level 2 Blas routine. */
00155 
00156 /*  -- Written on 22-October-1986. */
00157 /*     Jack Dongarra, Argonne National Lab. */
00158 /*     Jeremy Du Croz, Nag Central Office. */
00159 /*     Sven Hammarling, Nag Central Office. */
00160 /*     Richard Hanson, Sandia National Labs. */
00161 
00162 
00163 /*     .. Parameters .. */
00164 /*     .. */
00165 /*     .. Local Scalars .. */
00166 /*     .. */
00167 /*     .. External Functions .. */
00168 /*     .. */
00169 /*     .. External Subroutines .. */
00170 /*     .. */
00171 /*     .. Intrinsic Functions .. */
00172 /*     .. */
00173 
00174 /*     Test the input parameters. */
00175 
00176     /* Parameter adjustments */
00177     a_dim1 = *lda;
00178     a_offset = 1 + a_dim1;
00179     a -= a_offset;
00180     --x;
00181 
00182     /* Function Body */
00183     info = 0;
00184     if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) {
00185         info = 1;
00186     } else if (! lsame_(trans, "N") && ! lsame_(trans, 
00187             "T") && ! lsame_(trans, "C")) {
00188         info = 2;
00189     } else if (! lsame_(diag, "U") && ! lsame_(diag, 
00190             "N")) {
00191         info = 3;
00192     } else if (*n < 0) {
00193         info = 4;
00194     } else if (*k < 0) {
00195         info = 5;
00196     } else if (*lda < *k + 1) {
00197         info = 7;
00198     } else if (*incx == 0) {
00199         info = 9;
00200     }
00201     if (info != 0) {
00202         xerbla_("DTBMV ", &info);
00203         return 0;
00204     }
00205 
00206 /*     Quick return if possible. */
00207 
00208     if (*n == 0) {
00209         return 0;
00210     }
00211 
00212     nounit = lsame_(diag, "N");
00213 
00214 /*     Set up the start point in X if the increment is not unity. This */
00215 /*     will be  ( N - 1 )*INCX   too small for descending loops. */
00216 
00217     if (*incx <= 0) {
00218         kx = 1 - (*n - 1) * *incx;
00219     } else if (*incx != 1) {
00220         kx = 1;
00221     }
00222 
00223 /*     Start the operations. In this version the elements of A are */
00224 /*     accessed sequentially with one pass through A. */
00225 
00226     if (lsame_(trans, "N")) {
00227 
00228 /*         Form  x := A*x. */
00229 
00230         if (lsame_(uplo, "U")) {
00231             kplus1 = *k + 1;
00232             if (*incx == 1) {
00233                 i__1 = *n;
00234                 for (j = 1; j <= i__1; ++j) {
00235                     if (x[j] != 0.) {
00236                         temp = x[j];
00237                         l = kplus1 - j;
00238 /* Computing MAX */
00239                         i__2 = 1, i__3 = j - *k;
00240                         i__4 = j - 1;
00241                         for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
00242                             x[i__] += temp * a[l + i__ + j * a_dim1];
00243 /* L10: */
00244                         }
00245                         if (nounit) {
00246                             x[j] *= a[kplus1 + j * a_dim1];
00247                         }
00248                     }
00249 /* L20: */
00250                 }
00251             } else {
00252                 jx = kx;
00253                 i__1 = *n;
00254                 for (j = 1; j <= i__1; ++j) {
00255                     if (x[jx] != 0.) {
00256                         temp = x[jx];
00257                         ix = kx;
00258                         l = kplus1 - j;
00259 /* Computing MAX */
00260                         i__4 = 1, i__2 = j - *k;
00261                         i__3 = j - 1;
00262                         for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
00263                             x[ix] += temp * a[l + i__ + j * a_dim1];
00264                             ix += *incx;
00265 /* L30: */
00266                         }
00267                         if (nounit) {
00268                             x[jx] *= a[kplus1 + j * a_dim1];
00269                         }
00270                     }
00271                     jx += *incx;
00272                     if (j > *k) {
00273                         kx += *incx;
00274                     }
00275 /* L40: */
00276                 }
00277             }
00278         } else {
00279             if (*incx == 1) {
00280                 for (j = *n; j >= 1; --j) {
00281                     if (x[j] != 0.) {
00282                         temp = x[j];
00283                         l = 1 - j;
00284 /* Computing MIN */
00285                         i__1 = *n, i__3 = j + *k;
00286                         i__4 = j + 1;
00287                         for (i__ = min(i__1,i__3); i__ >= i__4; --i__) {
00288                             x[i__] += temp * a[l + i__ + j * a_dim1];
00289 /* L50: */
00290                         }
00291                         if (nounit) {
00292                             x[j] *= a[j * a_dim1 + 1];
00293                         }
00294                     }
00295 /* L60: */
00296                 }
00297             } else {
00298                 kx += (*n - 1) * *incx;
00299                 jx = kx;
00300                 for (j = *n; j >= 1; --j) {
00301                     if (x[jx] != 0.) {
00302                         temp = x[jx];
00303                         ix = kx;
00304                         l = 1 - j;
00305 /* Computing MIN */
00306                         i__4 = *n, i__1 = j + *k;
00307                         i__3 = j + 1;
00308                         for (i__ = min(i__4,i__1); i__ >= i__3; --i__) {
00309                             x[ix] += temp * a[l + i__ + j * a_dim1];
00310                             ix -= *incx;
00311 /* L70: */
00312                         }
00313                         if (nounit) {
00314                             x[jx] *= a[j * a_dim1 + 1];
00315                         }
00316                     }
00317                     jx -= *incx;
00318                     if (*n - j >= *k) {
00319                         kx -= *incx;
00320                     }
00321 /* L80: */
00322                 }
00323             }
00324         }
00325     } else {
00326 
00327 /*        Form  x := A'*x. */
00328 
00329         if (lsame_(uplo, "U")) {
00330             kplus1 = *k + 1;
00331             if (*incx == 1) {
00332                 for (j = *n; j >= 1; --j) {
00333                     temp = x[j];
00334                     l = kplus1 - j;
00335                     if (nounit) {
00336                         temp *= a[kplus1 + j * a_dim1];
00337                     }
00338 /* Computing MAX */
00339                     i__4 = 1, i__1 = j - *k;
00340                     i__3 = max(i__4,i__1);
00341                     for (i__ = j - 1; i__ >= i__3; --i__) {
00342                         temp += a[l + i__ + j * a_dim1] * x[i__];
00343 /* L90: */
00344                     }
00345                     x[j] = temp;
00346 /* L100: */
00347                 }
00348             } else {
00349                 kx += (*n - 1) * *incx;
00350                 jx = kx;
00351                 for (j = *n; j >= 1; --j) {
00352                     temp = x[jx];
00353                     kx -= *incx;
00354                     ix = kx;
00355                     l = kplus1 - j;
00356                     if (nounit) {
00357                         temp *= a[kplus1 + j * a_dim1];
00358                     }
00359 /* Computing MAX */
00360                     i__4 = 1, i__1 = j - *k;
00361                     i__3 = max(i__4,i__1);
00362                     for (i__ = j - 1; i__ >= i__3; --i__) {
00363                         temp += a[l + i__ + j * a_dim1] * x[ix];
00364                         ix -= *incx;
00365 /* L110: */
00366                     }
00367                     x[jx] = temp;
00368                     jx -= *incx;
00369 /* L120: */
00370                 }
00371             }
00372         } else {
00373             if (*incx == 1) {
00374                 i__3 = *n;
00375                 for (j = 1; j <= i__3; ++j) {
00376                     temp = x[j];
00377                     l = 1 - j;
00378                     if (nounit) {
00379                         temp *= a[j * a_dim1 + 1];
00380                     }
00381 /* Computing MIN */
00382                     i__1 = *n, i__2 = j + *k;
00383                     i__4 = min(i__1,i__2);
00384                     for (i__ = j + 1; i__ <= i__4; ++i__) {
00385                         temp += a[l + i__ + j * a_dim1] * x[i__];
00386 /* L130: */
00387                     }
00388                     x[j] = temp;
00389 /* L140: */
00390                 }
00391             } else {
00392                 jx = kx;
00393                 i__3 = *n;
00394                 for (j = 1; j <= i__3; ++j) {
00395                     temp = x[jx];
00396                     kx += *incx;
00397                     ix = kx;
00398                     l = 1 - j;
00399                     if (nounit) {
00400                         temp *= a[j * a_dim1 + 1];
00401                     }
00402 /* Computing MIN */
00403                     i__1 = *n, i__2 = j + *k;
00404                     i__4 = min(i__1,i__2);
00405                     for (i__ = j + 1; i__ <= i__4; ++i__) {
00406                         temp += a[l + i__ + j * a_dim1] * x[ix];
00407                         ix += *incx;
00408 /* L150: */
00409                     }
00410                     x[jx] = temp;
00411                     jx += *incx;
00412 /* L160: */
00413                 }
00414             }
00415         }
00416     }
00417 
00418     return 0;
00419 
00420 /*     End of DTBMV . */
00421 
00422 } /* dtbmv_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:49