![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZSTEMR 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZSTEMR + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, 00022 * M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, 00023 * IWORK, LIWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER JOBZ, RANGE 00027 * LOGICAL TRYRAC 00028 * INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N 00029 * DOUBLE PRECISION VL, VU 00030 * .. 00031 * .. Array Arguments .. 00032 * INTEGER ISUPPZ( * ), IWORK( * ) 00033 * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 00034 * COMPLEX*16 Z( LDZ, * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> ZSTEMR computes selected eigenvalues and, optionally, eigenvectors 00044 *> of a real symmetric tridiagonal matrix T. Any such unreduced matrix has 00045 *> a well defined set of pairwise different real eigenvalues, the corresponding 00046 *> real eigenvectors are pairwise orthogonal. 00047 *> 00048 *> The spectrum may be computed either completely or partially by specifying 00049 *> either an interval (VL,VU] or a range of indices IL:IU for the desired 00050 *> eigenvalues. 00051 *> 00052 *> Depending on the number of desired eigenvalues, these are computed either 00053 *> by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are 00054 *> computed by the use of various suitable L D L^T factorizations near clusters 00055 *> of close eigenvalues (referred to as RRRs, Relatively Robust 00056 *> Representations). An informal sketch of the algorithm follows. 00057 *> 00058 *> For each unreduced block (submatrix) of T, 00059 *> (a) Compute T - sigma I = L D L^T, so that L and D 00060 *> define all the wanted eigenvalues to high relative accuracy. 00061 *> This means that small relative changes in the entries of D and L 00062 *> cause only small relative changes in the eigenvalues and 00063 *> eigenvectors. The standard (unfactored) representation of the 00064 *> tridiagonal matrix T does not have this property in general. 00065 *> (b) Compute the eigenvalues to suitable accuracy. 00066 *> If the eigenvectors are desired, the algorithm attains full 00067 *> accuracy of the computed eigenvalues only right before 00068 *> the corresponding vectors have to be computed, see steps c) and d). 00069 *> (c) For each cluster of close eigenvalues, select a new 00070 *> shift close to the cluster, find a new factorization, and refine 00071 *> the shifted eigenvalues to suitable accuracy. 00072 *> (d) For each eigenvalue with a large enough relative separation compute 00073 *> the corresponding eigenvector by forming a rank revealing twisted 00074 *> factorization. Go back to (c) for any clusters that remain. 00075 *> 00076 *> For more details, see: 00077 *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations 00078 *> to compute orthogonal eigenvectors of symmetric tridiagonal matrices," 00079 *> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. 00080 *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and 00081 *> Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, 00082 *> 2004. Also LAPACK Working Note 154. 00083 *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric 00084 *> tridiagonal eigenvalue/eigenvector problem", 00085 *> Computer Science Division Technical Report No. UCB/CSD-97-971, 00086 *> UC Berkeley, May 1997. 00087 *> 00088 *> Further Details 00089 *> 1.ZSTEMR works only on machines which follow IEEE-754 00090 *> floating-point standard in their handling of infinities and NaNs. 00091 *> This permits the use of efficient inner loops avoiding a check for 00092 *> zero divisors. 00093 *> 00094 *> 2. LAPACK routines can be used to reduce a complex Hermitean matrix to 00095 *> real symmetric tridiagonal form. 00096 *> 00097 *> (Any complex Hermitean tridiagonal matrix has real values on its diagonal 00098 *> and potentially complex numbers on its off-diagonals. By applying a 00099 *> similarity transform with an appropriate diagonal matrix 00100 *> diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean 00101 *> matrix can be transformed into a real symmetric matrix and complex 00102 *> arithmetic can be entirely avoided.) 00103 *> 00104 *> While the eigenvectors of the real symmetric tridiagonal matrix are real, 00105 *> the eigenvectors of original complex Hermitean matrix have complex entries 00106 *> in general. 00107 *> Since LAPACK drivers overwrite the matrix data with the eigenvectors, 00108 *> ZSTEMR accepts complex workspace to facilitate interoperability 00109 *> with ZUNMTR or ZUPMTR. 00110 *> \endverbatim 00111 * 00112 * Arguments: 00113 * ========== 00114 * 00115 *> \param[in] JOBZ 00116 *> \verbatim 00117 *> JOBZ is CHARACTER*1 00118 *> = 'N': Compute eigenvalues only; 00119 *> = 'V': Compute eigenvalues and eigenvectors. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] RANGE 00123 *> \verbatim 00124 *> RANGE is CHARACTER*1 00125 *> = 'A': all eigenvalues will be found. 00126 *> = 'V': all eigenvalues in the half-open interval (VL,VU] 00127 *> will be found. 00128 *> = 'I': the IL-th through IU-th eigenvalues will be found. 00129 *> \endverbatim 00130 *> 00131 *> \param[in] N 00132 *> \verbatim 00133 *> N is INTEGER 00134 *> The order of the matrix. N >= 0. 00135 *> \endverbatim 00136 *> 00137 *> \param[in,out] D 00138 *> \verbatim 00139 *> D is DOUBLE PRECISION array, dimension (N) 00140 *> On entry, the N diagonal elements of the tridiagonal matrix 00141 *> T. On exit, D is overwritten. 00142 *> \endverbatim 00143 *> 00144 *> \param[in,out] E 00145 *> \verbatim 00146 *> E is DOUBLE PRECISION array, dimension (N) 00147 *> On entry, the (N-1) subdiagonal elements of the tridiagonal 00148 *> matrix T in elements 1 to N-1 of E. E(N) need not be set on 00149 *> input, but is used internally as workspace. 00150 *> On exit, E is overwritten. 00151 *> \endverbatim 00152 *> 00153 *> \param[in] VL 00154 *> \verbatim 00155 *> VL is DOUBLE PRECISION 00156 *> \endverbatim 00157 *> 00158 *> \param[in] VU 00159 *> \verbatim 00160 *> VU is DOUBLE PRECISION 00161 *> 00162 *> If RANGE='V', the lower and upper bounds of the interval to 00163 *> be searched for eigenvalues. VL < VU. 00164 *> Not referenced if RANGE = 'A' or 'I'. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] IL 00168 *> \verbatim 00169 *> IL is INTEGER 00170 *> \endverbatim 00171 *> 00172 *> \param[in] IU 00173 *> \verbatim 00174 *> IU is INTEGER 00175 *> 00176 *> If RANGE='I', the indices (in ascending order) of the 00177 *> smallest and largest eigenvalues to be returned. 00178 *> 1 <= IL <= IU <= N, if N > 0. 00179 *> Not referenced if RANGE = 'A' or 'V'. 00180 *> \endverbatim 00181 *> 00182 *> \param[out] M 00183 *> \verbatim 00184 *> M is INTEGER 00185 *> The total number of eigenvalues found. 0 <= M <= N. 00186 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00187 *> \endverbatim 00188 *> 00189 *> \param[out] W 00190 *> \verbatim 00191 *> W is DOUBLE PRECISION array, dimension (N) 00192 *> The first M elements contain the selected eigenvalues in 00193 *> ascending order. 00194 *> \endverbatim 00195 *> 00196 *> \param[out] Z 00197 *> \verbatim 00198 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M) ) 00199 *> If JOBZ = 'V', and if INFO = 0, then the first M columns of Z 00200 *> contain the orthonormal eigenvectors of the matrix T 00201 *> corresponding to the selected eigenvalues, with the i-th 00202 *> column of Z holding the eigenvector associated with W(i). 00203 *> If JOBZ = 'N', then Z is not referenced. 00204 *> Note: the user must ensure that at least max(1,M) columns are 00205 *> supplied in the array Z; if RANGE = 'V', the exact value of M 00206 *> is not known in advance and can be computed with a workspace 00207 *> query by setting NZC = -1, see below. 00208 *> \endverbatim 00209 *> 00210 *> \param[in] LDZ 00211 *> \verbatim 00212 *> LDZ is INTEGER 00213 *> The leading dimension of the array Z. LDZ >= 1, and if 00214 *> JOBZ = 'V', then LDZ >= max(1,N). 00215 *> \endverbatim 00216 *> 00217 *> \param[in] NZC 00218 *> \verbatim 00219 *> NZC is INTEGER 00220 *> The number of eigenvectors to be held in the array Z. 00221 *> If RANGE = 'A', then NZC >= max(1,N). 00222 *> If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU]. 00223 *> If RANGE = 'I', then NZC >= IU-IL+1. 00224 *> If NZC = -1, then a workspace query is assumed; the 00225 *> routine calculates the number of columns of the array Z that 00226 *> are needed to hold the eigenvectors. 00227 *> This value is returned as the first entry of the Z array, and 00228 *> no error message related to NZC is issued by XERBLA. 00229 *> \endverbatim 00230 *> 00231 *> \param[out] ISUPPZ 00232 *> \verbatim 00233 *> ISUPPZ is INTEGER ARRAY, dimension ( 2*max(1,M) ) 00234 *> The support of the eigenvectors in Z, i.e., the indices 00235 *> indicating the nonzero elements in Z. The i-th computed eigenvector 00236 *> is nonzero only in elements ISUPPZ( 2*i-1 ) through 00237 *> ISUPPZ( 2*i ). This is relevant in the case when the matrix 00238 *> is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0. 00239 *> \endverbatim 00240 *> 00241 *> \param[in,out] TRYRAC 00242 *> \verbatim 00243 *> TRYRAC is LOGICAL 00244 *> If TRYRAC.EQ..TRUE., indicates that the code should check whether 00245 *> the tridiagonal matrix defines its eigenvalues to high relative 00246 *> accuracy. If so, the code uses relative-accuracy preserving 00247 *> algorithms that might be (a bit) slower depending on the matrix. 00248 *> If the matrix does not define its eigenvalues to high relative 00249 *> accuracy, the code can uses possibly faster algorithms. 00250 *> If TRYRAC.EQ..FALSE., the code is not required to guarantee 00251 *> relatively accurate eigenvalues and can use the fastest possible 00252 *> techniques. 00253 *> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix 00254 *> does not define its eigenvalues to high relative accuracy. 00255 *> \endverbatim 00256 *> 00257 *> \param[out] WORK 00258 *> \verbatim 00259 *> WORK is DOUBLE PRECISION array, dimension (LWORK) 00260 *> On exit, if INFO = 0, WORK(1) returns the optimal 00261 *> (and minimal) LWORK. 00262 *> \endverbatim 00263 *> 00264 *> \param[in] LWORK 00265 *> \verbatim 00266 *> LWORK is INTEGER 00267 *> The dimension of the array WORK. LWORK >= max(1,18*N) 00268 *> if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'. 00269 *> If LWORK = -1, then a workspace query is assumed; the routine 00270 *> only calculates the optimal size of the WORK array, returns 00271 *> this value as the first entry of the WORK array, and no error 00272 *> message related to LWORK is issued by XERBLA. 00273 *> \endverbatim 00274 *> 00275 *> \param[out] IWORK 00276 *> \verbatim 00277 *> IWORK is INTEGER array, dimension (LIWORK) 00278 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 00279 *> \endverbatim 00280 *> 00281 *> \param[in] LIWORK 00282 *> \verbatim 00283 *> LIWORK is INTEGER 00284 *> The dimension of the array IWORK. LIWORK >= max(1,10*N) 00285 *> if the eigenvectors are desired, and LIWORK >= max(1,8*N) 00286 *> if only the eigenvalues are to be computed. 00287 *> If LIWORK = -1, then a workspace query is assumed; the 00288 *> routine only calculates the optimal size of the IWORK array, 00289 *> returns this value as the first entry of the IWORK array, and 00290 *> no error message related to LIWORK is issued by XERBLA. 00291 *> \endverbatim 00292 *> 00293 *> \param[out] INFO 00294 *> \verbatim 00295 *> INFO is INTEGER 00296 *> On exit, INFO 00297 *> = 0: successful exit 00298 *> < 0: if INFO = -i, the i-th argument had an illegal value 00299 *> > 0: if INFO = 1X, internal error in DLARRE, 00300 *> if INFO = 2X, internal error in ZLARRV. 00301 *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is 00302 *> the nonzero error code returned by DLARRE or 00303 *> ZLARRV, respectively. 00304 *> \endverbatim 00305 * 00306 * Authors: 00307 * ======== 00308 * 00309 *> \author Univ. of Tennessee 00310 *> \author Univ. of California Berkeley 00311 *> \author Univ. of Colorado Denver 00312 *> \author NAG Ltd. 00313 * 00314 *> \date November 2011 00315 * 00316 *> \ingroup complex16OTHERcomputational 00317 * 00318 *> \par Contributors: 00319 * ================== 00320 *> 00321 *> Beresford Parlett, University of California, Berkeley, USA \n 00322 *> Jim Demmel, University of California, Berkeley, USA \n 00323 *> Inderjit Dhillon, University of Texas, Austin, USA \n 00324 *> Osni Marques, LBNL/NERSC, USA \n 00325 *> Christof Voemel, University of California, Berkeley, USA \n 00326 * 00327 * ===================================================================== 00328 SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, 00329 $ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, 00330 $ IWORK, LIWORK, INFO ) 00331 * 00332 * -- LAPACK computational routine (version 3.4.0) -- 00333 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00334 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00335 * November 2011 00336 * 00337 * .. Scalar Arguments .. 00338 CHARACTER JOBZ, RANGE 00339 LOGICAL TRYRAC 00340 INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N 00341 DOUBLE PRECISION VL, VU 00342 * .. 00343 * .. Array Arguments .. 00344 INTEGER ISUPPZ( * ), IWORK( * ) 00345 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ) 00346 COMPLEX*16 Z( LDZ, * ) 00347 * .. 00348 * 00349 * ===================================================================== 00350 * 00351 * .. Parameters .. 00352 DOUBLE PRECISION ZERO, ONE, FOUR, MINRGP 00353 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, 00354 $ FOUR = 4.0D0, 00355 $ MINRGP = 1.0D-3 ) 00356 * .. 00357 * .. Local Scalars .. 00358 LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY 00359 INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW, 00360 $ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD, 00361 $ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP, 00362 $ ITMP2, J, JBLK, JJ, LIWMIN, LWMIN, NSPLIT, 00363 $ NZCMIN, OFFSET, WBEGIN, WEND 00364 DOUBLE PRECISION BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN, 00365 $ RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN, 00366 $ THRESH, TMP, TNRM, WL, WU 00367 * .. 00368 * .. 00369 * .. External Functions .. 00370 LOGICAL LSAME 00371 DOUBLE PRECISION DLAMCH, DLANST 00372 EXTERNAL LSAME, DLAMCH, DLANST 00373 * .. 00374 * .. External Subroutines .. 00375 EXTERNAL DCOPY, DLAE2, DLAEV2, DLARRC, DLARRE, DLARRJ, 00376 $ DLARRR, DLASRT, DSCAL, XERBLA, ZLARRV, ZSWAP 00377 * .. 00378 * .. Intrinsic Functions .. 00379 INTRINSIC MAX, MIN, SQRT 00380 00381 00382 * .. 00383 * .. Executable Statements .. 00384 * 00385 * Test the input parameters. 00386 * 00387 WANTZ = LSAME( JOBZ, 'V' ) 00388 ALLEIG = LSAME( RANGE, 'A' ) 00389 VALEIG = LSAME( RANGE, 'V' ) 00390 INDEIG = LSAME( RANGE, 'I' ) 00391 * 00392 LQUERY = ( ( LWORK.EQ.-1 ).OR.( LIWORK.EQ.-1 ) ) 00393 ZQUERY = ( NZC.EQ.-1 ) 00394 00395 * DSTEMR needs WORK of size 6*N, IWORK of size 3*N. 00396 * In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N. 00397 * Furthermore, ZLARRV needs WORK of size 12*N, IWORK of size 7*N. 00398 IF( WANTZ ) THEN 00399 LWMIN = 18*N 00400 LIWMIN = 10*N 00401 ELSE 00402 * need less workspace if only the eigenvalues are wanted 00403 LWMIN = 12*N 00404 LIWMIN = 8*N 00405 ENDIF 00406 00407 WL = ZERO 00408 WU = ZERO 00409 IIL = 0 00410 IIU = 0 00411 00412 IF( VALEIG ) THEN 00413 * We do not reference VL, VU in the cases RANGE = 'I','A' 00414 * The interval (WL, WU] contains all the wanted eigenvalues. 00415 * It is either given by the user or computed in DLARRE. 00416 WL = VL 00417 WU = VU 00418 ELSEIF( INDEIG ) THEN 00419 * We do not reference IL, IU in the cases RANGE = 'V','A' 00420 IIL = IL 00421 IIU = IU 00422 ENDIF 00423 * 00424 INFO = 0 00425 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00426 INFO = -1 00427 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00428 INFO = -2 00429 ELSE IF( N.LT.0 ) THEN 00430 INFO = -3 00431 ELSE IF( VALEIG .AND. N.GT.0 .AND. WU.LE.WL ) THEN 00432 INFO = -7 00433 ELSE IF( INDEIG .AND. ( IIL.LT.1 .OR. IIL.GT.N ) ) THEN 00434 INFO = -8 00435 ELSE IF( INDEIG .AND. ( IIU.LT.IIL .OR. IIU.GT.N ) ) THEN 00436 INFO = -9 00437 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00438 INFO = -13 00439 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00440 INFO = -17 00441 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00442 INFO = -19 00443 END IF 00444 * 00445 * Get machine constants. 00446 * 00447 SAFMIN = DLAMCH( 'Safe minimum' ) 00448 EPS = DLAMCH( 'Precision' ) 00449 SMLNUM = SAFMIN / EPS 00450 BIGNUM = ONE / SMLNUM 00451 RMIN = SQRT( SMLNUM ) 00452 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00453 * 00454 IF( INFO.EQ.0 ) THEN 00455 WORK( 1 ) = LWMIN 00456 IWORK( 1 ) = LIWMIN 00457 * 00458 IF( WANTZ .AND. ALLEIG ) THEN 00459 NZCMIN = N 00460 ELSE IF( WANTZ .AND. VALEIG ) THEN 00461 CALL DLARRC( 'T', N, VL, VU, D, E, SAFMIN, 00462 $ NZCMIN, ITMP, ITMP2, INFO ) 00463 ELSE IF( WANTZ .AND. INDEIG ) THEN 00464 NZCMIN = IIU-IIL+1 00465 ELSE 00466 * WANTZ .EQ. FALSE. 00467 NZCMIN = 0 00468 ENDIF 00469 IF( ZQUERY .AND. INFO.EQ.0 ) THEN 00470 Z( 1,1 ) = NZCMIN 00471 ELSE IF( NZC.LT.NZCMIN .AND. .NOT.ZQUERY ) THEN 00472 INFO = -14 00473 END IF 00474 END IF 00475 00476 IF( INFO.NE.0 ) THEN 00477 * 00478 CALL XERBLA( 'ZSTEMR', -INFO ) 00479 * 00480 RETURN 00481 ELSE IF( LQUERY .OR. ZQUERY ) THEN 00482 RETURN 00483 END IF 00484 * 00485 * Handle N = 0, 1, and 2 cases immediately 00486 * 00487 M = 0 00488 IF( N.EQ.0 ) 00489 $ RETURN 00490 * 00491 IF( N.EQ.1 ) THEN 00492 IF( ALLEIG .OR. INDEIG ) THEN 00493 M = 1 00494 W( 1 ) = D( 1 ) 00495 ELSE 00496 IF( WL.LT.D( 1 ) .AND. WU.GE.D( 1 ) ) THEN 00497 M = 1 00498 W( 1 ) = D( 1 ) 00499 END IF 00500 END IF 00501 IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN 00502 Z( 1, 1 ) = ONE 00503 ISUPPZ(1) = 1 00504 ISUPPZ(2) = 1 00505 END IF 00506 RETURN 00507 END IF 00508 * 00509 IF( N.EQ.2 ) THEN 00510 IF( .NOT.WANTZ ) THEN 00511 CALL DLAE2( D(1), E(1), D(2), R1, R2 ) 00512 ELSE IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN 00513 CALL DLAEV2( D(1), E(1), D(2), R1, R2, CS, SN ) 00514 END IF 00515 IF( ALLEIG.OR. 00516 $ (VALEIG.AND.(R2.GT.WL).AND. 00517 $ (R2.LE.WU)).OR. 00518 $ (INDEIG.AND.(IIL.EQ.1)) ) THEN 00519 M = M+1 00520 W( M ) = R2 00521 IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN 00522 Z( 1, M ) = -SN 00523 Z( 2, M ) = CS 00524 * Note: At most one of SN and CS can be zero. 00525 IF (SN.NE.ZERO) THEN 00526 IF (CS.NE.ZERO) THEN 00527 ISUPPZ(2*M-1) = 1 00528 ISUPPZ(2*M-1) = 2 00529 ELSE 00530 ISUPPZ(2*M-1) = 1 00531 ISUPPZ(2*M-1) = 1 00532 END IF 00533 ELSE 00534 ISUPPZ(2*M-1) = 2 00535 ISUPPZ(2*M) = 2 00536 END IF 00537 ENDIF 00538 ENDIF 00539 IF( ALLEIG.OR. 00540 $ (VALEIG.AND.(R1.GT.WL).AND. 00541 $ (R1.LE.WU)).OR. 00542 $ (INDEIG.AND.(IIU.EQ.2)) ) THEN 00543 M = M+1 00544 W( M ) = R1 00545 IF( WANTZ.AND.(.NOT.ZQUERY) ) THEN 00546 Z( 1, M ) = CS 00547 Z( 2, M ) = SN 00548 * Note: At most one of SN and CS can be zero. 00549 IF (SN.NE.ZERO) THEN 00550 IF (CS.NE.ZERO) THEN 00551 ISUPPZ(2*M-1) = 1 00552 ISUPPZ(2*M-1) = 2 00553 ELSE 00554 ISUPPZ(2*M-1) = 1 00555 ISUPPZ(2*M-1) = 1 00556 END IF 00557 ELSE 00558 ISUPPZ(2*M-1) = 2 00559 ISUPPZ(2*M) = 2 00560 END IF 00561 ENDIF 00562 ENDIF 00563 RETURN 00564 END IF 00565 00566 * Continue with general N 00567 00568 INDGRS = 1 00569 INDERR = 2*N + 1 00570 INDGP = 3*N + 1 00571 INDD = 4*N + 1 00572 INDE2 = 5*N + 1 00573 INDWRK = 6*N + 1 00574 * 00575 IINSPL = 1 00576 IINDBL = N + 1 00577 IINDW = 2*N + 1 00578 IINDWK = 3*N + 1 00579 * 00580 * Scale matrix to allowable range, if necessary. 00581 * The allowable range is related to the PIVMIN parameter; see the 00582 * comments in DLARRD. The preference for scaling small values 00583 * up is heuristic; we expect users' matrices not to be close to the 00584 * RMAX threshold. 00585 * 00586 SCALE = ONE 00587 TNRM = DLANST( 'M', N, D, E ) 00588 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN 00589 SCALE = RMIN / TNRM 00590 ELSE IF( TNRM.GT.RMAX ) THEN 00591 SCALE = RMAX / TNRM 00592 END IF 00593 IF( SCALE.NE.ONE ) THEN 00594 CALL DSCAL( N, SCALE, D, 1 ) 00595 CALL DSCAL( N-1, SCALE, E, 1 ) 00596 TNRM = TNRM*SCALE 00597 IF( VALEIG ) THEN 00598 * If eigenvalues in interval have to be found, 00599 * scale (WL, WU] accordingly 00600 WL = WL*SCALE 00601 WU = WU*SCALE 00602 ENDIF 00603 END IF 00604 * 00605 * Compute the desired eigenvalues of the tridiagonal after splitting 00606 * into smaller subblocks if the corresponding off-diagonal elements 00607 * are small 00608 * THRESH is the splitting parameter for DLARRE 00609 * A negative THRESH forces the old splitting criterion based on the 00610 * size of the off-diagonal. A positive THRESH switches to splitting 00611 * which preserves relative accuracy. 00612 * 00613 IF( TRYRAC ) THEN 00614 * Test whether the matrix warrants the more expensive relative approach. 00615 CALL DLARRR( N, D, E, IINFO ) 00616 ELSE 00617 * The user does not care about relative accurately eigenvalues 00618 IINFO = -1 00619 ENDIF 00620 * Set the splitting criterion 00621 IF (IINFO.EQ.0) THEN 00622 THRESH = EPS 00623 ELSE 00624 THRESH = -EPS 00625 * relative accuracy is desired but T does not guarantee it 00626 TRYRAC = .FALSE. 00627 ENDIF 00628 * 00629 IF( TRYRAC ) THEN 00630 * Copy original diagonal, needed to guarantee relative accuracy 00631 CALL DCOPY(N,D,1,WORK(INDD),1) 00632 ENDIF 00633 * Store the squares of the offdiagonal values of T 00634 DO 5 J = 1, N-1 00635 WORK( INDE2+J-1 ) = E(J)**2 00636 5 CONTINUE 00637 00638 * Set the tolerance parameters for bisection 00639 IF( .NOT.WANTZ ) THEN 00640 * DLARRE computes the eigenvalues to full precision. 00641 RTOL1 = FOUR * EPS 00642 RTOL2 = FOUR * EPS 00643 ELSE 00644 * DLARRE computes the eigenvalues to less than full precision. 00645 * ZLARRV will refine the eigenvalue approximations, and we only 00646 * need less accurate initial bisection in DLARRE. 00647 * Note: these settings do only affect the subset case and DLARRE 00648 RTOL1 = SQRT(EPS) 00649 RTOL2 = MAX( SQRT(EPS)*5.0D-3, FOUR * EPS ) 00650 ENDIF 00651 CALL DLARRE( RANGE, N, WL, WU, IIL, IIU, D, E, 00652 $ WORK(INDE2), RTOL1, RTOL2, THRESH, NSPLIT, 00653 $ IWORK( IINSPL ), M, W, WORK( INDERR ), 00654 $ WORK( INDGP ), IWORK( IINDBL ), 00655 $ IWORK( IINDW ), WORK( INDGRS ), PIVMIN, 00656 $ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 00657 IF( IINFO.NE.0 ) THEN 00658 INFO = 10 + ABS( IINFO ) 00659 RETURN 00660 END IF 00661 * Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired 00662 * part of the spectrum. All desired eigenvalues are contained in 00663 * (WL,WU] 00664 00665 00666 IF( WANTZ ) THEN 00667 * 00668 * Compute the desired eigenvectors corresponding to the computed 00669 * eigenvalues 00670 * 00671 CALL ZLARRV( N, WL, WU, D, E, 00672 $ PIVMIN, IWORK( IINSPL ), M, 00673 $ 1, M, MINRGP, RTOL1, RTOL2, 00674 $ W, WORK( INDERR ), WORK( INDGP ), IWORK( IINDBL ), 00675 $ IWORK( IINDW ), WORK( INDGRS ), Z, LDZ, 00676 $ ISUPPZ, WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 00677 IF( IINFO.NE.0 ) THEN 00678 INFO = 20 + ABS( IINFO ) 00679 RETURN 00680 END IF 00681 ELSE 00682 * DLARRE computes eigenvalues of the (shifted) root representation 00683 * ZLARRV returns the eigenvalues of the unshifted matrix. 00684 * However, if the eigenvectors are not desired by the user, we need 00685 * to apply the corresponding shifts from DLARRE to obtain the 00686 * eigenvalues of the original matrix. 00687 DO 20 J = 1, M 00688 ITMP = IWORK( IINDBL+J-1 ) 00689 W( J ) = W( J ) + E( IWORK( IINSPL+ITMP-1 ) ) 00690 20 CONTINUE 00691 END IF 00692 * 00693 00694 IF ( TRYRAC ) THEN 00695 * Refine computed eigenvalues so that they are relatively accurate 00696 * with respect to the original matrix T. 00697 IBEGIN = 1 00698 WBEGIN = 1 00699 DO 39 JBLK = 1, IWORK( IINDBL+M-1 ) 00700 IEND = IWORK( IINSPL+JBLK-1 ) 00701 IN = IEND - IBEGIN + 1 00702 WEND = WBEGIN - 1 00703 * check if any eigenvalues have to be refined in this block 00704 36 CONTINUE 00705 IF( WEND.LT.M ) THEN 00706 IF( IWORK( IINDBL+WEND ).EQ.JBLK ) THEN 00707 WEND = WEND + 1 00708 GO TO 36 00709 END IF 00710 END IF 00711 IF( WEND.LT.WBEGIN ) THEN 00712 IBEGIN = IEND + 1 00713 GO TO 39 00714 END IF 00715 00716 OFFSET = IWORK(IINDW+WBEGIN-1)-1 00717 IFIRST = IWORK(IINDW+WBEGIN-1) 00718 ILAST = IWORK(IINDW+WEND-1) 00719 RTOL2 = FOUR * EPS 00720 CALL DLARRJ( IN, 00721 $ WORK(INDD+IBEGIN-1), WORK(INDE2+IBEGIN-1), 00722 $ IFIRST, ILAST, RTOL2, OFFSET, W(WBEGIN), 00723 $ WORK( INDERR+WBEGIN-1 ), 00724 $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, 00725 $ TNRM, IINFO ) 00726 IBEGIN = IEND + 1 00727 WBEGIN = WEND + 1 00728 39 CONTINUE 00729 ENDIF 00730 * 00731 * If matrix was scaled, then rescale eigenvalues appropriately. 00732 * 00733 IF( SCALE.NE.ONE ) THEN 00734 CALL DSCAL( M, ONE / SCALE, W, 1 ) 00735 END IF 00736 * 00737 * If eigenvalues are not in increasing order, then sort them, 00738 * possibly along with eigenvectors. 00739 * 00740 IF( NSPLIT.GT.1 ) THEN 00741 IF( .NOT. WANTZ ) THEN 00742 CALL DLASRT( 'I', M, W, IINFO ) 00743 IF( IINFO.NE.0 ) THEN 00744 INFO = 3 00745 RETURN 00746 END IF 00747 ELSE 00748 DO 60 J = 1, M - 1 00749 I = 0 00750 TMP = W( J ) 00751 DO 50 JJ = J + 1, M 00752 IF( W( JJ ).LT.TMP ) THEN 00753 I = JJ 00754 TMP = W( JJ ) 00755 END IF 00756 50 CONTINUE 00757 IF( I.NE.0 ) THEN 00758 W( I ) = W( J ) 00759 W( J ) = TMP 00760 IF( WANTZ ) THEN 00761 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00762 ITMP = ISUPPZ( 2*I-1 ) 00763 ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 ) 00764 ISUPPZ( 2*J-1 ) = ITMP 00765 ITMP = ISUPPZ( 2*I ) 00766 ISUPPZ( 2*I ) = ISUPPZ( 2*J ) 00767 ISUPPZ( 2*J ) = ITMP 00768 END IF 00769 END IF 00770 60 CONTINUE 00771 END IF 00772 ENDIF 00773 * 00774 * 00775 WORK( 1 ) = LWMIN 00776 IWORK( 1 ) = LIWMIN 00777 RETURN 00778 * 00779 * End of ZSTEMR 00780 * 00781 END