1 SUBROUTINE SSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, 2 $ VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, 3 $ LWORK, IWORK, IFAIL, INFO ) 4* 5* -- LAPACK driver routine (version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* June 30, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER JOBZ, RANGE, UPLO 12 INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N 13 REAL ABSTOL, VL, VU 14* .. 15* .. Array Arguments .. 16 INTEGER IFAIL( * ), IWORK( * ) 17 REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * ), 18 $ Z( LDZ, * ) 19* .. 20* 21* Purpose 22* ======= 23* 24* SSYGVX computes selected eigenvalues, and optionally, eigenvectors 25* of a real generalized symmetric-definite eigenproblem, of the form 26* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A 27* and B are assumed to be symmetric and B is also positive definite. 28* Eigenvalues and eigenvectors can be selected by specifying either a 29* range of values or a range of indices for the desired eigenvalues. 30* 31* Arguments 32* ========= 33* 34* ITYPE (input) INTEGER 35* Specifies the problem type to be solved: 36* = 1: A*x = (lambda)*B*x 37* = 2: A*B*x = (lambda)*x 38* = 3: B*A*x = (lambda)*x 39* 40* JOBZ (input) CHARACTER*1 41* = 'N': Compute eigenvalues only; 42* = 'V': Compute eigenvalues and eigenvectors. 43* 44* RANGE (input) CHARACTER*1 45* = 'A': all eigenvalues will be found. 46* = 'V': all eigenvalues in the half-open interval (VL,VU] 47* will be found. 48* = 'I': the IL-th through IU-th eigenvalues will be found. 49* 50* UPLO (input) CHARACTER*1 51* = 'U': Upper triangle of A and B are stored; 52* = 'L': Lower triangle of A and B are stored. 53* 54* N (input) INTEGER 55* The order of the matrix pencil (A,B). N >= 0. 56* 57* A (input/output) REAL array, dimension (LDA, N) 58* On entry, the symmetric matrix A. If UPLO = 'U', the 59* leading N-by-N upper triangular part of A contains the 60* upper triangular part of the matrix A. If UPLO = 'L', 61* the leading N-by-N lower triangular part of A contains 62* the lower triangular part of the matrix A. 63* 64* On exit, the lower triangle (if UPLO='L') or the upper 65* triangle (if UPLO='U') of A, including the diagonal, is 66* destroyed. 67* 68* LDA (input) INTEGER 69* The leading dimension of the array A. LDA >= max(1,N). 70* 71* B (input/output) REAL array, dimension (LDA, N) 72* On entry, the symmetric matrix B. If UPLO = 'U', the 73* leading N-by-N upper triangular part of B contains the 74* upper triangular part of the matrix B. If UPLO = 'L', 75* the leading N-by-N lower triangular part of B contains 76* the lower triangular part of the matrix B. 77* 78* On exit, if INFO <= N, the part of B containing the matrix is 79* overwritten by the triangular factor U or L from the Cholesky 80* factorization B = U**T*U or B = L*L**T. 81* 82* LDB (input) INTEGER 83* The leading dimension of the array B. LDB >= max(1,N). 84* 85* VL (input) REAL 86* VU (input) REAL 87* If RANGE='V', the lower and upper bounds of the interval to 88* be searched for eigenvalues. VL < VU. 89* Not referenced if RANGE = 'A' or 'I'. 90* 91* IL (input) INTEGER 92* IU (input) INTEGER 93* If RANGE='I', the indices (in ascending order) of the 94* smallest and largest eigenvalues to be returned. 95* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 96* Not referenced if RANGE = 'A' or 'V'. 97* 98* ABSTOL (input) REAL 99* The absolute error tolerance for the eigenvalues. 100* An approximate eigenvalue is accepted as converged 101* when it is determined to lie in an interval [a,b] 102* of width less than or equal to 103* 104* ABSTOL + EPS * max( |a|,|b| ) , 105* 106* where EPS is the machine precision. If ABSTOL is less than 107* or equal to zero, then EPS*|T| will be used in its place, 108* where |T| is the 1-norm of the tridiagonal matrix obtained 109* by reducing A to tridiagonal form. 110* 111* Eigenvalues will be computed most accurately when ABSTOL is 112* set to twice the underflow threshold 2*DLAMCH('S'), not zero. 113* If this routine returns with INFO>0, indicating that some 114* eigenvectors did not converge, try setting ABSTOL to 115* 2*SLAMCH('S'). 116* 117* M (output) INTEGER 118* The total number of eigenvalues found. 0 <= M <= N. 119* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 120* 121* W (output) REAL array, dimension (N) 122* On normal exit, the first M elements contain the selected 123* eigenvalues in ascending order. 124* 125* Z (output) REAL array, dimension (LDZ, max(1,M)) 126* If JOBZ = 'N', then Z is not referenced. 127* If JOBZ = 'V', then if INFO = 0, the first M columns of Z 128* contain the orthonormal eigenvectors of the matrix A 129* corresponding to the selected eigenvalues, with the i-th 130* column of Z holding the eigenvector associated with W(i). 131* The eigenvectors are normalized as follows: 132* if ITYPE = 1 or 2, Z**T*B*Z = I; 133* if ITYPE = 3, Z**T*inv(B)*Z = I. 134* 135* If an eigenvector fails to converge, then that column of Z 136* contains the latest approximation to the eigenvector, and the 137* index of the eigenvector is returned in IFAIL. 138* Note: the user must ensure that at least max(1,M) columns are 139* supplied in the array Z; if RANGE = 'V', the exact value of M 140* is not known in advance and an upper bound must be used. 141* 142* LDZ (input) INTEGER 143* The leading dimension of the array Z. LDZ >= 1, and if 144* JOBZ = 'V', LDZ >= max(1,N). 145* 146* WORK (workspace/output) REAL array, dimension (LWORK) 147* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 148* 149* LWORK (input) INTEGER 150* The length of the array WORK. LWORK >= max(1,8*N). 151* For optimal efficiency, LWORK >= (NB+3)*N, 152* where NB is the blocksize for SSYTRD returned by ILAENV. 153* 154* If LWORK = -1, then a workspace query is assumed; the routine 155* only calculates the optimal size of the WORK array, returns 156* this value as the first entry of the WORK array, and no error 157* message related to LWORK is issued by XERBLA. 158* 159* IWORK (workspace) INTEGER array, dimension (5*N) 160* 161* IFAIL (output) INTEGER array, dimension (N) 162* If JOBZ = 'V', then if INFO = 0, the first M elements of 163* IFAIL are zero. If INFO > 0, then IFAIL contains the 164* indices of the eigenvectors that failed to converge. 165* If JOBZ = 'N', then IFAIL is not referenced. 166* 167* INFO (output) INTEGER 168* = 0: successful exit 169* < 0: if INFO = -i, the i-th argument had an illegal value 170* > 0: SPOTRF or SSYEVX returned an error code: 171* <= N: if INFO = i, SSYEVX failed to converge; 172* i eigenvectors failed to converge. Their indices 173* are stored in array IFAIL. 174* > N: if INFO = N + i, for 1 <= i <= N, then the leading 175* minor of order i of B is not positive definite. 176* The factorization of B could not be completed and 177* no eigenvalues or eigenvectors were computed. 178* 179* Further Details 180* =============== 181* 182* Based on contributions by 183* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 184* 185* ===================================================================== 186* 187* .. Parameters .. 188 REAL ONE 189 PARAMETER ( ONE = 1.0E+0 ) 190* .. 191* .. Local Scalars .. 192 LOGICAL ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ 193 CHARACTER TRANS 194 INTEGER LOPT, LWKOPT, NB 195* .. 196* .. External Functions .. 197 LOGICAL LSAME 198 INTEGER ILAENV 199 EXTERNAL ILAENV, LSAME 200* .. 201* .. External Subroutines .. 202 EXTERNAL SPOTRF, SSYEVX, SSYGST, STRMM, STRSM, XERBLA 203* .. 204* .. Intrinsic Functions .. 205 INTRINSIC MAX, MIN 206* .. 207* .. Executable Statements .. 208* 209* Test the input parameters. 210* 211 UPPER = LSAME( UPLO, 'U' ) 212 WANTZ = LSAME( JOBZ, 'V' ) 213 ALLEIG = LSAME( RANGE, 'A' ) 214 VALEIG = LSAME( RANGE, 'V' ) 215 INDEIG = LSAME( RANGE, 'I' ) 216 LQUERY = ( LWORK.EQ.-1 ) 217* 218 INFO = 0 219 IF( ITYPE.LT.0 .OR. ITYPE.GT.3 ) THEN 220 INFO = -1 221 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 222 INFO = -2 223 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 224 INFO = -3 225 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 226 INFO = -4 227 ELSE IF( N.LT.0 ) THEN 228 INFO = -5 229 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 230 INFO = -7 231 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 232 INFO = -9 233 ELSE IF( VALEIG .AND. N.GT.0 ) THEN 234 IF( VU.LE.VL ) INFO = -11 235 ELSE IF( INDEIG .AND. IL.LT.1 ) THEN 236 INFO = -12 237 ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN 238 INFO = -13 239 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 240 INFO = -18 241 ELSE IF( LWORK.LT.MAX( 1, 8*N ) .AND. .NOT.LQUERY ) THEN 242 INFO = -20 243 END IF 244* 245 IF( INFO.EQ.0 ) THEN 246 NB = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 ) 247 LWKOPT = ( NB+3 )*N 248 WORK( 1 ) = LWKOPT 249 END IF 250* 251 IF( INFO.NE.0 ) THEN 252 CALL XERBLA( 'SSYGVX', -INFO ) 253 RETURN 254 ELSE IF( LQUERY ) THEN 255 RETURN 256 END IF 257* 258* Quick return if possible 259* 260 M = 0 261 IF( N.EQ.0 ) THEN 262 WORK( 1 ) = 1 263 RETURN 264 END IF 265* 266* Form a Cholesky factorization of B. 267* 268 CALL SPOTRF( UPLO, N, B, LDB, INFO ) 269 IF( INFO.NE.0 ) THEN 270 INFO = N + INFO 271 RETURN 272 END IF 273* 274* Transform problem to standard eigenvalue problem and solve. 275* 276 CALL SSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 277 CALL SSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL, 278 $ M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO ) 279 LOPT = WORK( 1 ) 280* 281 IF( WANTZ ) THEN 282* 283* Backtransform eigenvectors to the original problem. 284* 285 IF( INFO.GT.0 ) 286 $ M = INFO - 1 287 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN 288* 289* For A*x=(lambda)*B*x and A*B*x=(lambda)*x; 290* backtransform eigenvectors: x = inv(L)'*y or inv(U)*y 291* 292 IF( UPPER ) THEN 293 TRANS = 'N' 294 ELSE 295 TRANS = 'T' 296 END IF 297* 298 CALL STRSM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B, 299 $ LDB, Z, LDZ ) 300* 301 ELSE IF( ITYPE.EQ.3 ) THEN 302* 303* For B*A*x=(lambda)*x; 304* backtransform eigenvectors: x = L*y or U'*y 305* 306 IF( UPPER ) THEN 307 TRANS = 'T' 308 ELSE 309 TRANS = 'N' 310 END IF 311* 312 CALL STRMM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B, 313 $ LDB, Z, LDZ ) 314 END IF 315 END IF 316* 317* Set WORK(1) to optimal workspace size. 318* 319 WORK( 1 ) = LWKOPT 320* 321 RETURN 322* 323* End of SSYGVX 324* 325 END 326