Home | History | Annotate | Download | only in fortran
      1       SUBROUTINE ZHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
      2 *     .. Scalar Arguments ..
      3       DOUBLE COMPLEX ALPHA,BETA
      4       INTEGER INCX,INCY,K,LDA,N
      5       CHARACTER UPLO
      6 *     ..
      7 *     .. Array Arguments ..
      8       DOUBLE COMPLEX 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 *  Arguments
     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 *  Further Details
    124 *  ===============
    125 *
    126 *  Level 2 Blas routine.
    127 *
    128 *  -- Written on 22-October-1986.
    129 *     Jack Dongarra, Argonne National Lab.
    130 *     Jeremy Du Croz, Nag Central Office.
    131 *     Sven Hammarling, Nag Central Office.
    132 *     Richard Hanson, Sandia National Labs.
    133 *
    134 *  =====================================================================
    135 *
    136 *     .. Parameters ..
    137       DOUBLE COMPLEX ONE
    138       PARAMETER (ONE= (1.0D+0,0.0D+0))
    139       DOUBLE COMPLEX ZERO
    140       PARAMETER (ZERO= (0.0D+0,0.0D+0))
    141 *     ..
    142 *     .. Local Scalars ..
    143       DOUBLE COMPLEX TEMP1,TEMP2
    144       INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
    145 *     ..
    146 *     .. External Functions ..
    147       LOGICAL LSAME
    148       EXTERNAL LSAME
    149 *     ..
    150 *     .. External Subroutines ..
    151       EXTERNAL XERBLA
    152 *     ..
    153 *     .. Intrinsic Functions ..
    154       INTRINSIC DBLE,DCONJG,MAX,MIN
    155 *     ..
    156 *
    157 *     Test the input parameters.
    158 *
    159       INFO = 0
    160       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
    161           INFO = 1
    162       ELSE IF (N.LT.0) THEN
    163           INFO = 2
    164       ELSE IF (K.LT.0) THEN
    165           INFO = 3
    166       ELSE IF (LDA.LT. (K+1)) THEN
    167           INFO = 6
    168       ELSE IF (INCX.EQ.0) THEN
    169           INFO = 8
    170       ELSE IF (INCY.EQ.0) THEN
    171           INFO = 11
    172       END IF
    173       IF (INFO.NE.0) THEN
    174           CALL XERBLA('ZHBMV ',INFO)
    175           RETURN
    176       END IF
    177 *
    178 *     Quick return if possible.
    179 *
    180       IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
    181 *
    182 *     Set up the start points in  X  and  Y.
    183 *
    184       IF (INCX.GT.0) THEN
    185           KX = 1
    186       ELSE
    187           KX = 1 - (N-1)*INCX
    188       END IF
    189       IF (INCY.GT.0) THEN
    190           KY = 1
    191       ELSE
    192           KY = 1 - (N-1)*INCY
    193       END IF
    194 *
    195 *     Start the operations. In this version the elements of the array A
    196 *     are accessed sequentially with one pass through A.
    197 *
    198 *     First form  y := beta*y.
    199 *
    200       IF (BETA.NE.ONE) THEN
    201           IF (INCY.EQ.1) THEN
    202               IF (BETA.EQ.ZERO) THEN
    203                   DO 10 I = 1,N
    204                       Y(I) = ZERO
    205    10             CONTINUE
    206               ELSE
    207                   DO 20 I = 1,N
    208                       Y(I) = BETA*Y(I)
    209    20             CONTINUE
    210               END IF
    211           ELSE
    212               IY = KY
    213               IF (BETA.EQ.ZERO) THEN
    214                   DO 30 I = 1,N
    215                       Y(IY) = ZERO
    216                       IY = IY + INCY
    217    30             CONTINUE
    218               ELSE
    219                   DO 40 I = 1,N
    220                       Y(IY) = BETA*Y(IY)
    221                       IY = IY + INCY
    222    40             CONTINUE
    223               END IF
    224           END IF
    225       END IF
    226       IF (ALPHA.EQ.ZERO) RETURN
    227       IF (LSAME(UPLO,'U')) THEN
    228 *
    229 *        Form  y  when upper triangle of A is stored.
    230 *
    231           KPLUS1 = K + 1
    232           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
    233               DO 60 J = 1,N
    234                   TEMP1 = ALPHA*X(J)
    235                   TEMP2 = ZERO
    236                   L = KPLUS1 - J
    237                   DO 50 I = MAX(1,J-K),J - 1
    238                       Y(I) = Y(I) + TEMP1*A(L+I,J)
    239                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(I)
    240    50             CONTINUE
    241                   Y(J) = Y(J) + TEMP1*DBLE(A(KPLUS1,J)) + ALPHA*TEMP2
    242    60         CONTINUE
    243           ELSE
    244               JX = KX
    245               JY = KY
    246               DO 80 J = 1,N
    247                   TEMP1 = ALPHA*X(JX)
    248                   TEMP2 = ZERO
    249                   IX = KX
    250                   IY = KY
    251                   L = KPLUS1 - J
    252                   DO 70 I = MAX(1,J-K),J - 1
    253                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
    254                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(IX)
    255                       IX = IX + INCX
    256                       IY = IY + INCY
    257    70             CONTINUE
    258                   Y(JY) = Y(JY) + TEMP1*DBLE(A(KPLUS1,J)) + ALPHA*TEMP2
    259                   JX = JX + INCX
    260                   JY = JY + INCY
    261                   IF (J.GT.K) THEN
    262                       KX = KX + INCX
    263                       KY = KY + INCY
    264                   END IF
    265    80         CONTINUE
    266           END IF
    267       ELSE
    268 *
    269 *        Form  y  when lower triangle of A is stored.
    270 *
    271           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
    272               DO 100 J = 1,N
    273                   TEMP1 = ALPHA*X(J)
    274                   TEMP2 = ZERO
    275                   Y(J) = Y(J) + TEMP1*DBLE(A(1,J))
    276                   L = 1 - J
    277                   DO 90 I = J + 1,MIN(N,J+K)
    278                       Y(I) = Y(I) + TEMP1*A(L+I,J)
    279                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(I)
    280    90             CONTINUE
    281                   Y(J) = Y(J) + ALPHA*TEMP2
    282   100         CONTINUE
    283           ELSE
    284               JX = KX
    285               JY = KY
    286               DO 120 J = 1,N
    287                   TEMP1 = ALPHA*X(JX)
    288                   TEMP2 = ZERO
    289                   Y(JY) = Y(JY) + TEMP1*DBLE(A(1,J))
    290                   L = 1 - J
    291                   IX = JX
    292                   IY = JY
    293                   DO 110 I = J + 1,MIN(N,J+K)
    294                       IX = IX + INCX
    295                       IY = IY + INCY
    296                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
    297                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(IX)
    298   110             CONTINUE
    299                   Y(JY) = Y(JY) + ALPHA*TEMP2
    300                   JX = JX + INCX
    301                   JY = JY + INCY
    302   120         CONTINUE
    303           END IF
    304       END IF
    305 *
    306       RETURN
    307 *
    308 *     End of ZHBMV .
    309 *
    310       END
    311