share/lapack/blas/fortran/zhbmv.f

   1       SUBROUTINE ZHBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
   2      $                   BETA, Y, INCY )
   3 *     .. Scalar Arguments ..
   4       COMPLEX*16         ALPHA, BETA
   5       INTEGER            INCX, INCY, K, LDA, N
   6       CHARACTER*1        UPLO
   7 *     .. Array Arguments ..
   8       COMPLEX*16         A( LDA, * ), X( * ), Y( * )
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  ZHBMV  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 hermitian 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  - COMPLEX*16      .
  48 *           On entry, ALPHA specifies the scalar alpha.
  49 *           Unchanged on exit.
  50 *
  51 *  A      - COMPLEX*16       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 hermitian 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 hermitian 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 hermitian 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 hermitian 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 *           Note that the imaginary parts of the diagonal elements need
  89 *           not be set and are assumed to be zero.
  90 *           Unchanged on exit.
  91 *
  92 *  LDA    - INTEGER.
  93 *           On entry, LDA specifies the first dimension of A as declared
  94 *           in the calling (sub) program. LDA must be at least
  95 *           ( k + 1 ).
  96 *           Unchanged on exit.
  97 *
  98 *  X      - COMPLEX*16       array of DIMENSION at least
  99 *           ( 1 + ( n - 1 )*abs( INCX ) ).
 100 *           Before entry, the incremented array X must contain the
 101 *           vector x.
 102 *           Unchanged on exit.
 103 *
 104 *  INCX   - INTEGER.
 105 *           On entry, INCX specifies the increment for the elements of
 106 *           X. INCX must not be zero.
 107 *           Unchanged on exit.
 108 *
 109 *  BETA   - COMPLEX*16      .
 110 *           On entry, BETA specifies the scalar beta.
 111 *           Unchanged on exit.
 112 *
 113 *  Y      - COMPLEX*16       array of DIMENSION at least
 114 *           ( 1 + ( n - 1 )*abs( INCY ) ).
 115 *           Before entry, the incremented array Y must contain the
 116 *           vector y. On exit, Y is overwritten by the updated vector y.
 117 *
 118 *  INCY   - INTEGER.
 119 *           On entry, INCY specifies the increment for the elements of
 120 *           Y. INCY must not be zero.
 121 *           Unchanged on exit.
 122 *
 123 *
 124 *  Level 2 Blas routine.
 125 *
 126 *  -- Written on 22-October-1986.
 127 *     Jack Dongarra, Argonne National Lab.
 128 *     Jeremy Du Croz, Nag Central Office.
 129 *     Sven Hammarling, Nag Central Office.
 130 *     Richard Hanson, Sandia National Labs.
 131 *
 132 *
 133 *     .. Parameters ..
 134       COMPLEX*16         ONE
 135       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
 136       COMPLEX*16         ZERO
 137       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 138 *     .. Local Scalars ..
 139       COMPLEX*16         TEMP1, TEMP2
 140       INTEGER            I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L
 141 *     .. External Functions ..
 142       LOGICAL            LSAME
 143       EXTERNAL           LSAME
 144 *     .. External Subroutines ..
 145       EXTERNAL           XERBLA
 146 *     .. Intrinsic Functions ..
 147       INTRINSIC          DCONJG, MAX, MIN, DBLE
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Test the input parameters.
 152 *
 153       INFO = 0
 154       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
 155      $         .NOT.LSAME( UPLO, 'L' )      )THEN
 156          INFO = 1
 157       ELSE IF( N.LT.0 )THEN
 158          INFO = 2
 159       ELSE IF( K.LT.0 )THEN
 160          INFO = 3
 161       ELSE IF( LDA.LT.( K + 1 ) )THEN
 162          INFO = 6
 163       ELSE IF( INCX.EQ.0 )THEN
 164          INFO = 8
 165       ELSE IF( INCY.EQ.0 )THEN
 166          INFO = 11
 167       END IF
 168       IF( INFO.NE.0 )THEN
 169          CALL XERBLA( 'ZHBMV ', INFO )
 170          RETURN
 171       END IF
 172 *
 173 *     Quick return if possible.
 174 *
 175       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 176      $   RETURN
 177 *
 178 *     Set up the start points in  X  and  Y.
 179 *
 180       IF( INCX.GT.0 )THEN
 181          KX = 1
 182       ELSE
 183          KX = 1 - ( N - 1 )*INCX
 184       END IF
 185       IF( INCY.GT.0 )THEN
 186          KY = 1
 187       ELSE
 188          KY = 1 - ( N - 1 )*INCY
 189       END IF
 190 *
 191 *     Start the operations. In this version the elements of the array A
 192 *     are accessed sequentially with one pass through A.
 193 *
 194 *     First form  y := beta*y.
 195 *
 196       IF( BETA.NE.ONE )THEN
 197          IF( INCY.EQ.1 )THEN
 198             IF( BETA.EQ.ZERO )THEN
 199                DO 10, I = 1, N
 200                   Y( I ) = ZERO
 201    10          CONTINUE
 202             ELSE
 203                DO 20, I = 1, N
 204                   Y( I ) = BETA*Y( I )
 205    20          CONTINUE
 206             END IF
 207          ELSE
 208             IY = KY
 209             IF( BETA.EQ.ZERO )THEN
 210                DO 30, I = 1, N
 211                   Y( IY ) = ZERO
 212                   IY      = IY   + INCY
 213    30          CONTINUE
 214             ELSE
 215                DO 40, I = 1, N
 216                   Y( IY ) = BETA*Y( IY )
 217                   IY      = IY           + INCY
 218    40          CONTINUE
 219             END IF
 220          END IF
 221       END IF
 222       IF( ALPHA.EQ.ZERO )
 223      $   RETURN
 224       IF( LSAME( UPLO, 'U' ) )THEN
 225 *
 226 *        Form  y  when upper triangle of A is stored.
 227 *
 228          KPLUS1 = K + 1
 229          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
 230             DO 60, J = 1, N
 231                TEMP1 = ALPHA*X( J )
 232                TEMP2 = ZERO
 233                L     = KPLUS1 - J
 234                DO 50, I = MAX( 1, J - K ), J - 1
 235                   Y( I ) = Y( I ) + TEMP1*A( L + I, J )
 236                   TEMP2  = TEMP2  + DCONJG( A( L + I, J ) )*X( I )
 237    50          CONTINUE
 238                Y( J ) = Y( J ) + TEMP1*DBLE( A( KPLUS1, J ) )
 239      $                         + ALPHA*TEMP2
 240    60       CONTINUE
 241          ELSE
 242             JX = KX
 243             JY = KY
 244             DO 80, J = 1, N
 245                TEMP1 = ALPHA*X( JX )
 246                TEMP2 = ZERO
 247                IX    = KX
 248                IY    = KY
 249                L     = KPLUS1 - J
 250                DO 70, I = MAX( 1, J - K ), J - 1
 251                   Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
 252                   TEMP2   = TEMP2   + DCONJG( A( L + I, J ) )*X( IX )
 253                   IX      = IX      + INCX
 254                   IY      = IY      + INCY
 255    70          CONTINUE
 256                Y( JY ) = Y( JY ) + TEMP1*DBLE( A( KPLUS1, J ) )
 257      $                           + ALPHA*TEMP2
 258                JX      = JX      + INCX
 259                JY      = JY      + INCY
 260                IF( J.GT.K )THEN
 261                   KX = KX + INCX
 262                   KY = KY + INCY
 263                END IF
 264    80       CONTINUE
 265          END IF
 266       ELSE
 267 *
 268 *        Form  y  when lower triangle of A is stored.
 269 *
 270          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
 271             DO 100, J = 1, N
 272                TEMP1  = ALPHA*X( J )
 273                TEMP2  = ZERO
 274                Y( J ) = Y( J ) + TEMP1*DBLE( A( 1, J ) )
 275                L      = 1      - J
 276                DO 90, I = J + 1, MIN( N, J + K )
 277                   Y( I ) = Y( I ) + TEMP1*A( L + I, J )
 278                   TEMP2  = TEMP2  + DCONJG( A( L + I, J ) )*X( I )
 279    90          CONTINUE
 280                Y( J ) = Y( J ) + ALPHA*TEMP2
 281   100       CONTINUE
 282          ELSE
 283             JX = KX
 284             JY = KY
 285             DO 120, J = 1, N
 286                TEMP1   = ALPHA*X( JX )
 287                TEMP2   = ZERO
 288                Y( JY ) = Y( JY ) + TEMP1*DBLE( A( 1, J ) )
 289                L       = 1       - J
 290                IX      = JX
 291                IY      = JY
 292                DO 110, I = J + 1, MIN( N, J + K )
 293                   IX      = IX      + INCX
 294                   IY      = IY      + INCY
 295                   Y( IY ) = Y( IY ) + TEMP1*A( L + I, J )
 296                   TEMP2   = TEMP2   + DCONJG( A( L + I, J ) )*X( IX )
 297   110          CONTINUE
 298                Y( JY ) = Y( JY ) + ALPHA*TEMP2
 299                JX      = JX      + INCX
 300                JY      = JY      + INCY
 301   120       CONTINUE
 302          END IF
 303       END IF
 304 *
 305       RETURN
 306 *
 307 *     End of ZHBMV .
 308 *
 309       END