1*> \brief \b DSYGVD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DSYGVD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygvd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygvd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygvd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, 22* LWORK, IWORK, LIWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER JOBZ, UPLO 26* INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IWORK( * ) 30* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> DSYGVD computes all the eigenvalues, and optionally, the eigenvectors 40*> of a real generalized symmetric-definite eigenproblem, of the form 41*> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and 42*> B are assumed to be symmetric and B is also positive definite. 43*> If eigenvectors are desired, it uses a divide and conquer algorithm. 44*> 45*> The divide and conquer algorithm makes very mild assumptions about 46*> floating point arithmetic. It will work on machines with a guard 47*> digit in add/subtract, or on those binary machines without guard 48*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 49*> Cray-2. It could conceivably fail on hexadecimal or decimal machines 50*> without guard digits, but we know of none. 51*> \endverbatim 52* 53* Arguments: 54* ========== 55* 56*> \param[in] ITYPE 57*> \verbatim 58*> ITYPE is INTEGER 59*> Specifies the problem type to be solved: 60*> = 1: A*x = (lambda)*B*x 61*> = 2: A*B*x = (lambda)*x 62*> = 3: B*A*x = (lambda)*x 63*> \endverbatim 64*> 65*> \param[in] JOBZ 66*> \verbatim 67*> JOBZ is CHARACTER*1 68*> = 'N': Compute eigenvalues only; 69*> = 'V': Compute eigenvalues and eigenvectors. 70*> \endverbatim 71*> 72*> \param[in] UPLO 73*> \verbatim 74*> UPLO is CHARACTER*1 75*> = 'U': Upper triangles of A and B are stored; 76*> = 'L': Lower triangles of A and B are stored. 77*> \endverbatim 78*> 79*> \param[in] N 80*> \verbatim 81*> N is INTEGER 82*> The order of the matrices A and B. N >= 0. 83*> \endverbatim 84*> 85*> \param[in,out] A 86*> \verbatim 87*> A is DOUBLE PRECISION array, dimension (LDA, N) 88*> On entry, the symmetric matrix A. If UPLO = 'U', the 89*> leading N-by-N upper triangular part of A contains the 90*> upper triangular part of the matrix A. If UPLO = 'L', 91*> the leading N-by-N lower triangular part of A contains 92*> the lower triangular part of the matrix A. 93*> 94*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the 95*> matrix Z of eigenvectors. The eigenvectors are normalized 96*> as follows: 97*> if ITYPE = 1 or 2, Z**T*B*Z = I; 98*> if ITYPE = 3, Z**T*inv(B)*Z = I. 99*> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') 100*> or the lower triangle (if UPLO='L') of A, including the 101*> diagonal, is destroyed. 102*> \endverbatim 103*> 104*> \param[in] LDA 105*> \verbatim 106*> LDA is INTEGER 107*> The leading dimension of the array A. LDA >= max(1,N). 108*> \endverbatim 109*> 110*> \param[in,out] B 111*> \verbatim 112*> B is DOUBLE PRECISION array, dimension (LDB, N) 113*> On entry, the symmetric matrix B. If UPLO = 'U', the 114*> leading N-by-N upper triangular part of B contains the 115*> upper triangular part of the matrix B. If UPLO = 'L', 116*> the leading N-by-N lower triangular part of B contains 117*> the lower triangular part of the matrix B. 118*> 119*> On exit, if INFO <= N, the part of B containing the matrix is 120*> overwritten by the triangular factor U or L from the Cholesky 121*> factorization B = U**T*U or B = L*L**T. 122*> \endverbatim 123*> 124*> \param[in] LDB 125*> \verbatim 126*> LDB is INTEGER 127*> The leading dimension of the array B. LDB >= max(1,N). 128*> \endverbatim 129*> 130*> \param[out] W 131*> \verbatim 132*> W is DOUBLE PRECISION array, dimension (N) 133*> If INFO = 0, the eigenvalues in ascending order. 134*> \endverbatim 135*> 136*> \param[out] WORK 137*> \verbatim 138*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 139*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 140*> \endverbatim 141*> 142*> \param[in] LWORK 143*> \verbatim 144*> LWORK is INTEGER 145*> The dimension of the array WORK. 146*> If N <= 1, LWORK >= 1. 147*> If JOBZ = 'N' and N > 1, LWORK >= 2*N+1. 148*> If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2. 149*> 150*> If LWORK = -1, then a workspace query is assumed; the routine 151*> only calculates the optimal sizes of the WORK and IWORK 152*> arrays, returns these values as the first entries of the WORK 153*> and IWORK arrays, and no error message related to LWORK or 154*> LIWORK is issued by XERBLA. 155*> \endverbatim 156*> 157*> \param[out] IWORK 158*> \verbatim 159*> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 160*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 161*> \endverbatim 162*> 163*> \param[in] LIWORK 164*> \verbatim 165*> LIWORK is INTEGER 166*> The dimension of the array IWORK. 167*> If N <= 1, LIWORK >= 1. 168*> If JOBZ = 'N' and N > 1, LIWORK >= 1. 169*> If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 170*> 171*> If LIWORK = -1, then a workspace query is assumed; the 172*> routine only calculates the optimal sizes of the WORK and 173*> IWORK arrays, returns these values as the first entries of 174*> the WORK and IWORK arrays, and no error message related to 175*> LWORK or LIWORK is issued by XERBLA. 176*> \endverbatim 177*> 178*> \param[out] INFO 179*> \verbatim 180*> INFO is INTEGER 181*> = 0: successful exit 182*> < 0: if INFO = -i, the i-th argument had an illegal value 183*> > 0: DPOTRF or DSYEVD returned an error code: 184*> <= N: if INFO = i and JOBZ = 'N', then the algorithm 185*> failed to converge; i off-diagonal elements of an 186*> intermediate tridiagonal form did not converge to 187*> zero; 188*> if INFO = i and JOBZ = 'V', then the algorithm 189*> failed to compute an eigenvalue while working on 190*> the submatrix lying in rows and columns INFO/(N+1) 191*> through mod(INFO,N+1); 192*> > N: if INFO = N + i, for 1 <= i <= N, then the leading 193*> minor of order i of B is not positive definite. 194*> The factorization of B could not be completed and 195*> no eigenvalues or eigenvectors were computed. 196*> \endverbatim 197* 198* Authors: 199* ======== 200* 201*> \author Univ. of Tennessee 202*> \author Univ. of California Berkeley 203*> \author Univ. of Colorado Denver 204*> \author NAG Ltd. 205* 206*> \ingroup doubleSYeigen 207* 208*> \par Further Details: 209* ===================== 210*> 211*> \verbatim 212*> 213*> Modified so that no backsubstitution is performed if DSYEVD fails to 214*> converge (NEIG in old code could be greater than N causing out of 215*> bounds reference to A - reported by Ralf Meyer). Also corrected the 216*> description of INFO and the test on ITYPE. Sven, 16 Feb 05. 217*> \endverbatim 218* 219*> \par Contributors: 220* ================== 221*> 222*> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 223*> 224* ===================================================================== 225 SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, 226 $ LWORK, IWORK, LIWORK, INFO ) 227* 228* -- LAPACK driver routine -- 229* -- LAPACK is a software package provided by Univ. of Tennessee, -- 230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 231* 232* .. Scalar Arguments .. 233 CHARACTER JOBZ, UPLO 234 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N 235* .. 236* .. Array Arguments .. 237 INTEGER IWORK( * ) 238 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) 239* .. 240* 241* ===================================================================== 242* 243* .. Parameters .. 244 DOUBLE PRECISION ONE 245 PARAMETER ( ONE = 1.0D+0 ) 246* .. 247* .. Local Scalars .. 248 LOGICAL LQUERY, UPPER, WANTZ 249 CHARACTER TRANS 250 INTEGER LIOPT, LIWMIN, LOPT, LWMIN 251* .. 252* .. External Functions .. 253 LOGICAL LSAME 254 EXTERNAL LSAME 255* .. 256* .. External Subroutines .. 257 EXTERNAL DPOTRF, DSYEVD, DSYGST, DTRMM, DTRSM, XERBLA 258* .. 259* .. Intrinsic Functions .. 260 INTRINSIC DBLE, MAX 261* .. 262* .. Executable Statements .. 263* 264* Test the input parameters. 265* 266 WANTZ = LSAME( JOBZ, 'V' ) 267 UPPER = LSAME( UPLO, 'U' ) 268 LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 269* 270 INFO = 0 271 IF( N.LE.1 ) THEN 272 LIWMIN = 1 273 LWMIN = 1 274 ELSE IF( WANTZ ) THEN 275 LIWMIN = 3 + 5*N 276 LWMIN = 1 + 6*N + 2*N**2 277 ELSE 278 LIWMIN = 1 279 LWMIN = 2*N + 1 280 END IF 281 LOPT = LWMIN 282 LIOPT = LIWMIN 283 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 284 INFO = -1 285 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 286 INFO = -2 287 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 288 INFO = -3 289 ELSE IF( N.LT.0 ) THEN 290 INFO = -4 291 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 292 INFO = -6 293 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 294 INFO = -8 295 END IF 296* 297 IF( INFO.EQ.0 ) THEN 298 WORK( 1 ) = LOPT 299 IWORK( 1 ) = LIOPT 300* 301 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 302 INFO = -11 303 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 304 INFO = -13 305 END IF 306 END IF 307* 308 IF( INFO.NE.0 ) THEN 309 CALL XERBLA( 'DSYGVD', -INFO ) 310 RETURN 311 ELSE IF( LQUERY ) THEN 312 RETURN 313 END IF 314* 315* Quick return if possible 316* 317 IF( N.EQ.0 ) 318 $ RETURN 319* 320* Form a Cholesky factorization of B. 321* 322 CALL DPOTRF( UPLO, N, B, LDB, INFO ) 323 IF( INFO.NE.0 ) THEN 324 INFO = N + INFO 325 RETURN 326 END IF 327* 328* Transform problem to standard eigenvalue problem and solve. 329* 330 CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 331 CALL DSYEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, IWORK, LIWORK, 332 $ INFO ) 333 LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) 334 LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) 335* 336 IF( WANTZ .AND. INFO.EQ.0 ) THEN 337* 338* Backtransform eigenvectors to the original problem. 339* 340 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN 341* 342* For A*x=(lambda)*B*x and A*B*x=(lambda)*x; 343* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y 344* 345 IF( UPPER ) THEN 346 TRANS = 'N' 347 ELSE 348 TRANS = 'T' 349 END IF 350* 351 CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE, 352 $ B, LDB, A, LDA ) 353* 354 ELSE IF( ITYPE.EQ.3 ) THEN 355* 356* For B*A*x=(lambda)*x; 357* backtransform eigenvectors: x = L*y or U**T*y 358* 359 IF( UPPER ) THEN 360 TRANS = 'T' 361 ELSE 362 TRANS = 'N' 363 END IF 364* 365 CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, ONE, 366 $ B, LDB, A, LDA ) 367 END IF 368 END IF 369* 370 WORK( 1 ) = LOPT 371 IWORK( 1 ) = LIOPT 372* 373 RETURN 374* 375* End of DSYGVD 376* 377 END 378