1 SUBROUTINE DSPGVD( ITYPE, JOBZ, UPLO, N, AP, BP, W, Z, LDZ, WORK, 2 $ LWORK, 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, ITYPE, LDZ, LIWORK, LWORK, N 12* .. 13* .. Array Arguments .. 14 INTEGER IWORK( * ) 15 DOUBLE PRECISION AP( * ), BP( * ), W( * ), WORK( * ), 16 $ Z( LDZ, * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* DSPGVD computes all the eigenvalues, and optionally, the eigenvectors 23* of a real generalized symmetric-definite eigenproblem, of the form 24* A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and 25* B are assumed to be symmetric, stored in packed format, and B is also 26* positive definite. 27* If eigenvectors are desired, it uses a divide and conquer algorithm. 28* 29* The divide and conquer algorithm makes very mild assumptions about 30* floating point arithmetic. It will work on machines with a guard 31* digit in add/subtract, or on those binary machines without guard 32* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 33* Cray-2. It could conceivably fail on hexadecimal or decimal machines 34* without guard digits, but we know of none. 35* 36* Arguments 37* ========= 38* 39* ITYPE (input) INTEGER 40* Specifies the problem type to be solved: 41* = 1: A*x = (lambda)*B*x 42* = 2: A*B*x = (lambda)*x 43* = 3: B*A*x = (lambda)*x 44* 45* JOBZ (input) CHARACTER*1 46* = 'N': Compute eigenvalues only; 47* = 'V': Compute eigenvalues and eigenvectors. 48* 49* UPLO (input) CHARACTER*1 50* = 'U': Upper triangles of A and B are stored; 51* = 'L': Lower triangles of A and B are stored. 52* 53* N (input) INTEGER 54* The order of the matrices A and B. N >= 0. 55* 56* AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 57* On entry, the upper or lower triangle of the symmetric matrix 58* A, packed columnwise in a linear array. The j-th column of A 59* is stored in the array AP as follows: 60* if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 61* if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 62* 63* On exit, the contents of AP are destroyed. 64* 65* BP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2) 66* On entry, the upper or lower triangle of the symmetric matrix 67* B, packed columnwise in a linear array. The j-th column of B 68* is stored in the array BP as follows: 69* if UPLO = 'U', BP(i + (j-1)*j/2) = B(i,j) for 1<=i<=j; 70* if UPLO = 'L', BP(i + (j-1)*(2*n-j)/2) = B(i,j) for j<=i<=n. 71* 72* On exit, the triangular factor U or L from the Cholesky 73* factorization B = U**T*U or B = L*L**T, in the same storage 74* format as B. 75* 76* W (output) DOUBLE PRECISION array, dimension (N) 77* If INFO = 0, the eigenvalues in ascending order. 78* 79* Z (output) DOUBLE PRECISION array, dimension (LDZ, N) 80* If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of 81* eigenvectors. The eigenvectors are normalized as follows: 82* if ITYPE = 1 or 2, Z**T*B*Z = I; 83* if ITYPE = 3, Z**T*inv(B)*Z = I. 84* If JOBZ = 'N', then Z is not referenced. 85* 86* LDZ (input) INTEGER 87* The leading dimension of the array Z. LDZ >= 1, and if 88* JOBZ = 'V', LDZ >= max(1,N). 89* 90* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 91* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 92* 93* LWORK (input) INTEGER 94* The dimension of the array WORK. 95* If N <= 1, LWORK >= 1. 96* If JOBZ = 'N' and N > 1, LWORK >= 2*N. 97* If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2. 98* 99* If LWORK = -1, then a workspace query is assumed; the routine 100* only calculates the optimal size of the WORK array, returns 101* this value as the first entry of the WORK array, and no error 102* message related to LWORK is issued by XERBLA. 103* 104* IWORK (workspace/output) INTEGER array, dimension (LIWORK) 105* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 106* 107* LIWORK (input) INTEGER 108* The dimension of the array IWORK. 109* If JOBZ = 'N' or N <= 1, LIWORK >= 1. 110* If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 111* 112* If LIWORK = -1, then a workspace query is assumed; the 113* routine only calculates the optimal size of the IWORK array, 114* returns this value as the first entry of the IWORK array, and 115* no error message related to LIWORK is issued by XERBLA. 116* 117* INFO (output) INTEGER 118* = 0: successful exit 119* < 0: if INFO = -i, the i-th argument had an illegal value 120* > 0: DPPTRF or DSPEVD returned an error code: 121* <= N: if INFO = i, DSPEVD failed to converge; 122* i off-diagonal elements of an intermediate 123* tridiagonal form did not converge to zero; 124* > N: if INFO = N + i, for 1 <= i <= N, then the leading 125* minor of order i of B is not positive definite. 126* The factorization of B could not be completed and 127* no eigenvalues or eigenvectors were computed. 128* 129* Further Details 130* =============== 131* 132* Based on contributions by 133* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 134* 135* ===================================================================== 136* 137* .. Parameters .. 138 DOUBLE PRECISION TWO 139 PARAMETER ( TWO = 2.0D+0 ) 140* .. 141* .. Local Scalars .. 142 LOGICAL LQUERY, UPPER, WANTZ 143 CHARACTER TRANS 144 INTEGER J, LGN, LIWMIN, LWMIN, NEIG 145* .. 146* .. External Functions .. 147 LOGICAL LSAME 148 EXTERNAL LSAME 149* .. 150* .. External Subroutines .. 151 EXTERNAL DPPTRF, DSPEVD, DSPGST, DTPMV, DTPSV, XERBLA 152* .. 153* .. Intrinsic Functions .. 154 INTRINSIC DBLE, INT, LOG, MAX 155* .. 156* .. Executable Statements .. 157* 158* Test the input parameters. 159* 160 WANTZ = LSAME( JOBZ, 'V' ) 161 UPPER = LSAME( UPLO, 'U' ) 162 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 163* 164 INFO = 0 165 IF( N.LE.1 ) THEN 166 LGN = 0 167 LIWMIN = 1 168 LWMIN = 1 169 ELSE 170 LGN = INT( LOG( DBLE( N ) ) / LOG( TWO ) ) 171 IF( 2**LGN.LT.N ) 172 $ LGN = LGN + 1 173 IF( 2**LGN.LT.N ) 174 $ LGN = LGN + 1 175 IF( WANTZ ) THEN 176 LIWMIN = 3 + 5*N 177 LWMIN = 1 + 5*N + 2*N*LGN + 2*N**2 178 ELSE 179 LIWMIN = 1 180 LWMIN = 2*N 181 END IF 182 END IF 183* 184 IF( ITYPE.LT.0 .OR. ITYPE.GT.3 ) THEN 185 INFO = -1 186 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 187 INFO = -2 188 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 189 INFO = -3 190 ELSE IF( N.LT.0 ) THEN 191 INFO = -4 192 ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN 193 INFO = -9 194 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 195 INFO = -11 196 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 197 INFO = -13 198 END IF 199* 200 IF( INFO.EQ.0 ) THEN 201 WORK( 1 ) = LWMIN 202 IWORK( 1 ) = LIWMIN 203 END IF 204* 205 IF( INFO.NE.0 ) THEN 206 CALL XERBLA( 'DSPGVD', -INFO ) 207 RETURN 208 ELSE IF( LQUERY ) THEN 209 RETURN 210 END IF 211* 212* Quick return if possible 213* 214 IF( N.EQ.0 ) 215 $ RETURN 216* 217* Form a Cholesky factorization of BP. 218* 219 CALL DPPTRF( UPLO, N, BP, INFO ) 220 IF( INFO.NE.0 ) THEN 221 INFO = N + INFO 222 RETURN 223 END IF 224* 225* Transform problem to standard eigenvalue problem and solve. 226* 227 CALL DSPGST( ITYPE, UPLO, N, AP, BP, INFO ) 228 CALL DSPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK, IWORK, 229 $ LIWORK, INFO ) 230 LWMIN = MAX( DBLE( LWMIN ), DBLE( WORK( 1 ) ) ) 231 LIWMIN = MAX( DBLE( LIWMIN ), DBLE( IWORK( 1 ) ) ) 232* 233 IF( WANTZ ) THEN 234* 235* Backtransform eigenvectors to the original problem. 236* 237 NEIG = N 238 IF( INFO.GT.0 ) 239 $ NEIG = INFO - 1 240 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN 241* 242* For A*x=(lambda)*B*x and A*B*x=(lambda)*x; 243* backtransform eigenvectors: x = inv(L)'*y or inv(U)*y 244* 245 IF( UPPER ) THEN 246 TRANS = 'N' 247 ELSE 248 TRANS = 'T' 249 END IF 250* 251 DO 10 J = 1, NEIG 252 CALL DTPSV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ), 253 $ 1 ) 254 10 CONTINUE 255* 256 ELSE IF( ITYPE.EQ.3 ) THEN 257* 258* For B*A*x=(lambda)*x; 259* backtransform eigenvectors: x = L*y or U'*y 260* 261 IF( UPPER ) THEN 262 TRANS = 'T' 263 ELSE 264 TRANS = 'N' 265 END IF 266* 267 DO 20 J = 1, NEIG 268 CALL DTPMV( UPLO, TRANS, 'Non-unit', N, BP, Z( 1, J ), 269 $ 1 ) 270 20 CONTINUE 271 END IF 272 END IF 273* 274 WORK( 1 ) = LWMIN 275 IWORK( 1 ) = LIWMIN 276* 277 RETURN 278* 279* End of DSPGVD 280* 281 END 282