LAPACK: BLAS/SRC/zhbmv.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
00002 *     .. Scalar Arguments ..
00003       DOUBLE COMPLEX ALPHA,BETA
00004       INTEGER INCX,INCY,K,LDA,N
00005       CHARACTER UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE COMPLEX A(LDA,*),X(*),Y(*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  ZHBMV  performs the matrix-vector  operation
00015 *
00016 *     y := alpha*A*x + beta*y,
00017 *
00018 *  where alpha and beta are scalars, x and y are n element vectors and
00019 *  A is an n by n hermitian band matrix, with k super-diagonals.
00020 *
00021 *  Arguments
00022 *  ==========
00023 *
00024 *  UPLO   - CHARACTER*1.
00025 *           On entry, UPLO specifies whether the upper or lower
00026 *           triangular part of the band matrix A is being supplied as
00027 *           follows:
00028 *
00029 *              UPLO = 'U' or 'u'   The upper triangular part of A is
00030 *                                  being supplied.
00031 *
00032 *              UPLO = 'L' or 'l'   The lower triangular part of A is
00033 *                                  being supplied.
00034 *
00035 *           Unchanged on exit.
00036 *
00037 *  N      - INTEGER.
00038 *           On entry, N specifies the order of the matrix A.
00039 *           N must be at least zero.
00040 *           Unchanged on exit.
00041 *
00042 *  K      - INTEGER.
00043 *           On entry, K specifies the number of super-diagonals of the
00044 *           matrix A. K must satisfy  0 .le. K.
00045 *           Unchanged on exit.
00046 *
00047 *  ALPHA  - COMPLEX*16      .
00048 *           On entry, ALPHA specifies the scalar alpha.
00049 *           Unchanged on exit.
00050 *
00051 *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
00052 *           Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
00053 *           by n part of the array A must contain the upper triangular
00054 *           band part of the hermitian matrix, supplied column by
00055 *           column, with the leading diagonal of the matrix in row
00056 *           ( k + 1 ) of the array, the first super-diagonal starting at
00057 *           position 2 in row k, and so on. The top left k by k triangle
00058 *           of the array A is not referenced.
00059 *           The following program segment will transfer the upper
00060 *           triangular part of a hermitian band matrix from conventional
00061 *           full matrix storage to band storage:
00062 *
00063 *                 DO 20, J = 1, N
00064 *                    M = K + 1 - J
00065 *                    DO 10, I = MAX( 1, J - K ), J
00066 *                       A( M + I, J ) = matrix( I, J )
00067 *              10    CONTINUE
00068 *              20 CONTINUE
00069 *
00070 *           Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
00071 *           by n part of the array A must contain the lower triangular
00072 *           band part of the hermitian matrix, supplied column by
00073 *           column, with the leading diagonal of the matrix in row 1 of
00074 *           the array, the first sub-diagonal starting at position 1 in
00075 *           row 2, and so on. The bottom right k by k triangle of the
00076 *           array A is not referenced.
00077 *           The following program segment will transfer the lower
00078 *           triangular part of a hermitian band matrix from conventional
00079 *           full matrix storage to band storage:
00080 *
00081 *                 DO 20, J = 1, N
00082 *                    M = 1 - J
00083 *                    DO 10, I = J, MIN( N, J + K )
00084 *                       A( M + I, J ) = matrix( I, J )
00085 *              10    CONTINUE
00086 *              20 CONTINUE
00087 *
00088 *           Note that the imaginary parts of the diagonal elements need
00089 *           not be set and are assumed to be zero.
00090 *           Unchanged on exit.
00091 *
00092 *  LDA    - INTEGER.
00093 *           On entry, LDA specifies the first dimension of A as declared
00094 *           in the calling (sub) program. LDA must be at least
00095 *           ( k + 1 ).
00096 *           Unchanged on exit.
00097 *
00098 *  X      - COMPLEX*16       array of DIMENSION at least
00099 *           ( 1 + ( n - 1 )*abs( INCX ) ).
00100 *           Before entry, the incremented array X must contain the
00101 *           vector x.
00102 *           Unchanged on exit.
00103 *
00104 *  INCX   - INTEGER.
00105 *           On entry, INCX specifies the increment for the elements of
00106 *           X. INCX must not be zero.
00107 *           Unchanged on exit.
00108 *
00109 *  BETA   - COMPLEX*16      .
00110 *           On entry, BETA specifies the scalar beta.
00111 *           Unchanged on exit.
00112 *
00113 *  Y      - COMPLEX*16       array of DIMENSION at least
00114 *           ( 1 + ( n - 1 )*abs( INCY ) ).
00115 *           Before entry, the incremented array Y must contain the
00116 *           vector y. On exit, Y is overwritten by the updated vector y.
00117 *
00118 *  INCY   - INTEGER.
00119 *           On entry, INCY specifies the increment for the elements of
00120 *           Y. INCY must not be zero.
00121 *           Unchanged on exit.
00122 *
00123 *  Further Details
00124 *  ===============
00125 *
00126 *  Level 2 Blas routine.
00127 *
00128 *  -- Written on 22-October-1986.
00129 *     Jack Dongarra, Argonne National Lab.
00130 *     Jeremy Du Croz, Nag Central Office.
00131 *     Sven Hammarling, Nag Central Office.
00132 *     Richard Hanson, Sandia National Labs.
00133 *
00134 *  =====================================================================
00135 *
00136 *     .. Parameters ..
00137       DOUBLE COMPLEX ONE
00138       PARAMETER (ONE= (1.0D+0,0.0D+0))
00139       DOUBLE COMPLEX ZERO
00140       PARAMETER (ZERO= (0.0D+0,0.0D+0))
00141 *     ..
00142 *     .. Local Scalars ..
00143       DOUBLE COMPLEX TEMP1,TEMP2
00144       INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
00145 *     ..
00146 *     .. External Functions ..
00147       LOGICAL LSAME
00148       EXTERNAL LSAME
00149 *     ..
00150 *     .. External Subroutines ..
00151       EXTERNAL XERBLA
00152 *     ..
00153 *     .. Intrinsic Functions ..
00154       INTRINSIC DBLE,DCONJG,MAX,MIN
00155 *     ..
00156 *
00157 *     Test the input parameters.
00158 *
00159       INFO = 0
00160       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
00161           INFO = 1
00162       ELSE IF (N.LT.0) THEN
00163           INFO = 2
00164       ELSE IF (K.LT.0) THEN
00165           INFO = 3
00166       ELSE IF (LDA.LT. (K+1)) THEN
00167           INFO = 6
00168       ELSE IF (INCX.EQ.0) THEN
00169           INFO = 8
00170       ELSE IF (INCY.EQ.0) THEN
00171           INFO = 11
00172       END IF
00173       IF (INFO.NE.0) THEN
00174           CALL XERBLA('ZHBMV ',INFO)
00175           RETURN
00176       END IF
00177 *
00178 *     Quick return if possible.
00179 *
00180       IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
00181 *
00182 *     Set up the start points in  X  and  Y.
00183 *
00184       IF (INCX.GT.0) THEN
00185           KX = 1
00186       ELSE
00187           KX = 1 - (N-1)*INCX
00188       END IF
00189       IF (INCY.GT.0) THEN
00190           KY = 1
00191       ELSE
00192           KY = 1 - (N-1)*INCY
00193       END IF
00194 *
00195 *     Start the operations. In this version the elements of the array A
00196 *     are accessed sequentially with one pass through A.
00197 *
00198 *     First form  y := beta*y.
00199 *
00200       IF (BETA.NE.ONE) THEN
00201           IF (INCY.EQ.1) THEN
00202               IF (BETA.EQ.ZERO) THEN
00203                   DO 10 I = 1,N
00204                       Y(I) = ZERO
00205    10             CONTINUE
00206               ELSE
00207                   DO 20 I = 1,N
00208                       Y(I) = BETA*Y(I)
00209    20             CONTINUE
00210               END IF
00211           ELSE
00212               IY = KY
00213               IF (BETA.EQ.ZERO) THEN
00214                   DO 30 I = 1,N
00215                       Y(IY) = ZERO
00216                       IY = IY + INCY
00217    30             CONTINUE
00218               ELSE
00219                   DO 40 I = 1,N
00220                       Y(IY) = BETA*Y(IY)
00221                       IY = IY + INCY
00222    40             CONTINUE
00223               END IF
00224           END IF
00225       END IF
00226       IF (ALPHA.EQ.ZERO) RETURN
00227       IF (LSAME(UPLO,'U')) THEN
00228 *
00229 *        Form  y  when upper triangle of A is stored.
00230 *
00231           KPLUS1 = K + 1
00232           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00233               DO 60 J = 1,N
00234                   TEMP1 = ALPHA*X(J)
00235                   TEMP2 = ZERO
00236                   L = KPLUS1 - J
00237                   DO 50 I = MAX(1,J-K),J - 1
00238                       Y(I) = Y(I) + TEMP1*A(L+I,J)
00239                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(I)
00240    50             CONTINUE
00241                   Y(J) = Y(J) + TEMP1*DBLE(A(KPLUS1,J)) + ALPHA*TEMP2
00242    60         CONTINUE
00243           ELSE
00244               JX = KX
00245               JY = KY
00246               DO 80 J = 1,N
00247                   TEMP1 = ALPHA*X(JX)
00248                   TEMP2 = ZERO
00249                   IX = KX
00250                   IY = KY
00251                   L = KPLUS1 - J
00252                   DO 70 I = MAX(1,J-K),J - 1
00253                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
00254                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(IX)
00255                       IX = IX + INCX
00256                       IY = IY + INCY
00257    70             CONTINUE
00258                   Y(JY) = Y(JY) + TEMP1*DBLE(A(KPLUS1,J)) + ALPHA*TEMP2
00259                   JX = JX + INCX
00260                   JY = JY + INCY
00261                   IF (J.GT.K) THEN
00262                       KX = KX + INCX
00263                       KY = KY + INCY
00264                   END IF
00265    80         CONTINUE
00266           END IF
00267       ELSE
00268 *
00269 *        Form  y  when lower triangle of A is stored.
00270 *
00271           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00272               DO 100 J = 1,N
00273                   TEMP1 = ALPHA*X(J)
00274                   TEMP2 = ZERO
00275                   Y(J) = Y(J) + TEMP1*DBLE(A(1,J))
00276                   L = 1 - J
00277                   DO 90 I = J + 1,MIN(N,J+K)
00278                       Y(I) = Y(I) + TEMP1*A(L+I,J)
00279                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(I)
00280    90             CONTINUE
00281                   Y(J) = Y(J) + ALPHA*TEMP2
00282   100         CONTINUE
00283           ELSE
00284               JX = KX
00285               JY = KY
00286               DO 120 J = 1,N
00287                   TEMP1 = ALPHA*X(JX)
00288                   TEMP2 = ZERO
00289                   Y(JY) = Y(JY) + TEMP1*DBLE(A(1,J))
00290                   L = 1 - J
00291                   IX = JX
00292                   IY = JY
00293                   DO 110 I = J + 1,MIN(N,J+K)
00294                       IX = IX + INCX
00295                       IY = IY + INCY
00296                       Y(IY) = Y(IY) + TEMP1*A(L+I,J)
00297                       TEMP2 = TEMP2 + DCONJG(A(L+I,J))*X(IX)
00298   110             CONTINUE
00299                   Y(JY) = Y(JY) + ALPHA*TEMP2
00300                   JX = JX + INCX
00301                   JY = JY + INCY
00302   120         CONTINUE
00303           END IF
00304       END IF
00305 *
00306       RETURN
00307 *
00308 *     End of ZHBMV .
00309 *
00310       END