share/lapack/blas/fortran/dgbmv.f

   1       SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
   2      $                   BETA, Y, INCY )
   3 *     .. Scalar Arguments ..
   4       DOUBLE PRECISION   ALPHA, BETA
   5       INTEGER            INCX, INCY, KL, KU, LDA, M, N
   6       CHARACTER*1        TRANS
   7 *     .. Array Arguments ..
   8       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  DGBMV  performs one of the matrix-vector operations
  15 *
  16 *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
  17 *
  18 *  where alpha and beta are scalars, x and y are vectors and A is an
  19 *  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
  20 *
  21 *  Parameters
  22 *  ==========
  23 *
  24 *  TRANS  - CHARACTER*1.
  25 *           On entry, TRANS specifies the operation to be performed as
  26 *           follows:
  27 *
  28 *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  29 *
  30 *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  31 *
  32 *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
  33 *
  34 *           Unchanged on exit.
  35 *
  36 *  M      - INTEGER.
  37 *           On entry, M specifies the number of rows of the matrix A.
  38 *           M must be at least zero.
  39 *           Unchanged on exit.
  40 *
  41 *  N      - INTEGER.
  42 *           On entry, N specifies the number of columns of the matrix A.
  43 *           N must be at least zero.
  44 *           Unchanged on exit.
  45 *
  46 *  KL     - INTEGER.
  47 *           On entry, KL specifies the number of sub-diagonals of the
  48 *           matrix A. KL must satisfy  0 .le. KL.
  49 *           Unchanged on exit.
  50 *
  51 *  KU     - INTEGER.
  52 *           On entry, KU specifies the number of super-diagonals of the
  53 *           matrix A. KU must satisfy  0 .le. KU.
  54 *           Unchanged on exit.
  55 *
  56 *  ALPHA  - DOUBLE PRECISION.
  57 *           On entry, ALPHA specifies the scalar alpha.
  58 *           Unchanged on exit.
  59 *
  60 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  61 *           Before entry, the leading ( kl + ku + 1 ) by n part of the
  62 *           array A must contain the matrix of coefficients, supplied
  63 *           column by column, with the leading diagonal of the matrix in
  64 *           row ( ku + 1 ) of the array, the first super-diagonal
  65 *           starting at position 2 in row ku, the first sub-diagonal
  66 *           starting at position 1 in row ( ku + 2 ), and so on.
  67 *           Elements in the array A that do not correspond to elements
  68 *           in the band matrix (such as the top left ku by ku triangle)
  69 *           are not referenced.
  70 *           The following program segment will transfer a band matrix
  71 *           from conventional full matrix storage to band storage:
  72 *
  73 *                 DO 20, J = 1, N
  74 *                    K = KU + 1 - J
  75 *                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
  76 *                       A( K + I, J ) = matrix( I, J )
  77 *              10    CONTINUE
  78 *              20 CONTINUE
  79 *
  80 *           Unchanged on exit.
  81 *
  82 *  LDA    - INTEGER.
  83 *           On entry, LDA specifies the first dimension of A as declared
  84 *           in the calling (sub) program. LDA must be at least
  85 *           ( kl + ku + 1 ).
  86 *           Unchanged on exit.
  87 *
  88 *  X      - DOUBLE PRECISION array of DIMENSION at least
  89 *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  90 *           and at least
  91 *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  92 *           Before entry, the incremented array X must contain the
  93 *           vector x.
  94 *           Unchanged on exit.
  95 *
  96 *  INCX   - INTEGER.
  97 *           On entry, INCX specifies the increment for the elements of
  98 *           X. INCX must not be zero.
  99 *           Unchanged on exit.
 100 *
 101 *  BETA   - DOUBLE PRECISION.
 102 *           On entry, BETA specifies the scalar beta. When BETA is
 103 *           supplied as zero then Y need not be set on input.
 104 *           Unchanged on exit.
 105 *
 106 *  Y      - DOUBLE PRECISION array of DIMENSION at least
 107 *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
 108 *           and at least
 109 *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
 110 *           Before entry, the incremented array Y must contain the
 111 *           vector y. On exit, Y is overwritten by the updated vector y.
 112 *
 113 *  INCY   - INTEGER.
 114 *           On entry, INCY specifies the increment for the elements of
 115 *           Y. INCY must not be zero.
 116 *           Unchanged on exit.
 117 *
 118 *
 119 *  Level 2 Blas routine.
 120 *
 121 *  -- Written on 22-October-1986.
 122 *     Jack Dongarra, Argonne National Lab.
 123 *     Jeremy Du Croz, Nag Central Office.
 124 *     Sven Hammarling, Nag Central Office.
 125 *     Richard Hanson, Sandia National Labs.
 126 *
 127 *     .. Parameters ..
 128       DOUBLE PRECISION   ONE         , ZERO
 129       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 130 *     .. Local Scalars ..
 131       DOUBLE PRECISION   TEMP
 132       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY,
 133      $                   LENX, LENY
 134 *     .. External Functions ..
 135       LOGICAL            LSAME
 136       EXTERNAL           LSAME
 137 *     .. External Subroutines ..
 138       EXTERNAL           XERBLA
 139 *     .. Intrinsic Functions ..
 140       INTRINSIC          MAX, MIN
 141 *     ..
 142 *     .. Executable Statements ..
 143 *
 144 *     Test the input parameters.
 145 *
 146       INFO = 0
 147       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
 148      $         .NOT.LSAME( TRANS, 'T' ).AND.
 149      $         .NOT.LSAME( TRANS, 'C' )      )THEN
 150          INFO = 1
 151       ELSE IF( M.LT.0 )THEN
 152          INFO = 2
 153       ELSE IF( N.LT.0 )THEN
 154          INFO = 3
 155       ELSE IF( KL.LT.0 )THEN
 156          INFO = 4
 157       ELSE IF( KU.LT.0 )THEN
 158          INFO = 5
 159       ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN
 160          INFO = 8
 161       ELSE IF( INCX.EQ.0 )THEN
 162          INFO = 10
 163       ELSE IF( INCY.EQ.0 )THEN
 164          INFO = 13
 165       END IF
 166       IF( INFO.NE.0 )THEN
 167          CALL XERBLA( 'DGBMV ', INFO )
 168          RETURN
 169       END IF
 170 *
 171 *     Quick return if possible.
 172 *
 173       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
 174      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 175      $   RETURN
 176 *
 177 *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
 178 *     up the start points in  X  and  Y.
 179 *
 180       IF( LSAME( TRANS, 'N' ) )THEN
 181          LENX = N
 182          LENY = M
 183       ELSE
 184          LENX = M
 185          LENY = N
 186       END IF
 187       IF( INCX.GT.0 )THEN
 188          KX = 1
 189       ELSE
 190          KX = 1 - ( LENX - 1 )*INCX
 191       END IF
 192       IF( INCY.GT.0 )THEN
 193          KY = 1
 194       ELSE
 195          KY = 1 - ( LENY - 1 )*INCY
 196       END IF
 197 *
 198 *     Start the operations. In this version the elements of A are
 199 *     accessed sequentially with one pass through the band part of A.
 200 *
 201 *     First form  y := beta*y.
 202 *
 203       IF( BETA.NE.ONE )THEN
 204          IF( INCY.EQ.1 )THEN
 205             IF( BETA.EQ.ZERO )THEN
 206                DO 10, I = 1, LENY
 207                   Y( I ) = ZERO
 208    10          CONTINUE
 209             ELSE
 210                DO 20, I = 1, LENY
 211                   Y( I ) = BETA*Y( I )
 212    20          CONTINUE
 213             END IF
 214          ELSE
 215             IY = KY
 216             IF( BETA.EQ.ZERO )THEN
 217                DO 30, I = 1, LENY
 218                   Y( IY ) = ZERO
 219                   IY      = IY   + INCY
 220    30          CONTINUE
 221             ELSE
 222                DO 40, I = 1, LENY
 223                   Y( IY ) = BETA*Y( IY )
 224                   IY      = IY           + INCY
 225    40          CONTINUE
 226             END IF
 227          END IF
 228       END IF
 229       IF( ALPHA.EQ.ZERO )
 230      $   RETURN
 231       KUP1 = KU + 1
 232       IF( LSAME( TRANS, 'N' ) )THEN
 233 *
 234 *        Form  y := alpha*A*x + y.
 235 *
 236          JX = KX
 237          IF( INCY.EQ.1 )THEN
 238             DO 60, J = 1, N
 239                IF( X( JX ).NE.ZERO )THEN
 240                   TEMP = ALPHA*X( JX )
 241                   K    = KUP1 - J
 242                   DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL )
 243                      Y( I ) = Y( I ) + TEMP*A( K + I, J )
 244    50             CONTINUE
 245                END IF
 246                JX = JX + INCX
 247    60       CONTINUE
 248          ELSE
 249             DO 80, J = 1, N
 250                IF( X( JX ).NE.ZERO )THEN
 251                   TEMP = ALPHA*X( JX )
 252                   IY   = KY
 253                   K    = KUP1 - J
 254                   DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL )
 255                      Y( IY ) = Y( IY ) + TEMP*A( K + I, J )
 256                      IY      = IY      + INCY
 257    70             CONTINUE
 258                END IF
 259                JX = JX + INCX
 260                IF( J.GT.KU )
 261      $            KY = KY + INCY
 262    80       CONTINUE
 263          END IF
 264       ELSE
 265 *
 266 *        Form  y := alpha*A'*x + y.
 267 *
 268          JY = KY
 269          IF( INCX.EQ.1 )THEN
 270             DO 100, J = 1, N
 271                TEMP = ZERO
 272                K    = KUP1 - J
 273                DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL )
 274                   TEMP = TEMP + A( K + I, J )*X( I )
 275    90          CONTINUE
 276                Y( JY ) = Y( JY ) + ALPHA*TEMP
 277                JY      = JY      + INCY
 278   100       CONTINUE
 279          ELSE
 280             DO 120, J = 1, N
 281                TEMP = ZERO
 282                IX   = KX
 283                K    = KUP1 - J
 284                DO 110, I = MAX( 1, J - KU ), MIN( M, J + KL )
 285                   TEMP = TEMP + A( K + I, J )*X( IX )
 286                   IX   = IX   + INCX
 287   110          CONTINUE
 288                Y( JY ) = Y( JY ) + ALPHA*TEMP
 289                JY      = JY      + INCY
 290                IF( J.GT.KU )
 291      $            KX = KX + INCX
 292   120       CONTINUE
 293          END IF
 294       END IF
 295 *
 296       RETURN
 297 *
 298 *     End of DGBMV .
 299 *
 300       END