![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZHEEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b> 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHEEVD + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevd.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevd.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevd.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, 00022 * LRWORK, IWORK, LIWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER JOBZ, UPLO 00026 * INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * INTEGER IWORK( * ) 00030 * DOUBLE PRECISION RWORK( * ), W( * ) 00031 * COMPLEX*16 A( LDA, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a 00041 *> complex Hermitian matrix A. If eigenvectors are desired, it uses a 00042 *> divide and conquer algorithm. 00043 *> 00044 *> The divide and conquer algorithm makes very mild assumptions about 00045 *> floating point arithmetic. It will work on machines with a guard 00046 *> digit in add/subtract, or on those binary machines without guard 00047 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or 00048 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines 00049 *> without guard digits, but we know of none. 00050 *> \endverbatim 00051 * 00052 * Arguments: 00053 * ========== 00054 * 00055 *> \param[in] JOBZ 00056 *> \verbatim 00057 *> JOBZ is CHARACTER*1 00058 *> = 'N': Compute eigenvalues only; 00059 *> = 'V': Compute eigenvalues and eigenvectors. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] UPLO 00063 *> \verbatim 00064 *> UPLO is CHARACTER*1 00065 *> = 'U': Upper triangle of A is stored; 00066 *> = 'L': Lower triangle of A is stored. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] N 00070 *> \verbatim 00071 *> N is INTEGER 00072 *> The order of the matrix A. N >= 0. 00073 *> \endverbatim 00074 *> 00075 *> \param[in,out] A 00076 *> \verbatim 00077 *> A is COMPLEX*16 array, dimension (LDA, N) 00078 *> On entry, the Hermitian matrix A. If UPLO = 'U', the 00079 *> leading N-by-N upper triangular part of A contains the 00080 *> upper triangular part of the matrix A. If UPLO = 'L', 00081 *> the leading N-by-N lower triangular part of A contains 00082 *> the lower triangular part of the matrix A. 00083 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the 00084 *> orthonormal eigenvectors of the matrix A. 00085 *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') 00086 *> or the upper triangle (if UPLO='U') of A, including the 00087 *> diagonal, is destroyed. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] LDA 00091 *> \verbatim 00092 *> LDA is INTEGER 00093 *> The leading dimension of the array A. LDA >= max(1,N). 00094 *> \endverbatim 00095 *> 00096 *> \param[out] W 00097 *> \verbatim 00098 *> W is DOUBLE PRECISION array, dimension (N) 00099 *> If INFO = 0, the eigenvalues in ascending order. 00100 *> \endverbatim 00101 *> 00102 *> \param[out] WORK 00103 *> \verbatim 00104 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) 00105 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] LWORK 00109 *> \verbatim 00110 *> LWORK is INTEGER 00111 *> The length of the array WORK. 00112 *> If N <= 1, LWORK must be at least 1. 00113 *> If JOBZ = 'N' and N > 1, LWORK must be at least N + 1. 00114 *> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2. 00115 *> 00116 *> If LWORK = -1, then a workspace query is assumed; the routine 00117 *> only calculates the optimal sizes of the WORK, RWORK and 00118 *> IWORK arrays, returns these values as the first entries of 00119 *> the WORK, RWORK and IWORK arrays, and no error message 00120 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00121 *> \endverbatim 00122 *> 00123 *> \param[out] RWORK 00124 *> \verbatim 00125 *> RWORK is DOUBLE PRECISION array, 00126 *> dimension (LRWORK) 00127 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. 00128 *> \endverbatim 00129 *> 00130 *> \param[in] LRWORK 00131 *> \verbatim 00132 *> LRWORK is INTEGER 00133 *> The dimension of the array RWORK. 00134 *> If N <= 1, LRWORK must be at least 1. 00135 *> If JOBZ = 'N' and N > 1, LRWORK must be at least N. 00136 *> If JOBZ = 'V' and N > 1, LRWORK must be at least 00137 *> 1 + 5*N + 2*N**2. 00138 *> 00139 *> If LRWORK = -1, then a workspace query is assumed; the 00140 *> routine only calculates the optimal sizes of the WORK, RWORK 00141 *> and IWORK arrays, returns these values as the first entries 00142 *> of the WORK, RWORK and IWORK arrays, and no error message 00143 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00144 *> \endverbatim 00145 *> 00146 *> \param[out] IWORK 00147 *> \verbatim 00148 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) 00149 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00150 *> \endverbatim 00151 *> 00152 *> \param[in] LIWORK 00153 *> \verbatim 00154 *> LIWORK is INTEGER 00155 *> The dimension of the array IWORK. 00156 *> If N <= 1, LIWORK must be at least 1. 00157 *> If JOBZ = 'N' and N > 1, LIWORK must be at least 1. 00158 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N. 00159 *> 00160 *> If LIWORK = -1, then a workspace query is assumed; the 00161 *> routine only calculates the optimal sizes of the WORK, RWORK 00162 *> and IWORK arrays, returns these values as the first entries 00163 *> of the WORK, RWORK and IWORK arrays, and no error message 00164 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00165 *> \endverbatim 00166 *> 00167 *> \param[out] INFO 00168 *> \verbatim 00169 *> INFO is INTEGER 00170 *> = 0: successful exit 00171 *> < 0: if INFO = -i, the i-th argument had an illegal value 00172 *> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed 00173 *> to converge; i off-diagonal elements of an intermediate 00174 *> tridiagonal form did not converge to zero; 00175 *> if INFO = i and JOBZ = 'V', then the algorithm failed 00176 *> to compute an eigenvalue while working on the submatrix 00177 *> lying in rows and columns INFO/(N+1) through 00178 *> mod(INFO,N+1). 00179 *> \endverbatim 00180 * 00181 * Authors: 00182 * ======== 00183 * 00184 *> \author Univ. of Tennessee 00185 *> \author Univ. of California Berkeley 00186 *> \author Univ. of Colorado Denver 00187 *> \author NAG Ltd. 00188 * 00189 *> \date November 2011 00190 * 00191 *> \ingroup complex16HEeigen 00192 * 00193 *> \par Further Details: 00194 * ===================== 00195 *> 00196 *> Modified description of INFO. Sven, 16 Feb 05. 00197 * 00198 *> \par Contributors: 00199 * ================== 00200 *> 00201 *> Jeff Rutter, Computer Science Division, University of California 00202 *> at Berkeley, USA 00203 *> 00204 * ===================================================================== 00205 SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, 00206 $ LRWORK, IWORK, LIWORK, INFO ) 00207 * 00208 * -- LAPACK driver routine (version 3.4.0) -- 00209 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00211 * November 2011 00212 * 00213 * .. Scalar Arguments .. 00214 CHARACTER JOBZ, UPLO 00215 INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N 00216 * .. 00217 * .. Array Arguments .. 00218 INTEGER IWORK( * ) 00219 DOUBLE PRECISION RWORK( * ), W( * ) 00220 COMPLEX*16 A( LDA, * ), WORK( * ) 00221 * .. 00222 * 00223 * ===================================================================== 00224 * 00225 * .. Parameters .. 00226 DOUBLE PRECISION ZERO, ONE 00227 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00228 COMPLEX*16 CONE 00229 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 00230 * .. 00231 * .. Local Scalars .. 00232 LOGICAL LOWER, LQUERY, WANTZ 00233 INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2, 00234 $ INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK, 00235 $ LLWRK2, LOPT, LROPT, LRWMIN, LWMIN 00236 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, 00237 $ SMLNUM 00238 * .. 00239 * .. External Functions .. 00240 LOGICAL LSAME 00241 INTEGER ILAENV 00242 DOUBLE PRECISION DLAMCH, ZLANHE 00243 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE 00244 * .. 00245 * .. External Subroutines .. 00246 EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLACPY, ZLASCL, 00247 $ ZSTEDC, ZUNMTR 00248 * .. 00249 * .. Intrinsic Functions .. 00250 INTRINSIC MAX, SQRT 00251 * .. 00252 * .. Executable Statements .. 00253 * 00254 * Test the input parameters. 00255 * 00256 WANTZ = LSAME( JOBZ, 'V' ) 00257 LOWER = LSAME( UPLO, 'L' ) 00258 LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 ) 00259 * 00260 INFO = 0 00261 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00262 INFO = -1 00263 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00264 INFO = -2 00265 ELSE IF( N.LT.0 ) THEN 00266 INFO = -3 00267 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00268 INFO = -5 00269 END IF 00270 * 00271 IF( INFO.EQ.0 ) THEN 00272 IF( N.LE.1 ) THEN 00273 LWMIN = 1 00274 LRWMIN = 1 00275 LIWMIN = 1 00276 LOPT = LWMIN 00277 LROPT = LRWMIN 00278 LIOPT = LIWMIN 00279 ELSE 00280 IF( WANTZ ) THEN 00281 LWMIN = 2*N + N*N 00282 LRWMIN = 1 + 5*N + 2*N**2 00283 LIWMIN = 3 + 5*N 00284 ELSE 00285 LWMIN = N + 1 00286 LRWMIN = N 00287 LIWMIN = 1 00288 END IF 00289 LOPT = MAX( LWMIN, N + 00290 $ ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) ) 00291 LROPT = LRWMIN 00292 LIOPT = LIWMIN 00293 END IF 00294 WORK( 1 ) = LOPT 00295 RWORK( 1 ) = LROPT 00296 IWORK( 1 ) = LIOPT 00297 * 00298 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00299 INFO = -8 00300 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00301 INFO = -10 00302 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00303 INFO = -12 00304 END IF 00305 END IF 00306 * 00307 IF( INFO.NE.0 ) THEN 00308 CALL XERBLA( 'ZHEEVD', -INFO ) 00309 RETURN 00310 ELSE IF( LQUERY ) THEN 00311 RETURN 00312 END IF 00313 * 00314 * Quick return if possible 00315 * 00316 IF( N.EQ.0 ) 00317 $ RETURN 00318 * 00319 IF( N.EQ.1 ) THEN 00320 W( 1 ) = A( 1, 1 ) 00321 IF( WANTZ ) 00322 $ A( 1, 1 ) = CONE 00323 RETURN 00324 END IF 00325 * 00326 * Get machine constants. 00327 * 00328 SAFMIN = DLAMCH( 'Safe minimum' ) 00329 EPS = DLAMCH( 'Precision' ) 00330 SMLNUM = SAFMIN / EPS 00331 BIGNUM = ONE / SMLNUM 00332 RMIN = SQRT( SMLNUM ) 00333 RMAX = SQRT( BIGNUM ) 00334 * 00335 * Scale matrix to allowable range, if necessary. 00336 * 00337 ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK ) 00338 ISCALE = 0 00339 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00340 ISCALE = 1 00341 SIGMA = RMIN / ANRM 00342 ELSE IF( ANRM.GT.RMAX ) THEN 00343 ISCALE = 1 00344 SIGMA = RMAX / ANRM 00345 END IF 00346 IF( ISCALE.EQ.1 ) 00347 $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO ) 00348 * 00349 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form. 00350 * 00351 INDE = 1 00352 INDTAU = 1 00353 INDWRK = INDTAU + N 00354 INDRWK = INDE + N 00355 INDWK2 = INDWRK + N*N 00356 LLWORK = LWORK - INDWRK + 1 00357 LLWRK2 = LWORK - INDWK2 + 1 00358 LLRWK = LRWORK - INDRWK + 1 00359 CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ), 00360 $ WORK( INDWRK ), LLWORK, IINFO ) 00361 * 00362 * For eigenvalues only, call DSTERF. For eigenvectors, first call 00363 * ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the 00364 * tridiagonal matrix, then call ZUNMTR to multiply it to the 00365 * Householder transformations represented as Householder vectors in 00366 * A. 00367 * 00368 IF( .NOT.WANTZ ) THEN 00369 CALL DSTERF( N, W, RWORK( INDE ), INFO ) 00370 ELSE 00371 CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N, 00372 $ WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK, 00373 $ IWORK, LIWORK, INFO ) 00374 CALL ZUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ), 00375 $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO ) 00376 CALL ZLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA ) 00377 END IF 00378 * 00379 * If matrix was scaled, then rescale eigenvalues appropriately. 00380 * 00381 IF( ISCALE.EQ.1 ) THEN 00382 IF( INFO.EQ.0 ) THEN 00383 IMAX = N 00384 ELSE 00385 IMAX = INFO - 1 00386 END IF 00387 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00388 END IF 00389 * 00390 WORK( 1 ) = LOPT 00391 RWORK( 1 ) = LROPT 00392 IWORK( 1 ) = LIOPT 00393 * 00394 RETURN 00395 * 00396 * End of ZHEEVD 00397 * 00398 END