![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 ZHPEVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpevx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpevx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpevx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, 00022 * ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, 00023 * IFAIL, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBZ, RANGE, UPLO 00027 * INTEGER IL, INFO, IU, LDZ, M, N 00028 * DOUBLE PRECISION ABSTOL, VL, VU 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IFAIL( * ), IWORK( * ) 00032 * DOUBLE PRECISION RWORK( * ), W( * ) 00033 * COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * ) 00034 * .. 00035 * 00036 * 00037 *> \par Purpose: 00038 * ============= 00039 *> 00040 *> \verbatim 00041 *> 00042 *> ZHPEVX computes selected eigenvalues and, optionally, eigenvectors 00043 *> of a complex Hermitian matrix A in packed storage. 00044 *> Eigenvalues/vectors can be selected by specifying either a range of 00045 *> values or a range of indices for the desired eigenvalues. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] JOBZ 00052 *> \verbatim 00053 *> JOBZ is CHARACTER*1 00054 *> = 'N': Compute eigenvalues only; 00055 *> = 'V': Compute eigenvalues and eigenvectors. 00056 *> \endverbatim 00057 *> 00058 *> \param[in] RANGE 00059 *> \verbatim 00060 *> RANGE is CHARACTER*1 00061 *> = 'A': all eigenvalues will be found; 00062 *> = 'V': all eigenvalues in the half-open interval (VL,VU] 00063 *> will be found; 00064 *> = 'I': the IL-th through IU-th eigenvalues will be found. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] UPLO 00068 *> \verbatim 00069 *> UPLO is CHARACTER*1 00070 *> = 'U': Upper triangle of A is stored; 00071 *> = 'L': Lower triangle of A is stored. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] N 00075 *> \verbatim 00076 *> N is INTEGER 00077 *> The order of the matrix A. N >= 0. 00078 *> \endverbatim 00079 *> 00080 *> \param[in,out] AP 00081 *> \verbatim 00082 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2) 00083 *> On entry, the upper or lower triangle of the Hermitian matrix 00084 *> A, packed columnwise in a linear array. The j-th column of A 00085 *> is stored in the array AP as follows: 00086 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00087 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n. 00088 *> 00089 *> On exit, AP is overwritten by values generated during the 00090 *> reduction to tridiagonal form. If UPLO = 'U', the diagonal 00091 *> and first superdiagonal of the tridiagonal matrix T overwrite 00092 *> the corresponding elements of A, and if UPLO = 'L', the 00093 *> diagonal and first subdiagonal of T overwrite the 00094 *> corresponding elements of A. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] VL 00098 *> \verbatim 00099 *> VL is DOUBLE PRECISION 00100 *> \endverbatim 00101 *> 00102 *> \param[in] VU 00103 *> \verbatim 00104 *> VU is DOUBLE PRECISION 00105 *> If RANGE='V', the lower and upper bounds of the interval to 00106 *> be searched for eigenvalues. VL < VU. 00107 *> Not referenced if RANGE = 'A' or 'I'. 00108 *> \endverbatim 00109 *> 00110 *> \param[in] IL 00111 *> \verbatim 00112 *> IL is INTEGER 00113 *> \endverbatim 00114 *> 00115 *> \param[in] IU 00116 *> \verbatim 00117 *> IU is INTEGER 00118 *> If RANGE='I', the indices (in ascending order) of the 00119 *> smallest and largest eigenvalues to be returned. 00120 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00121 *> Not referenced if RANGE = 'A' or 'V'. 00122 *> \endverbatim 00123 *> 00124 *> \param[in] ABSTOL 00125 *> \verbatim 00126 *> ABSTOL is DOUBLE PRECISION 00127 *> The absolute error tolerance for the eigenvalues. 00128 *> An approximate eigenvalue is accepted as converged 00129 *> when it is determined to lie in an interval [a,b] 00130 *> of width less than or equal to 00131 *> 00132 *> ABSTOL + EPS * max( |a|,|b| ) , 00133 *> 00134 *> where EPS is the machine precision. If ABSTOL is less than 00135 *> or equal to zero, then EPS*|T| will be used in its place, 00136 *> where |T| is the 1-norm of the tridiagonal matrix obtained 00137 *> by reducing AP to tridiagonal form. 00138 *> 00139 *> Eigenvalues will be computed most accurately when ABSTOL is 00140 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero. 00141 *> If this routine returns with INFO>0, indicating that some 00142 *> eigenvectors did not converge, try setting ABSTOL to 00143 *> 2*DLAMCH('S'). 00144 *> 00145 *> See "Computing Small Singular Values of Bidiagonal Matrices 00146 *> with Guaranteed High Relative Accuracy," by Demmel and 00147 *> Kahan, LAPACK Working Note #3. 00148 *> \endverbatim 00149 *> 00150 *> \param[out] M 00151 *> \verbatim 00152 *> M is INTEGER 00153 *> The total number of eigenvalues found. 0 <= M <= N. 00154 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00155 *> \endverbatim 00156 *> 00157 *> \param[out] W 00158 *> \verbatim 00159 *> W is DOUBLE PRECISION array, dimension (N) 00160 *> If INFO = 0, the selected eigenvalues in ascending order. 00161 *> \endverbatim 00162 *> 00163 *> \param[out] Z 00164 *> \verbatim 00165 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M)) 00166 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00167 *> contain the orthonormal eigenvectors of the matrix A 00168 *> corresponding to the selected eigenvalues, with the i-th 00169 *> column of Z holding the eigenvector associated with W(i). 00170 *> If an eigenvector fails to converge, then that column of Z 00171 *> contains the latest approximation to the eigenvector, and 00172 *> the index of the eigenvector is returned in IFAIL. 00173 *> If JOBZ = 'N', then Z is not referenced. 00174 *> Note: the user must ensure that at least max(1,M) columns are 00175 *> supplied in the array Z; if RANGE = 'V', the exact value of M 00176 *> is not known in advance and an upper bound must be used. 00177 *> \endverbatim 00178 *> 00179 *> \param[in] LDZ 00180 *> \verbatim 00181 *> LDZ is INTEGER 00182 *> The leading dimension of the array Z. LDZ >= 1, and if 00183 *> JOBZ = 'V', LDZ >= max(1,N). 00184 *> \endverbatim 00185 *> 00186 *> \param[out] WORK 00187 *> \verbatim 00188 *> WORK is COMPLEX*16 array, dimension (2*N) 00189 *> \endverbatim 00190 *> 00191 *> \param[out] RWORK 00192 *> \verbatim 00193 *> RWORK is DOUBLE PRECISION array, dimension (7*N) 00194 *> \endverbatim 00195 *> 00196 *> \param[out] IWORK 00197 *> \verbatim 00198 *> IWORK is INTEGER array, dimension (5*N) 00199 *> \endverbatim 00200 *> 00201 *> \param[out] IFAIL 00202 *> \verbatim 00203 *> IFAIL is INTEGER array, dimension (N) 00204 *> If JOBZ = 'V', then if INFO = 0, the first M elements of 00205 *> IFAIL are zero. If INFO > 0, then IFAIL contains the 00206 *> indices of the eigenvectors that failed to converge. 00207 *> If JOBZ = 'N', then IFAIL is not referenced. 00208 *> \endverbatim 00209 *> 00210 *> \param[out] INFO 00211 *> \verbatim 00212 *> INFO is INTEGER 00213 *> = 0: successful exit 00214 *> < 0: if INFO = -i, the i-th argument had an illegal value 00215 *> > 0: if INFO = i, then i eigenvectors failed to converge. 00216 *> Their indices are stored in array IFAIL. 00217 *> \endverbatim 00218 * 00219 * Authors: 00220 * ======== 00221 * 00222 *> \author Univ. of Tennessee 00223 *> \author Univ. of California Berkeley 00224 *> \author Univ. of Colorado Denver 00225 *> \author NAG Ltd. 00226 * 00227 *> \date November 2011 00228 * 00229 *> \ingroup complex16OTHEReigen 00230 * 00231 * ===================================================================== 00232 SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, 00233 $ ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, 00234 $ IFAIL, INFO ) 00235 * 00236 * -- LAPACK driver routine (version 3.4.0) -- 00237 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00238 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00239 * November 2011 00240 * 00241 * .. Scalar Arguments .. 00242 CHARACTER JOBZ, RANGE, UPLO 00243 INTEGER IL, INFO, IU, LDZ, M, N 00244 DOUBLE PRECISION ABSTOL, VL, VU 00245 * .. 00246 * .. Array Arguments .. 00247 INTEGER IFAIL( * ), IWORK( * ) 00248 DOUBLE PRECISION RWORK( * ), W( * ) 00249 COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * ) 00250 * .. 00251 * 00252 * ===================================================================== 00253 * 00254 * .. Parameters .. 00255 DOUBLE PRECISION ZERO, ONE 00256 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00257 COMPLEX*16 CONE 00258 PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) ) 00259 * .. 00260 * .. Local Scalars .. 00261 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ 00262 CHARACTER ORDER 00263 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL, 00264 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE, 00265 $ ITMP1, J, JJ, NSPLIT 00266 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 00267 $ SIGMA, SMLNUM, TMP1, VLL, VUU 00268 * .. 00269 * .. External Functions .. 00270 LOGICAL LSAME 00271 DOUBLE PRECISION DLAMCH, ZLANHP 00272 EXTERNAL LSAME, DLAMCH, ZLANHP 00273 * .. 00274 * .. External Subroutines .. 00275 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTERF, XERBLA, ZDSCAL, 00276 $ ZHPTRD, ZSTEIN, ZSTEQR, ZSWAP, ZUPGTR, ZUPMTR 00277 * .. 00278 * .. Intrinsic Functions .. 00279 INTRINSIC DBLE, MAX, MIN, SQRT 00280 * .. 00281 * .. Executable Statements .. 00282 * 00283 * Test the input parameters. 00284 * 00285 WANTZ = LSAME( JOBZ, 'V' ) 00286 ALLEIG = LSAME( RANGE, 'A' ) 00287 VALEIG = LSAME( RANGE, 'V' ) 00288 INDEIG = LSAME( RANGE, 'I' ) 00289 * 00290 INFO = 0 00291 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00292 INFO = -1 00293 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00294 INFO = -2 00295 ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) ) 00296 $ THEN 00297 INFO = -3 00298 ELSE IF( N.LT.0 ) THEN 00299 INFO = -4 00300 ELSE 00301 IF( VALEIG ) THEN 00302 IF( N.GT.0 .AND. VU.LE.VL ) 00303 $ INFO = -7 00304 ELSE IF( INDEIG ) THEN 00305 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00306 INFO = -8 00307 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00308 INFO = -9 00309 END IF 00310 END IF 00311 END IF 00312 IF( INFO.EQ.0 ) THEN 00313 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) 00314 $ INFO = -14 00315 END IF 00316 * 00317 IF( INFO.NE.0 ) THEN 00318 CALL XERBLA( 'ZHPEVX', -INFO ) 00319 RETURN 00320 END IF 00321 * 00322 * Quick return if possible 00323 * 00324 M = 0 00325 IF( N.EQ.0 ) 00326 $ RETURN 00327 * 00328 IF( N.EQ.1 ) THEN 00329 IF( ALLEIG .OR. INDEIG ) THEN 00330 M = 1 00331 W( 1 ) = AP( 1 ) 00332 ELSE 00333 IF( VL.LT.DBLE( AP( 1 ) ) .AND. VU.GE.DBLE( AP( 1 ) ) ) THEN 00334 M = 1 00335 W( 1 ) = AP( 1 ) 00336 END IF 00337 END IF 00338 IF( WANTZ ) 00339 $ Z( 1, 1 ) = CONE 00340 RETURN 00341 END IF 00342 * 00343 * Get machine constants. 00344 * 00345 SAFMIN = DLAMCH( 'Safe minimum' ) 00346 EPS = DLAMCH( 'Precision' ) 00347 SMLNUM = SAFMIN / EPS 00348 BIGNUM = ONE / SMLNUM 00349 RMIN = SQRT( SMLNUM ) 00350 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00351 * 00352 * Scale matrix to allowable range, if necessary. 00353 * 00354 ISCALE = 0 00355 ABSTLL = ABSTOL 00356 IF( VALEIG ) THEN 00357 VLL = VL 00358 VUU = VU 00359 ELSE 00360 VLL = ZERO 00361 VUU = ZERO 00362 END IF 00363 ANRM = ZLANHP( 'M', UPLO, N, AP, RWORK ) 00364 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00365 ISCALE = 1 00366 SIGMA = RMIN / ANRM 00367 ELSE IF( ANRM.GT.RMAX ) THEN 00368 ISCALE = 1 00369 SIGMA = RMAX / ANRM 00370 END IF 00371 IF( ISCALE.EQ.1 ) THEN 00372 CALL ZDSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 ) 00373 IF( ABSTOL.GT.0 ) 00374 $ ABSTLL = ABSTOL*SIGMA 00375 IF( VALEIG ) THEN 00376 VLL = VL*SIGMA 00377 VUU = VU*SIGMA 00378 END IF 00379 END IF 00380 * 00381 * Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form. 00382 * 00383 INDD = 1 00384 INDE = INDD + N 00385 INDRWK = INDE + N 00386 INDTAU = 1 00387 INDWRK = INDTAU + N 00388 CALL ZHPTRD( UPLO, N, AP, RWORK( INDD ), RWORK( INDE ), 00389 $ WORK( INDTAU ), IINFO ) 00390 * 00391 * If all eigenvalues are desired and ABSTOL is less than or equal 00392 * to zero, then call DSTERF or ZUPGTR and ZSTEQR. If this fails 00393 * for some eigenvalue, then try DSTEBZ. 00394 * 00395 TEST = .FALSE. 00396 IF (INDEIG) THEN 00397 IF (IL.EQ.1 .AND. IU.EQ.N) THEN 00398 TEST = .TRUE. 00399 END IF 00400 END IF 00401 IF ((ALLEIG .OR. TEST) .AND. (ABSTOL.LE.ZERO)) THEN 00402 CALL DCOPY( N, RWORK( INDD ), 1, W, 1 ) 00403 INDEE = INDRWK + 2*N 00404 IF( .NOT.WANTZ ) THEN 00405 CALL DCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00406 CALL DSTERF( N, W, RWORK( INDEE ), INFO ) 00407 ELSE 00408 CALL ZUPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ, 00409 $ WORK( INDWRK ), IINFO ) 00410 CALL DCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 ) 00411 CALL ZSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ, 00412 $ RWORK( INDRWK ), INFO ) 00413 IF( INFO.EQ.0 ) THEN 00414 DO 10 I = 1, N 00415 IFAIL( I ) = 0 00416 10 CONTINUE 00417 END IF 00418 END IF 00419 IF( INFO.EQ.0 ) THEN 00420 M = N 00421 GO TO 20 00422 END IF 00423 INFO = 0 00424 END IF 00425 * 00426 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN. 00427 * 00428 IF( WANTZ ) THEN 00429 ORDER = 'B' 00430 ELSE 00431 ORDER = 'E' 00432 END IF 00433 INDIBL = 1 00434 INDISP = INDIBL + N 00435 INDIWK = INDISP + N 00436 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 00437 $ RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W, 00438 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00439 $ IWORK( INDIWK ), INFO ) 00440 * 00441 IF( WANTZ ) THEN 00442 CALL ZSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W, 00443 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00444 $ RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO ) 00445 * 00446 * Apply unitary matrix used in reduction to tridiagonal 00447 * form to eigenvectors returned by ZSTEIN. 00448 * 00449 INDWRK = INDTAU + N 00450 CALL ZUPMTR( 'L', UPLO, 'N', N, M, AP, WORK( INDTAU ), Z, LDZ, 00451 $ WORK( INDWRK ), IINFO ) 00452 END IF 00453 * 00454 * If matrix was scaled, then rescale eigenvalues appropriately. 00455 * 00456 20 CONTINUE 00457 IF( ISCALE.EQ.1 ) THEN 00458 IF( INFO.EQ.0 ) THEN 00459 IMAX = M 00460 ELSE 00461 IMAX = INFO - 1 00462 END IF 00463 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00464 END IF 00465 * 00466 * If eigenvalues are not in order, then sort them, along with 00467 * eigenvectors. 00468 * 00469 IF( WANTZ ) THEN 00470 DO 40 J = 1, M - 1 00471 I = 0 00472 TMP1 = W( J ) 00473 DO 30 JJ = J + 1, M 00474 IF( W( JJ ).LT.TMP1 ) THEN 00475 I = JJ 00476 TMP1 = W( JJ ) 00477 END IF 00478 30 CONTINUE 00479 * 00480 IF( I.NE.0 ) THEN 00481 ITMP1 = IWORK( INDIBL+I-1 ) 00482 W( I ) = W( J ) 00483 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00484 W( J ) = TMP1 00485 IWORK( INDIBL+J-1 ) = ITMP1 00486 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00487 IF( INFO.NE.0 ) THEN 00488 ITMP1 = IFAIL( I ) 00489 IFAIL( I ) = IFAIL( J ) 00490 IFAIL( J ) = ITMP1 00491 END IF 00492 END IF 00493 40 CONTINUE 00494 END IF 00495 * 00496 RETURN 00497 * 00498 * End of ZHPEVX 00499 * 00500 END