 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, 00002 $ LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBZ, UPLO 00011 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LRWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 DOUBLE PRECISION RWORK( * ), W( * ) 00016 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * ZHEGVD computes all the eigenvalues, and optionally, the eigenvectors 00023 * of a complex generalized Hermitian-definite eigenproblem, of the form 00024 * A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x. Here A and 00025 * B are assumed to be Hermitian and B is also positive definite. 00026 * If eigenvectors are desired, it uses a divide and conquer algorithm. 00027 * 00028 * The divide and conquer algorithm makes very mild assumptions about 00029 * floating point arithmetic. It will work on machines with a guard 00030 * digit in add/subtract, or on those binary machines without guard 00031 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00032 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00033 * without guard digits, but we know of none. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * ITYPE (input) INTEGER 00039 * Specifies the problem type to be solved: 00040 * = 1: A*x = (lambda)*B*x 00041 * = 2: A*B*x = (lambda)*x 00042 * = 3: B*A*x = (lambda)*x 00043 * 00044 * JOBZ (input) CHARACTER*1 00045 * = 'N': Compute eigenvalues only; 00046 * = 'V': Compute eigenvalues and eigenvectors. 00047 * 00048 * UPLO (input) CHARACTER*1 00049 * = 'U': Upper triangles of A and B are stored; 00050 * = 'L': Lower triangles of A and B are stored. 00051 * 00052 * N (input) INTEGER 00053 * The order of the matrices A and B. N >= 0. 00054 * 00055 * A (input/output) COMPLEX*16 array, dimension (LDA, N) 00056 * On entry, the Hermitian matrix A. If UPLO = 'U', the 00057 * leading N-by-N upper triangular part of A contains the 00058 * upper triangular part of the matrix A. If UPLO = 'L', 00059 * the leading N-by-N lower triangular part of A contains 00060 * the lower triangular part of the matrix A. 00061 * 00062 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00063 * matrix Z of eigenvectors. The eigenvectors are normalized 00064 * as follows: 00065 * if ITYPE = 1 or 2, Z**H*B*Z = I; 00066 * if ITYPE = 3, Z**H*inv(B)*Z = I. 00067 * If JOBZ = 'N', then on exit the upper triangle (if UPLO='U') 00068 * or the lower triangle (if UPLO='L') of A, including the 00069 * diagonal, is destroyed. 00070 * 00071 * LDA (input) INTEGER 00072 * The leading dimension of the array A. LDA >= max(1,N). 00073 * 00074 * B (input/output) COMPLEX*16 array, dimension (LDB, N) 00075 * On entry, the Hermitian matrix B. If UPLO = 'U', the 00076 * leading N-by-N upper triangular part of B contains the 00077 * upper triangular part of the matrix B. If UPLO = 'L', 00078 * the leading N-by-N lower triangular part of B contains 00079 * the lower triangular part of the matrix B. 00080 * 00081 * On exit, if INFO <= N, the part of B containing the matrix is 00082 * overwritten by the triangular factor U or L from the Cholesky 00083 * factorization B = U**H*U or B = L*L**H. 00084 * 00085 * LDB (input) INTEGER 00086 * The leading dimension of the array B. LDB >= max(1,N). 00087 * 00088 * W (output) DOUBLE PRECISION array, dimension (N) 00089 * If INFO = 0, the eigenvalues in ascending order. 00090 * 00091 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00092 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00093 * 00094 * LWORK (input) INTEGER 00095 * The length of the array WORK. 00096 * If N <= 1, LWORK >= 1. 00097 * If JOBZ = 'N' and N > 1, LWORK >= N + 1. 00098 * If JOBZ = 'V' and N > 1, LWORK >= 2*N + N**2. 00099 * 00100 * If LWORK = -1, then a workspace query is assumed; the routine 00101 * only calculates the optimal sizes of the WORK, RWORK and 00102 * IWORK arrays, returns these values as the first entries of 00103 * the WORK, RWORK and IWORK arrays, and no error message 00104 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00105 * 00106 * RWORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 00107 * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00108 * 00109 * LRWORK (input) INTEGER 00110 * The dimension of the array RWORK. 00111 * If N <= 1, LRWORK >= 1. 00112 * If JOBZ = 'N' and N > 1, LRWORK >= N. 00113 * If JOBZ = 'V' and N > 1, LRWORK >= 1 + 5*N + 2*N**2. 00114 * 00115 * If LRWORK = -1, then a workspace query is assumed; the 00116 * routine only calculates the optimal sizes of the WORK, RWORK 00117 * and IWORK arrays, returns these values as the first entries 00118 * of the WORK, RWORK and IWORK arrays, and no error message 00119 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00120 * 00121 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00122 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00123 * 00124 * LIWORK (input) INTEGER 00125 * The dimension of the array IWORK. 00126 * If N <= 1, LIWORK >= 1. 00127 * If JOBZ = 'N' and N > 1, LIWORK >= 1. 00128 * If JOBZ = 'V' and N > 1, LIWORK >= 3 + 5*N. 00129 * 00130 * If LIWORK = -1, then a workspace query is assumed; the 00131 * routine only calculates the optimal sizes of the WORK, RWORK 00132 * and IWORK arrays, returns these values as the first entries 00133 * of the WORK, RWORK and IWORK arrays, and no error message 00134 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00135 * 00136 * INFO (output) INTEGER 00137 * = 0: successful exit 00138 * < 0: if INFO = -i, the i-th argument had an illegal value 00139 * > 0: ZPOTRF or ZHEEVD returned an error code: 00140 * <= N: if INFO = i and JOBZ = 'N', then the algorithm 00141 * failed to converge; i off-diagonal elements of an 00142 * intermediate tridiagonal form did not converge to 00143 * zero; 00144 * if INFO = i and JOBZ = 'V', then the algorithm 00145 * failed to compute an eigenvalue while working on 00146 * the submatrix lying in rows and columns INFO/(N+1) 00147 * through mod(INFO,N+1); 00148 * > N: if INFO = N + i, for 1 <= i <= N, then the leading 00149 * minor of order i of B is not positive definite. 00150 * The factorization of B could not be completed and 00151 * no eigenvalues or eigenvectors were computed. 00152 * 00153 * Further Details 00154 * =============== 00155 * 00156 * Based on contributions by 00157 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 00158 * 00159 * Modified so that no backsubstitution is performed if ZHEEVD fails to 00160 * converge (NEIG in old code could be greater than N causing out of 00161 * bounds reference to A - reported by Ralf Meyer). Also corrected the 00162 * description of INFO and the test on ITYPE. Sven, 16 Feb 05. 00163 * ===================================================================== 00164 * 00165 * .. Parameters .. 00166 COMPLEX*16 CONE 00167 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) 00168 * .. 00169 * .. Local Scalars .. 00170 LOGICAL LQUERY, UPPER, WANTZ 00171 CHARACTER TRANS 00172 INTEGER LIOPT, LIWMIN, LOPT, LROPT, LRWMIN, LWMIN 00173 * .. 00174 * .. External Functions .. 00175 LOGICAL LSAME 00176 EXTERNAL LSAME 00177 * .. 00178 * .. External Subroutines .. 00179 EXTERNAL XERBLA, ZHEEVD, ZHEGST, ZPOTRF, ZTRMM, ZTRSM 00180 * .. 00181 * .. Intrinsic Functions .. 00182 INTRINSIC DBLE, MAX 00183 * .. 00184 * .. Executable Statements .. 00185 * 00186 * Test the input parameters. 00187 * 00188 WANTZ = LSAME( JOBZ, 'V' ) 00189 UPPER = LSAME( UPLO, 'U' ) 00190 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00191 * 00192 INFO = 0 00193 IF( N.LE.1 ) THEN 00194 LWMIN = 1 00195 LRWMIN = 1 00196 LIWMIN = 1 00197 ELSE IF( WANTZ ) THEN 00198 LWMIN = 2*N + N*N 00199 LRWMIN = 1 + 5*N + 2*N*N 00200 LIWMIN = 3 + 5*N 00201 ELSE 00202 LWMIN = N + 1 00203 LRWMIN = N 00204 LIWMIN = 1 00205 END IF 00206 LOPT = LWMIN 00207 LROPT = LRWMIN 00208 LIOPT = LIWMIN 00209 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00210 INFO = -1 00211 ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00212 INFO = -2 00213 ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN 00214 INFO = -3 00215 ELSE IF( N.LT.0 ) THEN 00216 INFO = -4 00217 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00218 INFO = -6 00219 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00220 INFO = -8 00221 END IF 00222 * 00223 IF( INFO.EQ.0 ) THEN 00224 WORK( 1 ) = LOPT 00225 RWORK( 1 ) = LROPT 00226 IWORK( 1 ) = LIOPT 00227 * 00228 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00229 INFO = -11 00230 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00231 INFO = -13 00232 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00233 INFO = -15 00234 END IF 00235 END IF 00236 * 00237 IF( INFO.NE.0 ) THEN 00238 CALL XERBLA( 'ZHEGVD', -INFO ) 00239 RETURN 00240 ELSE IF( LQUERY ) THEN 00241 RETURN 00242 END IF 00243 * 00244 * Quick return if possible 00245 * 00246 IF( N.EQ.0 ) 00247 $ RETURN 00248 * 00249 * Form a Cholesky factorization of B. 00250 * 00251 CALL ZPOTRF( UPLO, N, B, LDB, INFO ) 00252 IF( INFO.NE.0 ) THEN 00253 INFO = N + INFO 00254 RETURN 00255 END IF 00256 * 00257 * Transform problem to standard eigenvalue problem and solve. 00258 * 00259 CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00260 CALL ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, LRWORK, 00261 $ IWORK, LIWORK, INFO ) 00262 LOPT = MAX( DBLE( LOPT ), DBLE( WORK( 1 ) ) ) 00263 LROPT = MAX( DBLE( LROPT ), DBLE( RWORK( 1 ) ) ) 00264 LIOPT = MAX( DBLE( LIOPT ), DBLE( IWORK( 1 ) ) ) 00265 * 00266 IF( WANTZ .AND. INFO.EQ.0 ) THEN 00267 * 00268 * Backtransform eigenvectors to the original problem. 00269 * 00270 IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN 00271 * 00272 * For A*x=(lambda)*B*x and A*B*x=(lambda)*x; 00273 * backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y 00274 * 00275 IF( UPPER ) THEN 00276 TRANS = 'N' 00277 ELSE 00278 TRANS = 'C' 00279 END IF 00280 * 00281 CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE, 00282 $ B, LDB, A, LDA ) 00283 * 00284 ELSE IF( ITYPE.EQ.3 ) THEN 00285 * 00286 * For B*A*x=(lambda)*x; 00287 * backtransform eigenvectors: x = L*y or U**H *y 00288 * 00289 IF( UPPER ) THEN 00290 TRANS = 'C' 00291 ELSE 00292 TRANS = 'N' 00293 END IF 00294 * 00295 CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, N, CONE, 00296 $ B, LDB, A, LDA ) 00297 END IF 00298 END IF 00299 * 00300 WORK( 1 ) = LOPT 00301 RWORK( 1 ) = LROPT 00302 IWORK( 1 ) = LIOPT 00303 * 00304 RETURN 00305 * 00306 * End of ZHEGVD 00307 * 00308 END
1.7.3 
