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