|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, 00002 $ LRWORK, IWORK, LIWORK, 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, LDA, LIWORK, LRWORK, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER IWORK( * ) 00015 DOUBLE PRECISION RWORK( * ), W( * ) 00016 COMPLEX*16 A( LDA, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a 00023 * complex Hermitian matrix A. If eigenvectors are desired, it uses a 00024 * divide and conquer algorithm. 00025 * 00026 * The divide and conquer algorithm makes very mild assumptions about 00027 * floating point arithmetic. It will work on machines with a guard 00028 * digit in add/subtract, or on those binary machines without guard 00029 * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00030 * Cray-2. It could conceivably fail on hexadecimal or decimal machines 00031 * without guard digits, but we know of none. 00032 * 00033 * Arguments 00034 * ========= 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 triangle of A is stored; 00042 * = 'L': Lower triangle of A is stored. 00043 * 00044 * N (input) INTEGER 00045 * The order of the matrix A. 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 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00054 * orthonormal eigenvectors of the matrix A. 00055 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00056 * or the upper triangle (if UPLO='U') of A, including the 00057 * diagonal, is destroyed. 00058 * 00059 * LDA (input) INTEGER 00060 * The leading dimension of the array A. LDA >= max(1,N). 00061 * 00062 * W (output) DOUBLE PRECISION array, dimension (N) 00063 * If INFO = 0, the eigenvalues in ascending order. 00064 * 00065 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00066 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00067 * 00068 * LWORK (input) INTEGER 00069 * The length of the array WORK. 00070 * If N <= 1, LWORK must be at least 1. 00071 * If JOBZ = 'N' and N > 1, LWORK must be at least N + 1. 00072 * If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2. 00073 * 00074 * If LWORK = -1, then a workspace query is assumed; the routine 00075 * only calculates the optimal sizes of the WORK, RWORK and 00076 * IWORK arrays, returns these values as the first entries of 00077 * the WORK, RWORK and IWORK arrays, and no error message 00078 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00079 * 00080 * RWORK (workspace/output) DOUBLE PRECISION array, 00081 * dimension (LRWORK) 00082 * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00083 * 00084 * LRWORK (input) INTEGER 00085 * The dimension of the array RWORK. 00086 * If N <= 1, LRWORK must be at least 1. 00087 * If JOBZ = 'N' and N > 1, LRWORK must be at least N. 00088 * If JOBZ = 'V' and N > 1, LRWORK must be at least 00089 * 1 + 5*N + 2*N**2. 00090 * 00091 * If LRWORK = -1, then a workspace query is assumed; the 00092 * routine only calculates the optimal sizes of the WORK, RWORK 00093 * and IWORK arrays, returns these values as the first entries 00094 * of the WORK, RWORK and IWORK arrays, and no error message 00095 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00096 * 00097 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00098 * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00099 * 00100 * LIWORK (input) INTEGER 00101 * The dimension of the array IWORK. 00102 * If N <= 1, LIWORK must be at least 1. 00103 * If JOBZ = 'N' and N > 1, LIWORK must be at least 1. 00104 * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 00105 * 00106 * If LIWORK = -1, then a workspace query is assumed; the 00107 * routine only calculates the optimal sizes of the WORK, RWORK 00108 * and IWORK arrays, returns these values as the first entries 00109 * of the WORK, RWORK and IWORK arrays, and no error message 00110 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00111 * 00112 * INFO (output) INTEGER 00113 * = 0: successful exit 00114 * < 0: if INFO = -i, the i-th argument had an illegal value 00115 * > 0: if INFO = i and JOBZ = 'N', then the algorithm failed 00116 * to converge; i off-diagonal elements of an intermediate 00117 * tridiagonal form did not converge to zero; 00118 * if INFO = i and JOBZ = 'V', then the algorithm failed 00119 * to compute an eigenvalue while working on the submatrix 00120 * lying in rows and columns INFO/(N+1) through 00121 * mod(INFO,N+1). 00122 * 00123 * Further Details 00124 * =============== 00125 * 00126 * Based on contributions by 00127 * Jeff Rutter, Computer Science Division, University of California 00128 * at Berkeley, USA 00129 * 00130 * Modified description of INFO. Sven, 16 Feb 05. 00131 * ===================================================================== 00132 * 00133 * .. Parameters .. 00134 DOUBLE PRECISION ZERO, ONE 00135 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00136 COMPLEX*16 CONE 00137 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 00138 * .. 00139 * .. Local Scalars .. 00140 LOGICAL LOWER, LQUERY, WANTZ 00141 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2, 00142 $ INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK, 00143 $ LLWRK2, LOPT, LROPT, LRWMIN, LWMIN 00144 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00145 $ SMLNUM 00146 * .. 00147 * .. External Functions .. 00148 LOGICAL LSAME 00149 INTEGER ILAENV 00150 DOUBLE PRECISION DLAMCH, ZLANHE 00151 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE 00152 * .. 00153 * .. External Subroutines .. 00154 EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLACPY, ZLASCL, 00155 $ ZSTEDC, ZUNMTR 00156 * .. 00157 * .. Intrinsic Functions .. 00158 INTRINSIC MAX, SQRT 00159 * .. 00160 * .. Executable Statements .. 00161 * 00162 * Test the input parameters. 00163 * 00164 WANTZ = LSAME( JOBZ, 'V' ) 00165 LOWER = LSAME( UPLO, 'L' ) 00166 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00167 * 00168 INFO = 0 00169 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00170 INFO = -1 00171 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00172 INFO = -2 00173 ELSE IF( N.LT.0 ) THEN 00174 INFO = -3 00175 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00176 INFO = -5 00177 END IF 00178 * 00179 IF( INFO.EQ.0 ) THEN 00180 IF( N.LE.1 ) THEN 00181 LWMIN = 1 00182 LRWMIN = 1 00183 LIWMIN = 1 00184 LOPT = LWMIN 00185 LROPT = LRWMIN 00186 LIOPT = LIWMIN 00187 ELSE 00188 IF( WANTZ ) THEN 00189 LWMIN = 2*N + N*N 00190 LRWMIN = 1 + 5*N + 2*N**2 00191 LIWMIN = 3 + 5*N 00192 ELSE 00193 LWMIN = N + 1 00194 LRWMIN = N 00195 LIWMIN = 1 00196 END IF 00197 LOPT = MAX( LWMIN, N + 00198 $ ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) ) 00199 LROPT = LRWMIN 00200 LIOPT = LIWMIN 00201 END IF 00202 WORK( 1 ) = LOPT 00203 RWORK( 1 ) = LROPT 00204 IWORK( 1 ) = LIOPT 00205 * 00206 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00207 INFO = -8 00208 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00209 INFO = -10 00210 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00211 INFO = -12 00212 END IF 00213 END IF 00214 * 00215 IF( INFO.NE.0 ) THEN 00216 CALL XERBLA( 'ZHEEVD', -INFO ) 00217 RETURN 00218 ELSE IF( LQUERY ) THEN 00219 RETURN 00220 END IF 00221 * 00222 * Quick return if possible 00223 * 00224 IF( N.EQ.0 ) 00225 $ RETURN 00226 * 00227 IF( N.EQ.1 ) THEN 00228 W( 1 ) = A( 1, 1 ) 00229 IF( WANTZ ) 00230 $ A( 1, 1 ) = CONE 00231 RETURN 00232 END IF 00233 * 00234 * Get machine constants. 00235 * 00236 SAFMIN = DLAMCH( 'Safe minimum' ) 00237 EPS = DLAMCH( 'Precision' ) 00238 SMLNUM = SAFMIN / EPS 00239 BIGNUM = ONE / SMLNUM 00240 RMIN = SQRT( SMLNUM ) 00241 RMAX = SQRT( BIGNUM ) 00242 * 00243 * Scale matrix to allowable range, if necessary. 00244 * 00245 ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK ) 00246 ISCALE = 0 00247 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00248 ISCALE = 1 00249 SIGMA = RMIN / ANRM 00250 ELSE IF( ANRM.GT.RMAX ) THEN 00251 ISCALE = 1 00252 SIGMA = RMAX / ANRM 00253 END IF 00254 IF( ISCALE.EQ.1 ) 00255 $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 00256 * 00257 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form. 00258 * 00259 INDE = 1 00260 INDTAU = 1 00261 INDWRK = INDTAU + N 00262 INDRWK = INDE + N 00263 INDWK2 = INDWRK + N*N 00264 LLWORK = LWORK - INDWRK + 1 00265 LLWRK2 = LWORK - INDWK2 + 1 00266 LLRWK = LRWORK - INDRWK + 1 00267 CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ), 00268 $ WORK( INDWRK ), LLWORK, IINFO ) 00269 * 00270 * For eigenvalues only, call DSTERF. For eigenvectors, first call 00271 * ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 00272 * tridiagonal matrix, then call ZUNMTR to multiply it to the 00273 * Householder transformations represented as Householder vectors in 00274 * A. 00275 * 00276 IF( .NOT.WANTZ ) THEN 00277 CALL DSTERF( N, W, RWORK( INDE ), INFO ) 00278 ELSE 00279 CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N, 00280 $ WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK, 00281 $ IWORK, LIWORK, INFO ) 00282 CALL ZUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ), 00283 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO ) 00284 CALL ZLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA ) 00285 END IF 00286 * 00287 * If matrix was scaled, then rescale eigenvalues appropriately. 00288 * 00289 IF( ISCALE.EQ.1 ) THEN 00290 IF( INFO.EQ.0 ) THEN 00291 IMAX = N 00292 ELSE 00293 IMAX = INFO - 1 00294 END IF 00295 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00296 END IF 00297 * 00298 WORK( 1 ) = LOPT 00299 RWORK( 1 ) = LROPT 00300 IWORK( 1 ) = LIOPT 00301 * 00302 RETURN 00303 * 00304 * End of ZHEEVD 00305 * 00306 END
1.7.3 
