share/lapack/blas/fortran/dsbmv.f

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