1 SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, 2 $ INFO ) 3* 4* -- LAPACK driver routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER JOBZ, UPLO 11 INTEGER INFO, LDA, LWORK, N 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION RWORK( * ), W( * ) 15 COMPLEX*16 A( LDA, * ), WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* ZHEEV computes all eigenvalues and, optionally, eigenvectors of a 22* complex Hermitian matrix A. 23* 24* Arguments 25* ========= 26* 27* JOBZ (input) CHARACTER*1 28* = 'N': Compute eigenvalues only; 29* = 'V': Compute eigenvalues and eigenvectors. 30* 31* UPLO (input) CHARACTER*1 32* = 'U': Upper triangle of A is stored; 33* = 'L': Lower triangle of A is stored. 34* 35* N (input) INTEGER 36* The order of the matrix A. N >= 0. 37* 38* A (input/output) COMPLEX*16 array, dimension (LDA, N) 39* On entry, the Hermitian matrix A. If UPLO = 'U', the 40* leading N-by-N upper triangular part of A contains the 41* upper triangular part of the matrix A. If UPLO = 'L', 42* the leading N-by-N lower triangular part of A contains 43* the lower triangular part of the matrix A. 44* On exit, if JOBZ = 'V', then if INFO = 0, A contains the 45* orthonormal eigenvectors of the matrix A. 46* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 47* or the upper triangle (if UPLO='U') of A, including the 48* diagonal, is destroyed. 49* 50* LDA (input) INTEGER 51* The leading dimension of the array A. LDA >= max(1,N). 52* 53* W (output) DOUBLE PRECISION array, dimension (N) 54* If INFO = 0, the eigenvalues in ascending order. 55* 56* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) 57* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 58* 59* LWORK (input) INTEGER 60* The length of the array WORK. LWORK >= max(1,2*N-1). 61* For optimal efficiency, LWORK >= (NB+1)*N, 62* where NB is the blocksize for ZHETRD returned by ILAENV. 63* 64* If LWORK = -1, then a workspace query is assumed; the routine 65* only calculates the optimal size of the WORK array, returns 66* this value as the first entry of the WORK array, and no error 67* message related to LWORK is issued by XERBLA. 68* 69* RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2)) 70* 71* INFO (output) INTEGER 72* = 0: successful exit 73* < 0: if INFO = -i, the i-th argument had an illegal value 74* > 0: if INFO = i, the algorithm failed to converge; i 75* off-diagonal elements of an intermediate tridiagonal 76* form did not converge to zero. 77* 78* ===================================================================== 79* 80* .. Parameters .. 81 DOUBLE PRECISION ZERO, ONE 82 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 83 COMPLEX*16 CONE 84 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 85* .. 86* .. Local Scalars .. 87 LOGICAL LOWER, LQUERY, WANTZ 88 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE, 89 $ LLWORK, LOPT, LWKOPT, NB 90 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 91 $ SMLNUM 92* .. 93* .. External Functions .. 94 LOGICAL LSAME 95 INTEGER ILAENV 96 DOUBLE PRECISION DLAMCH, ZLANHE 97 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE 98* .. 99* .. External Subroutines .. 100 EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR, 101 $ ZUNGTR 102* .. 103* .. Intrinsic Functions .. 104 INTRINSIC MAX, SQRT 105* .. 106* .. Executable Statements .. 107* 108* Test the input parameters. 109* 110 WANTZ = LSAME( JOBZ, 'V' ) 111 LOWER = LSAME( UPLO, 'L' ) 112 LQUERY = ( LWORK.EQ.-1 ) 113* 114 INFO = 0 115 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 116 INFO = -1 117 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 118 INFO = -2 119 ELSE IF( N.LT.0 ) THEN 120 INFO = -3 121 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 122 INFO = -5 123 ELSE IF( LWORK.LT.MAX( 1, 2*N-1 ) .AND. .NOT.LQUERY ) THEN 124 INFO = -8 125 END IF 126* 127 IF( INFO.EQ.0 ) THEN 128 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) 129 LWKOPT = MAX( 1, ( NB+1 )*N ) 130 WORK( 1 ) = LWKOPT 131 END IF 132* 133 IF( INFO.NE.0 ) THEN 134 CALL XERBLA( 'ZHEEV ', -INFO ) 135 RETURN 136 ELSE IF( LQUERY ) THEN 137 RETURN 138 END IF 139* 140* Quick return if possible 141* 142 IF( N.EQ.0 ) THEN 143 WORK( 1 ) = 1 144 RETURN 145 END IF 146* 147 IF( N.EQ.1 ) THEN 148 W( 1 ) = A( 1, 1 ) 149 WORK( 1 ) = 3 150 IF( WANTZ ) 151 $ A( 1, 1 ) = CONE 152 RETURN 153 END IF 154* 155* Get machine constants. 156* 157 SAFMIN = DLAMCH( 'Safe minimum' ) 158 EPS = DLAMCH( 'Precision' ) 159 SMLNUM = SAFMIN / EPS 160 BIGNUM = ONE / SMLNUM 161 RMIN = SQRT( SMLNUM ) 162 RMAX = SQRT( BIGNUM ) 163* 164* Scale matrix to allowable range, if necessary. 165* 166 ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK ) 167 ISCALE = 0 168 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 169 ISCALE = 1 170 SIGMA = RMIN / ANRM 171 ELSE IF( ANRM.GT.RMAX ) THEN 172 ISCALE = 1 173 SIGMA = RMAX / ANRM 174 END IF 175 IF( ISCALE.EQ.1 ) 176 $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 177* 178* Call ZHETRD to reduce Hermitian matrix to tridiagonal form. 179* 180 INDE = 1 181 INDTAU = 1 182 INDWRK = INDTAU + N 183 LLWORK = LWORK - INDWRK + 1 184 CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ), 185 $ WORK( INDWRK ), LLWORK, IINFO ) 186 LOPT = N + WORK( INDWRK ) 187* 188* For eigenvalues only, call DSTERF. For eigenvectors, first call 189* ZUNGTR to generate the unitary matrix, then call ZSTEQR. 190* 191 IF( .NOT.WANTZ ) THEN 192 CALL DSTERF( N, W, RWORK( INDE ), INFO ) 193 ELSE 194 CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ), 195 $ LLWORK, IINFO ) 196 INDWRK = INDE + N 197 CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA, 198 $ RWORK( INDWRK ), INFO ) 199 END IF 200* 201* If matrix was scaled, then rescale eigenvalues appropriately. 202* 203 IF( ISCALE.EQ.1 ) THEN 204 IF( INFO.EQ.0 ) THEN 205 IMAX = N 206 ELSE 207 IMAX = INFO - 1 208 END IF 209 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 210 END IF 211* 212* Set WORK(1) to optimal complex workspace size. 213* 214 WORK( 1 ) = LWKOPT 215* 216 RETURN 217* 218* End of ZHEEV 219* 220 END 221