|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, 00002 $ LWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, UPLO 00011 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DSYGV computes all the eigenvalues, and optionally, the eigenvectors 00021 * of a real generalized symmetric-definite eigenproblem, of the form 00022 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. 00023 * Here A and B are assumed to be symmetric and B is also 00024 * positive definite. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * ITYPE (input) INTEGER 00030 * Specifies the problem type to be solved: 00031 * = 1: A*x = (lambda)*B*x 00032 * = 2: A*B*x = (lambda)*x 00033 * = 3: B*A*x = (lambda)*x 00034 * 00035 * JOBZ (input) CHARACTER*1 00036 * = 'N': Compute eigenvalues only; 00037 * = 'V': Compute eigenvalues and eigenvectors. 00038 * 00039 * UPLO (input) CHARACTER*1 00040 * = 'U': Upper triangles of A and B are stored; 00041 * = 'L': Lower triangles of A and B are stored. 00042 * 00043 * N (input) INTEGER 00044 * The order of the matrices A and B. N >= 0. 00045 * 00046 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00047 * On entry, the symmetric matrix A. If UPLO = 'U', the 00048 * leading N-by-N upper triangular part of A contains the 00049 * upper triangular part of the matrix A. If UPLO = 'L', 00050 * the leading N-by-N lower triangular part of A contains 00051 * the lower triangular part of the matrix A. 00052 * 00053 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00054 * matrix Z of eigenvectors. The eigenvectors are normalized 00055 * as follows: 00056 * if ITYPE = 1 or 2, Z**T*B*Z = I; 00057 * if ITYPE = 3, Z**T*inv(B)*Z = I. 00058 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') 00059 * or the lower triangle (if UPLO='L') of A, including the 00060 * diagonal, is destroyed. 00061 * 00062 * LDA (input) INTEGER 00063 * The leading dimension of the array A. LDA >= max(1,N). 00064 * 00065 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00066 * On entry, the symmetric positive definite matrix B. 00067 * If UPLO = 'U', the leading N-by-N upper triangular part of B 00068 * contains the upper triangular part of the matrix B. 00069 * If UPLO = 'L', the leading N-by-N lower triangular part of B 00070 * contains the lower triangular part of the matrix B. 00071 * 00072 * On exit, if INFO <= N, the part of B containing the matrix is 00073 * overwritten by the triangular factor U or L from the Cholesky 00074 * factorization B = U**T*U or B = L*L**T. 00075 * 00076 * LDB (input) INTEGER 00077 * The leading dimension of the array B. LDB >= max(1,N). 00078 * 00079 * W (output) DOUBLE PRECISION array, dimension (N) 00080 * If INFO = 0, the eigenvalues in ascending order. 00081 * 00082 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00083 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00084 * 00085 * LWORK (input) INTEGER 00086 * The length of the array WORK. LWORK >= max(1,3*N-1). 00087 * For optimal efficiency, LWORK >= (NB+2)*N, 00088 * where NB is the blocksize for DSYTRD returned by ILAENV. 00089 * 00090 * If LWORK = -1, then a workspace query is assumed; the routine 00091 * only calculates the optimal size of the WORK array, returns 00092 * this value as the first entry of the WORK array, and no error 00093 * message related to LWORK is issued by XERBLA. 00094 * 00095 * INFO (output) INTEGER 00096 * = 0: successful exit 00097 * < 0: if INFO = -i, the i-th argument had an illegal value 00098 * > 0: DPOTRF or DSYEV returned an error code: 00099 * <= N: if INFO = i, DSYEV failed to converge; 00100 * i off-diagonal elements of an intermediate 00101 * tridiagonal form did not converge to zero; 00102 * > N: if INFO = N + i, for 1 <= i <= N, then the leading 00103 * minor of order i of B is not positive definite. 00104 * The factorization of B could not be completed and 00105 * no eigenvalues or eigenvectors were computed. 00106 * 00107 * ===================================================================== 00108 * 00109 * .. Parameters .. 00110 DOUBLE PRECISION ONE 00111 PARAMETER ( ONE = 1.0D+0 ) 00112 * .. 00113 * .. Local Scalars .. 00114 LOGICAL LQUERY, UPPER, WANTZ 00115 CHARACTER TRANS 00116 INTEGER LWKMIN, LWKOPT, NB, NEIG 00117 * .. 00118 * .. External Functions .. 00119 LOGICAL LSAME 00120 INTEGER ILAENV 00121 EXTERNAL LSAME, ILAENV 00122 * .. 00123 * .. External Subroutines .. 00124 EXTERNAL DPOTRF, DSYEV, DSYGST, DTRMM, DTRSM, XERBLA 00125 * .. 00126 * .. Intrinsic Functions .. 00127 INTRINSIC MAX 00128 * .. 00129 * .. Executable Statements .. 00130 * 00131 * Test the input parameters. 00132 * 00133 WANTZ = LSAME( JOBZ, 'V' ) 00134 UPPER = LSAME( UPLO, 'U' ) 00135 LQUERY = ( LWORK.EQ.-1 ) 00136 * 00137 INFO = 0 00138 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00139 INFO = -1 00140 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00141 INFO = -2 00142 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 00143 INFO = -3 00144 ELSE IF( N.LT.0 ) THEN 00145 INFO = -4 00146 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00147 INFO = -6 00148 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00149 INFO = -8 00150 END IF 00151 * 00152 IF( INFO.EQ.0 ) THEN 00153 LWKMIN = MAX( 1, 3*N - 1 ) 00154 NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 ) 00155 LWKOPT = MAX( LWKMIN, ( NB + 2 )*N ) 00156 WORK( 1 ) = LWKOPT 00157 * 00158 IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00159 INFO = -11 00160 END IF 00161 END IF 00162 * 00163 IF( INFO.NE.0 ) THEN 00164 CALL XERBLA( 'DSYGV ', -INFO ) 00165 RETURN 00166 ELSE IF( LQUERY ) THEN 00167 RETURN 00168 END IF 00169 * 00170 * Quick return if possible 00171 * 00172 IF( N.EQ.0 ) 00173 $ RETURN 00174 * 00175 * Form a Cholesky factorization of B. 00176 * 00177 CALL DPOTRF( UPLO, N, B, LDB, INFO ) 00178 IF( INFO.NE.0 ) THEN 00179 INFO = N + INFO 00180 RETURN 00181 END IF 00182 * 00183 * Transform problem to standard eigenvalue problem and solve. 00184 * 00185 CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00186 CALL DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) 00187 * 00188 IF( WANTZ ) THEN 00189 * 00190 * Backtransform eigenvectors to the original problem. 00191 * 00192 NEIG = N 00193 IF( INFO.GT.0 ) 00194 $ NEIG = INFO - 1 00195 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN 00196 * 00197 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x; 00198 * backtransform eigenvectors: x = inv(L)'*y or inv(U)*y 00199 * 00200 IF( UPPER ) THEN 00201 TRANS = 'N' 00202 ELSE 00203 TRANS = 'T' 00204 END IF 00205 * 00206 CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE, 00207 $ B, LDB, A, LDA ) 00208 * 00209 ELSE IF( ITYPE.EQ.3 ) THEN 00210 * 00211 * For B*A*x=(lambda)*x; 00212 * backtransform eigenvectors: x = L*y or U'*y 00213 * 00214 IF( UPPER ) THEN 00215 TRANS = 'T' 00216 ELSE 00217 TRANS = 'N' 00218 END IF 00219 * 00220 CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE, 00221 $ B, LDB, A, LDA ) 00222 END IF 00223 END IF 00224 * 00225 WORK( 1 ) = LWKOPT 00226 RETURN 00227 * 00228 * End of DSYGV 00229 * 00230 END
1.7.3 
