Annotation of rpl/lapack/lapack/zhbevd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZHBEVD( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK,
! 2: $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER JOBZ, UPLO
! 11: INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: INTEGER IWORK( * )
! 15: DOUBLE PRECISION RWORK( * ), W( * )
! 16: COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * ZHBEVD computes all the eigenvalues and, optionally, eigenvectors of
! 23: * a complex Hermitian band matrix A. If eigenvectors are desired, it
! 24: * uses a divide and conquer algorithm.
! 25: *
! 26: * The divide and conquer algorithm makes very mild assumptions about
! 27: * floating point arithmetic. It will work on machines with a guard
! 28: * digit in add/subtract, or on those binary machines without guard
! 29: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
! 30: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
! 31: * without guard digits, but we know of none.
! 32: *
! 33: * Arguments
! 34: * =========
! 35: *
! 36: * JOBZ (input) CHARACTER*1
! 37: * = 'N': Compute eigenvalues only;
! 38: * = 'V': Compute eigenvalues and eigenvectors.
! 39: *
! 40: * UPLO (input) CHARACTER*1
! 41: * = 'U': Upper triangle of A is stored;
! 42: * = 'L': Lower triangle of A is stored.
! 43: *
! 44: * N (input) INTEGER
! 45: * The order of the matrix A. N >= 0.
! 46: *
! 47: * KD (input) INTEGER
! 48: * The number of superdiagonals of the matrix A if UPLO = 'U',
! 49: * or the number of subdiagonals if UPLO = 'L'. KD >= 0.
! 50: *
! 51: * AB (input/output) COMPLEX*16 array, dimension (LDAB, N)
! 52: * On entry, the upper or lower triangle of the Hermitian band
! 53: * matrix A, stored in the first KD+1 rows of the array. The
! 54: * j-th column of A is stored in the j-th column of the array AB
! 55: * as follows:
! 56: * if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
! 57: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
! 58: *
! 59: * On exit, AB is overwritten by values generated during the
! 60: * reduction to tridiagonal form. If UPLO = 'U', the first
! 61: * superdiagonal and the diagonal of the tridiagonal matrix T
! 62: * are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
! 63: * the diagonal and first subdiagonal of T are returned in the
! 64: * first two rows of AB.
! 65: *
! 66: * LDAB (input) INTEGER
! 67: * The leading dimension of the array AB. LDAB >= KD + 1.
! 68: *
! 69: * W (output) DOUBLE PRECISION array, dimension (N)
! 70: * If INFO = 0, the eigenvalues in ascending order.
! 71: *
! 72: * Z (output) COMPLEX*16 array, dimension (LDZ, N)
! 73: * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
! 74: * eigenvectors of the matrix A, with the i-th column of Z
! 75: * holding the eigenvector associated with W(i).
! 76: * If JOBZ = 'N', then Z is not referenced.
! 77: *
! 78: * LDZ (input) INTEGER
! 79: * The leading dimension of the array Z. LDZ >= 1, and if
! 80: * JOBZ = 'V', LDZ >= max(1,N).
! 81: *
! 82: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 83: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 84: *
! 85: * LWORK (input) INTEGER
! 86: * The dimension of the array WORK.
! 87: * If N <= 1, LWORK must be at least 1.
! 88: * If JOBZ = 'N' and N > 1, LWORK must be at least N.
! 89: * If JOBZ = 'V' and N > 1, LWORK must be at least 2*N**2.
! 90: *
! 91: * If LWORK = -1, then a workspace query is assumed; the routine
! 92: * only calculates the optimal sizes of the WORK, RWORK and
! 93: * IWORK arrays, returns these values as the first entries of
! 94: * the WORK, RWORK and IWORK arrays, and no error message
! 95: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
! 96: *
! 97: * RWORK (workspace/output) DOUBLE PRECISION array,
! 98: * dimension (LRWORK)
! 99: * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
! 100: *
! 101: * LRWORK (input) INTEGER
! 102: * The dimension of array RWORK.
! 103: * If N <= 1, LRWORK must be at least 1.
! 104: * If JOBZ = 'N' and N > 1, LRWORK must be at least N.
! 105: * If JOBZ = 'V' and N > 1, LRWORK must be at least
! 106: * 1 + 5*N + 2*N**2.
! 107: *
! 108: * If LRWORK = -1, then a workspace query is assumed; the
! 109: * routine only calculates the optimal sizes of the WORK, RWORK
! 110: * and IWORK arrays, returns these values as the first entries
! 111: * of the WORK, RWORK and IWORK arrays, and no error message
! 112: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
! 113: *
! 114: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
! 115: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
! 116: *
! 117: * LIWORK (input) INTEGER
! 118: * The dimension of array IWORK.
! 119: * If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
! 120: * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N .
! 121: *
! 122: * If LIWORK = -1, then a workspace query is assumed; the
! 123: * routine only calculates the optimal sizes of the WORK, RWORK
! 124: * and IWORK arrays, returns these values as the first entries
! 125: * of the WORK, RWORK and IWORK arrays, and no error message
! 126: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
! 127: *
! 128: * INFO (output) INTEGER
! 129: * = 0: successful exit.
! 130: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 131: * > 0: if INFO = i, the algorithm failed to converge; i
! 132: * off-diagonal elements of an intermediate tridiagonal
! 133: * form did not converge to zero.
! 134: *
! 135: * =====================================================================
! 136: *
! 137: * .. Parameters ..
! 138: DOUBLE PRECISION ZERO, ONE
! 139: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 140: COMPLEX*16 CZERO, CONE
! 141: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
! 142: $ CONE = ( 1.0D0, 0.0D0 ) )
! 143: * ..
! 144: * .. Local Scalars ..
! 145: LOGICAL LOWER, LQUERY, WANTZ
! 146: INTEGER IINFO, IMAX, INDE, INDWK2, INDWRK, ISCALE,
! 147: $ LIWMIN, LLRWK, LLWK2, LRWMIN, LWMIN
! 148: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
! 149: $ SMLNUM
! 150: * ..
! 151: * .. External Functions ..
! 152: LOGICAL LSAME
! 153: DOUBLE PRECISION DLAMCH, ZLANHB
! 154: EXTERNAL LSAME, DLAMCH, ZLANHB
! 155: * ..
! 156: * .. External Subroutines ..
! 157: EXTERNAL DSCAL, DSTERF, XERBLA, ZGEMM, ZHBTRD, ZLACPY,
! 158: $ ZLASCL, ZSTEDC
! 159: * ..
! 160: * .. Intrinsic Functions ..
! 161: INTRINSIC SQRT
! 162: * ..
! 163: * .. Executable Statements ..
! 164: *
! 165: * Test the input parameters.
! 166: *
! 167: WANTZ = LSAME( JOBZ, 'V' )
! 168: LOWER = LSAME( UPLO, 'L' )
! 169: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 .OR. LRWORK.EQ.-1 )
! 170: *
! 171: INFO = 0
! 172: IF( N.LE.1 ) THEN
! 173: LWMIN = 1
! 174: LRWMIN = 1
! 175: LIWMIN = 1
! 176: ELSE
! 177: IF( WANTZ ) THEN
! 178: LWMIN = 2*N**2
! 179: LRWMIN = 1 + 5*N + 2*N**2
! 180: LIWMIN = 3 + 5*N
! 181: ELSE
! 182: LWMIN = N
! 183: LRWMIN = N
! 184: LIWMIN = 1
! 185: END IF
! 186: END IF
! 187: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
! 188: INFO = -1
! 189: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
! 190: INFO = -2
! 191: ELSE IF( N.LT.0 ) THEN
! 192: INFO = -3
! 193: ELSE IF( KD.LT.0 ) THEN
! 194: INFO = -4
! 195: ELSE IF( LDAB.LT.KD+1 ) THEN
! 196: INFO = -6
! 197: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
! 198: INFO = -9
! 199: END IF
! 200: *
! 201: IF( INFO.EQ.0 ) THEN
! 202: WORK( 1 ) = LWMIN
! 203: RWORK( 1 ) = LRWMIN
! 204: IWORK( 1 ) = LIWMIN
! 205: *
! 206: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
! 207: INFO = -11
! 208: ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
! 209: INFO = -13
! 210: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
! 211: INFO = -15
! 212: END IF
! 213: END IF
! 214: *
! 215: IF( INFO.NE.0 ) THEN
! 216: CALL XERBLA( 'ZHBEVD', -INFO )
! 217: RETURN
! 218: ELSE IF( LQUERY ) THEN
! 219: RETURN
! 220: END IF
! 221: *
! 222: * Quick return if possible
! 223: *
! 224: IF( N.EQ.0 )
! 225: $ RETURN
! 226: *
! 227: IF( N.EQ.1 ) THEN
! 228: W( 1 ) = AB( 1, 1 )
! 229: IF( WANTZ )
! 230: $ Z( 1, 1 ) = CONE
! 231: RETURN
! 232: END IF
! 233: *
! 234: * Get machine constants.
! 235: *
! 236: SAFMIN = DLAMCH( 'Safe minimum' )
! 237: EPS = DLAMCH( 'Precision' )
! 238: SMLNUM = SAFMIN / EPS
! 239: BIGNUM = ONE / SMLNUM
! 240: RMIN = SQRT( SMLNUM )
! 241: RMAX = SQRT( BIGNUM )
! 242: *
! 243: * Scale matrix to allowable range, if necessary.
! 244: *
! 245: ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
! 246: ISCALE = 0
! 247: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
! 248: ISCALE = 1
! 249: SIGMA = RMIN / ANRM
! 250: ELSE IF( ANRM.GT.RMAX ) THEN
! 251: ISCALE = 1
! 252: SIGMA = RMAX / ANRM
! 253: END IF
! 254: IF( ISCALE.EQ.1 ) THEN
! 255: IF( LOWER ) THEN
! 256: CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
! 257: ELSE
! 258: CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
! 259: END IF
! 260: END IF
! 261: *
! 262: * Call ZHBTRD to reduce Hermitian band matrix to tridiagonal form.
! 263: *
! 264: INDE = 1
! 265: INDWRK = INDE + N
! 266: INDWK2 = 1 + N*N
! 267: LLWK2 = LWORK - INDWK2 + 1
! 268: LLRWK = LRWORK - INDWRK + 1
! 269: CALL ZHBTRD( JOBZ, UPLO, N, KD, AB, LDAB, W, RWORK( INDE ), Z,
! 270: $ LDZ, WORK, IINFO )
! 271: *
! 272: * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEDC.
! 273: *
! 274: IF( .NOT.WANTZ ) THEN
! 275: CALL DSTERF( N, W, RWORK( INDE ), INFO )
! 276: ELSE
! 277: CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK, N, WORK( INDWK2 ),
! 278: $ LLWK2, RWORK( INDWRK ), LLRWK, IWORK, LIWORK,
! 279: $ INFO )
! 280: CALL ZGEMM( 'N', 'N', N, N, N, CONE, Z, LDZ, WORK, N, CZERO,
! 281: $ WORK( INDWK2 ), N )
! 282: CALL ZLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
! 283: END IF
! 284: *
! 285: * If matrix was scaled, then rescale eigenvalues appropriately.
! 286: *
! 287: IF( ISCALE.EQ.1 ) THEN
! 288: IF( INFO.EQ.0 ) THEN
! 289: IMAX = N
! 290: ELSE
! 291: IMAX = INFO - 1
! 292: END IF
! 293: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
! 294: END IF
! 295: *
! 296: WORK( 1 ) = LWMIN
! 297: RWORK( 1 ) = LRWMIN
! 298: IWORK( 1 ) = LIWMIN
! 299: RETURN
! 300: *
! 301: * End of ZHBEVD
! 302: *
! 303: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhbevd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack