1 SUBROUTINE CHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, 2 $ LRWORK, IWORK, LIWORK, 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, LIWORK, LRWORK, LWORK, N 12* .. 13* .. Array Arguments .. 14 INTEGER IWORK( * ) 15 REAL RWORK( * ), W( * ) 16 COMPLEX A( LDA, * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* CHEEVD computes all eigenvalues and, optionally, eigenvectors of a 23* complex Hermitian matrix A. If eigenvectors are desired, it uses a 24* 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* A (input/output) COMPLEX array, dimension (LDA, N) 48* On entry, the Hermitian matrix A. If UPLO = 'U', the 49* leading N-by-N upper triangular part of A contains the 50* upper triangular part of the matrix A. If UPLO = 'L', 51* the leading N-by-N lower triangular part of A contains 52* the lower triangular part of the matrix A. 53* On exit, if JOBZ = 'V', then if INFO = 0, A contains the 54* orthonormal eigenvectors of the matrix A. 55* If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 56* or the upper triangle (if UPLO='U') of A, including the 57* diagonal, is destroyed. 58* 59* LDA (input) INTEGER 60* The leading dimension of the array A. LDA >= max(1,N). 61* 62* W (output) REAL array, dimension (N) 63* If INFO = 0, the eigenvalues in ascending order. 64* 65* WORK (workspace/output) COMPLEX array, dimension (LWORK) 66* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 67* 68* LWORK (input) INTEGER 69* The length of the array WORK. 70* If N <= 1, LWORK must be at least 1. 71* If JOBZ = 'N' and N > 1, LWORK must be at least N + 1. 72* If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2. 73* 74* If LWORK = -1, then a workspace query is assumed; the routine 75* only calculates the optimal size of the WORK array, returns 76* this value as the first entry of the WORK array, and no error 77* message related to LWORK is issued by XERBLA. 78* 79* RWORK (workspace/output) REAL array, 80* dimension (LRWORK) 81* On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 82* 83* LRWORK (input) INTEGER 84* The dimension of the array RWORK. 85* If N <= 1, LRWORK must be at least 1. 86* If JOBZ = 'N' and N > 1, LRWORK must be at least N. 87* If JOBZ = 'V' and N > 1, LRWORK must be at least 88* 1 + 5*N + 2*N**2. 89* 90* If LRWORK = -1, then a workspace query is assumed; the 91* routine only calculates the optimal size of the RWORK array, 92* returns this value as the first entry of the RWORK array, and 93* no error message related to LRWORK is issued by XERBLA. 94* 95* IWORK (workspace/output) INTEGER array, dimension (LIWORK) 96* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 97* 98* LIWORK (input) INTEGER 99* The dimension of the array IWORK. 100* If N <= 1, LIWORK must be at least 1. 101* If JOBZ = 'N' and N > 1, LIWORK must be at least 1. 102* If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 103* 104* If LIWORK = -1, then a workspace query is assumed; the 105* routine only calculates the optimal size of the IWORK array, 106* returns this value as the first entry of the IWORK array, and 107* no error message related to LIWORK is issued by XERBLA. 108* 109* INFO (output) INTEGER 110* = 0: successful exit 111* < 0: if INFO = -i, the i-th argument had an illegal value 112* > 0: if INFO = i, the algorithm failed to converge; i 113* off-diagonal elements of an intermediate tridiagonal 114* form did not converge to zero. 115* 116* Further Details 117* =============== 118* 119* Based on contributions by 120* Jeff Rutter, Computer Science Division, University of California 121* at Berkeley, USA 122* 123* ===================================================================== 124* 125* .. Parameters .. 126 REAL ZERO, ONE 127 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 128 COMPLEX CONE 129 PARAMETER ( CONE = ( 1.0E0, 0.0E0 ) ) 130* .. 131* .. Local Scalars .. 132 LOGICAL LOWER, LQUERY, WANTZ 133 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2, 134 $ INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK, 135 $ LLWRK2, LOPT, LROPT, LRWMIN, LWMIN 136 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 137 $ SMLNUM 138* .. 139* .. External Functions .. 140 LOGICAL LSAME 141 REAL CLANHE, SLAMCH 142 EXTERNAL LSAME, CLANHE, SLAMCH 143* .. 144* .. External Subroutines .. 145 EXTERNAL CHETRD, CLACPY, CLASCL, CSTEDC, CUNMTR, SSCAL, 146 $ SSTERF, XERBLA 147* .. 148* .. Intrinsic Functions .. 149 INTRINSIC INT, MAX, REAL, SQRT 150* .. 151* .. Executable Statements .. 152* 153* Test the input parameters. 154* 155 WANTZ = LSAME( JOBZ, 'V' ) 156 LOWER = LSAME( UPLO, 'L' ) 157 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 158* 159 INFO = 0 160 IF( N.LE.1 ) THEN 161 LWMIN = 1 162 LRWMIN = 1 163 LIWMIN = 1 164 LOPT = LWMIN 165 LROPT = LRWMIN 166 LIOPT = LIWMIN 167 ELSE 168 IF( WANTZ ) THEN 169 LWMIN = 2*N + N*N 170 LRWMIN = 1 + 5*N + 2*N**2 171 LIWMIN = 3 + 5*N 172 ELSE 173 LWMIN = N + 1 174 LRWMIN = N 175 LIWMIN = 1 176 END IF 177 LOPT = LWMIN 178 LROPT = LRWMIN 179 LIOPT = LIWMIN 180 END IF 181 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 182 INFO = -1 183 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 184 INFO = -2 185 ELSE IF( N.LT.0 ) THEN 186 INFO = -3 187 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 188 INFO = -5 189 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 190 INFO = -8 191 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 192 INFO = -10 193 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 194 INFO = -12 195 END IF 196* 197 IF( INFO.EQ.0 ) THEN 198 WORK( 1 ) = LOPT 199 RWORK( 1 ) = LROPT 200 IWORK( 1 ) = LIOPT 201 END IF 202* 203 IF( INFO.NE.0 ) THEN 204 CALL XERBLA( 'CHEEVD', -INFO ) 205 RETURN 206 ELSE IF( LQUERY ) THEN 207 RETURN 208 END IF 209* 210* Quick return if possible 211* 212 IF( N.EQ.0 ) 213 $ RETURN 214* 215 IF( N.EQ.1 ) THEN 216 W( 1 ) = A( 1, 1 ) 217 IF( WANTZ ) 218 $ A( 1, 1 ) = CONE 219 RETURN 220 END IF 221* 222* Get machine constants. 223* 224 SAFMIN = SLAMCH( 'Safe minimum' ) 225 EPS = SLAMCH( 'Precision' ) 226 SMLNUM = SAFMIN / EPS 227 BIGNUM = ONE / SMLNUM 228 RMIN = SQRT( SMLNUM ) 229 RMAX = SQRT( BIGNUM ) 230* 231* Scale matrix to allowable range, if necessary. 232* 233 ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK ) 234 ISCALE = 0 235 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 236 ISCALE = 1 237 SIGMA = RMIN / ANRM 238 ELSE IF( ANRM.GT.RMAX ) THEN 239 ISCALE = 1 240 SIGMA = RMAX / ANRM 241 END IF 242 IF( ISCALE.EQ.1 ) 243 $ CALL CLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 244* 245* Call CHETRD to reduce Hermitian matrix to tridiagonal form. 246* 247 INDE = 1 248 INDTAU = 1 249 INDWRK = INDTAU + N 250 INDRWK = INDE + N 251 INDWK2 = INDWRK + N*N 252 LLWORK = LWORK - INDWRK + 1 253 LLWRK2 = LWORK - INDWK2 + 1 254 LLRWK = LRWORK - INDRWK + 1 255 CALL CHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ), 256 $ WORK( INDWRK ), LLWORK, IINFO ) 257 LOPT = MAX( REAL( LOPT ), REAL( N )+REAL( WORK( INDWRK ) ) ) 258* 259* For eigenvalues only, call SSTERF. For eigenvectors, first call 260* CSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 261* tridiagonal matrix, then call CUNMTR to multiply it to the 262* Householder transformations represented as Householder vectors in 263* A. 264* 265 IF( .NOT.WANTZ ) THEN 266 CALL SSTERF( N, W, RWORK( INDE ), INFO ) 267 ELSE 268 CALL CSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N, 269 $ WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK, 270 $ IWORK, LIWORK, INFO ) 271 CALL CUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ), 272 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO ) 273 CALL CLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA ) 274 LOPT = MAX( LOPT, N+N**2+INT( WORK( INDWK2 ) ) ) 275 END IF 276* 277* If matrix was scaled, then rescale eigenvalues appropriately. 278* 279 IF( ISCALE.EQ.1 ) THEN 280 IF( INFO.EQ.0 ) THEN 281 IMAX = N 282 ELSE 283 IMAX = INFO - 1 284 END IF 285 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 286 END IF 287* 288 WORK( 1 ) = LOPT 289 RWORK( 1 ) = LROPT 290 IWORK( 1 ) = LIOPT 291* 292 RETURN 293* 294* End of CHEEVD 295* 296 END 297