 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZHEEVR( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, 00002 $ ABSTOL, M, W, Z, LDZ, ISUPPZ, WORK, LWORK, 00003 $ RWORK, LRWORK, IWORK, LIWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.3.1) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * -- April 2011 -- 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBZ, RANGE, UPLO 00012 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LRWORK, LWORK, 00013 $ M, N 00014 DOUBLE PRECISION ABSTOL, VL, VU 00015 * .. 00016 * .. Array Arguments .. 00017 INTEGER ISUPPZ( * ), IWORK( * ) 00018 DOUBLE PRECISION RWORK( * ), W( * ) 00019 COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * ZHEEVR computes selected eigenvalues and, optionally, eigenvectors 00026 * of a complex Hermitian matrix A. Eigenvalues and eigenvectors can 00027 * be selected by specifying either a range of values or a range of 00028 * indices for the desired eigenvalues. 00029 * 00030 * ZHEEVR first reduces the matrix A to tridiagonal form T with a call 00031 * to ZHETRD. Then, whenever possible, ZHEEVR calls ZSTEMR to compute 00032 * eigenspectrum using Relatively Robust Representations. ZSTEMR 00033 * computes eigenvalues by the dqds algorithm, while orthogonal 00034 * eigenvectors are computed from various "good" L D L^T representations 00035 * (also known as Relatively Robust Representations). Gram-Schmidt 00036 * orthogonalization is avoided as far as possible. More specifically, 00037 * the various steps of the algorithm are as follows. 00038 * 00039 * For each unreduced block (submatrix) of T, 00040 * (a) Compute T - sigma I = L D L^T, so that L and D 00041 * define all the wanted eigenvalues to high relative accuracy. 00042 * This means that small relative changes in the entries of D and L 00043 * cause only small relative changes in the eigenvalues and 00044 * eigenvectors. The standard (unfactored) representation of the 00045 * tridiagonal matrix T does not have this property in general. 00046 * (b) Compute the eigenvalues to suitable accuracy. 00047 * If the eigenvectors are desired, the algorithm attains full 00048 * accuracy of the computed eigenvalues only right before 00049 * the corresponding vectors have to be computed, see steps c) and d). 00050 * (c) For each cluster of close eigenvalues, select a new 00051 * shift close to the cluster, find a new factorization, and refine 00052 * the shifted eigenvalues to suitable accuracy. 00053 * (d) For each eigenvalue with a large enough relative separation compute 00054 * the corresponding eigenvector by forming a rank revealing twisted 00055 * factorization. Go back to (c) for any clusters that remain. 00056 * 00057 * The desired accuracy of the output can be specified by the input 00058 * parameter ABSTOL. 00059 * 00060 * For more details, see DSTEMR's documentation and: 00061 * - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations 00062 * to compute orthogonal eigenvectors of symmetric tridiagonal matrices," 00063 * Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. 00064 * - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and 00065 * Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25, 00066 * 2004. Also LAPACK Working Note 154. 00067 * - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric 00068 * tridiagonal eigenvalue/eigenvector problem", 00069 * Computer Science Division Technical Report No. UCB/CSD-97-971, 00070 * UC Berkeley, May 1997. 00071 * 00072 * 00073 * Note 1 : ZHEEVR calls ZSTEMR when the full spectrum is requested 00074 * on machines which conform to the ieee-754 floating point standard. 00075 * ZHEEVR calls DSTEBZ and ZSTEIN on non-ieee machines and 00076 * when partial spectrum requests are made. 00077 * 00078 * Normal execution of ZSTEMR may create NaNs and infinities and 00079 * hence may abort due to a floating point exception in environments 00080 * which do not handle NaNs and infinities in the ieee standard default 00081 * manner. 00082 * 00083 * Arguments 00084 * ========= 00085 * 00086 * JOBZ (input) CHARACTER*1 00087 * = 'N': Compute eigenvalues only; 00088 * = 'V': Compute eigenvalues and eigenvectors. 00089 * 00090 * RANGE (input) CHARACTER*1 00091 * = 'A': all eigenvalues will be found. 00092 * = 'V': all eigenvalues in the half-open interval (VL,VU] 00093 * will be found. 00094 * = 'I': the IL-th through IU-th eigenvalues will be found. 00095 * For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and 00096 * ZSTEIN are called 00097 * 00098 * UPLO (input) CHARACTER*1 00099 * = 'U': Upper triangle of A is stored; 00100 * = 'L': Lower triangle of A is stored. 00101 * 00102 * N (input) INTEGER 00103 * The order of the matrix A. N >= 0. 00104 * 00105 * A (input/output) COMPLEX*16 array, dimension (LDA, N) 00106 * On entry, the Hermitian matrix A. If UPLO = 'U', the 00107 * leading N-by-N upper triangular part of A contains the 00108 * upper triangular part of the matrix A. If UPLO = 'L', 00109 * the leading N-by-N lower triangular part of A contains 00110 * the lower triangular part of the matrix A. 00111 * On exit, the lower triangle (if UPLO='L') or the upper 00112 * triangle (if UPLO='U') of A, including the diagonal, is 00113 * destroyed. 00114 * 00115 * LDA (input) INTEGER 00116 * The leading dimension of the array A. LDA >= max(1,N). 00117 * 00118 * VL (input) DOUBLE PRECISION 00119 * VU (input) DOUBLE PRECISION 00120 * If RANGE='V', the lower and upper bounds of the interval to 00121 * be searched for eigenvalues. VL < VU. 00122 * Not referenced if RANGE = 'A' or 'I'. 00123 * 00124 * IL (input) INTEGER 00125 * IU (input) INTEGER 00126 * If RANGE='I', the indices (in ascending order) of the 00127 * smallest and largest eigenvalues to be returned. 00128 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00129 * Not referenced if RANGE = 'A' or 'V'. 00130 * 00131 * ABSTOL (input) DOUBLE PRECISION 00132 * The absolute error tolerance for the eigenvalues. 00133 * An approximate eigenvalue is accepted as converged 00134 * when it is determined to lie in an interval [a,b] 00135 * of width less than or equal to 00136 * 00137 * ABSTOL + EPS * max( |a|,|b| ) , 00138 * 00139 * where EPS is the machine precision. If ABSTOL is less than 00140 * or equal to zero, then EPS*|T| will be used in its place, 00141 * where |T| is the 1-norm of the tridiagonal matrix obtained 00142 * by reducing A to tridiagonal form. 00143 * 00144 * See "Computing Small Singular Values of Bidiagonal Matrices 00145 * with Guaranteed High Relative Accuracy," by Demmel and 00146 * Kahan, LAPACK Working Note #3. 00147 * 00148 * If high relative accuracy is important, set ABSTOL to 00149 * DLAMCH( 'Safe minimum' ). Doing so will guarantee that 00150 * eigenvalues are computed to high relative accuracy when 00151 * possible in future releases. The current code does not 00152 * make any guarantees about high relative accuracy, but 00153 * furutre releases will. See J. Barlow and J. Demmel, 00154 * "Computing Accurate Eigensystems of Scaled Diagonally 00155 * Dominant Matrices", LAPACK Working Note #7, for a discussion 00156 * of which matrices define their eigenvalues to high relative 00157 * accuracy. 00158 * 00159 * M (output) INTEGER 00160 * The total number of eigenvalues found. 0 <= M <= N. 00161 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00162 * 00163 * W (output) DOUBLE PRECISION array, dimension (N) 00164 * The first M elements contain the selected eigenvalues in 00165 * ascending order. 00166 * 00167 * Z (output) COMPLEX*16 array, dimension (LDZ, max(1,M)) 00168 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00169 * contain the orthonormal eigenvectors of the matrix A 00170 * corresponding to the selected eigenvalues, with the i-th 00171 * column of Z holding the eigenvector associated with W(i). 00172 * If JOBZ = 'N', then Z is not referenced. 00173 * Note: the user must ensure that at least max(1,M) columns are 00174 * supplied in the array Z; if RANGE = 'V', the exact value of M 00175 * is not known in advance and an upper bound must be used. 00176 * 00177 * LDZ (input) INTEGER 00178 * The leading dimension of the array Z. LDZ >= 1, and if 00179 * JOBZ = 'V', LDZ >= max(1,N). 00180 * 00181 * ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) 00182 * The support of the eigenvectors in Z, i.e., the indices 00183 * indicating the nonzero elements in Z. The i-th eigenvector 00184 * is nonzero only in elements ISUPPZ( 2*i-1 ) through 00185 * ISUPPZ( 2*i ). 00186 * Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 00187 * 00188 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00189 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00190 * 00191 * LWORK (input) INTEGER 00192 * The length of the array WORK. LWORK >= max(1,2*N). 00193 * For optimal efficiency, LWORK >= (NB+1)*N, 00194 * where NB is the max of the blocksize for ZHETRD and for 00195 * ZUNMTR as returned by ILAENV. 00196 * 00197 * If LWORK = -1, then a workspace query is assumed; the routine 00198 * only calculates the optimal sizes of the WORK, RWORK and 00199 * IWORK arrays, returns these values as the first entries of 00200 * the WORK, RWORK and IWORK arrays, and no error message 00201 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00202 * 00203 * RWORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LRWORK)) 00204 * On exit, if INFO = 0, RWORK(1) returns the optimal 00205 * (and minimal) LRWORK. 00206 * 00207 * LRWORK (input) INTEGER 00208 * The length of the array RWORK. LRWORK >= max(1,24*N). 00209 * 00210 * If LRWORK = -1, then a workspace query is assumed; the 00211 * routine only calculates the optimal sizes of the WORK, RWORK 00212 * and IWORK arrays, returns these values as the first entries 00213 * of the WORK, RWORK and IWORK arrays, and no error message 00214 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00215 * 00216 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00217 * On exit, if INFO = 0, IWORK(1) returns the optimal 00218 * (and minimal) LIWORK. 00219 * 00220 * LIWORK (input) INTEGER 00221 * The dimension of the array IWORK. LIWORK >= max(1,10*N). 00222 * 00223 * If LIWORK = -1, then a workspace query is assumed; the 00224 * routine only calculates the optimal sizes of the WORK, RWORK 00225 * and IWORK arrays, returns these values as the first entries 00226 * of the WORK, RWORK and IWORK arrays, and no error message 00227 * related to LWORK or LRWORK or LIWORK is issued by XERBLA. 00228 * 00229 * INFO (output) INTEGER 00230 * = 0: successful exit 00231 * < 0: if INFO = -i, the i-th argument had an illegal value 00232 * > 0: Internal error 00233 * 00234 * Further Details 00235 * =============== 00236 * 00237 * Based on contributions by 00238 * Inderjit Dhillon, IBM Almaden, USA 00239 * Osni Marques, LBNL/NERSC, USA 00240 * Ken Stanley, Computer Science Division, University of 00241 * California at Berkeley, USA 00242 * Jason Riedy, Computer Science Division, University of 00243 * California at Berkeley, USA 00244 * 00245 * ===================================================================== 00246 * 00247 * .. Parameters .. 00248 DOUBLE PRECISION ZERO, ONE, TWO 00249 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00250 * .. 00251 * .. Local Scalars .. 00252 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG, 00253 $ WANTZ, TRYRAC 00254 CHARACTER ORDER 00255 INTEGER I, IEEEOK, IINFO, IMAX, INDIBL, INDIFL, INDISP, 00256 $ INDIWO, INDRD, INDRDD, INDRE, INDREE, INDRWK, 00257 $ INDTAU, INDWK, INDWKN, ISCALE, ITMP1, J, JJ, 00258 $ LIWMIN, LLWORK, LLRWORK, LLWRKN, LRWMIN, 00259 $ LWKOPT, LWMIN, NB, NSPLIT 00260 DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, 00261 $ SIGMA, SMLNUM, TMP1, VLL, VUU 00262 * .. 00263 * .. External Functions .. 00264 LOGICAL LSAME 00265 INTEGER ILAENV 00266 DOUBLE PRECISION DLAMCH, ZLANSY 00267 EXTERNAL LSAME, ILAENV, DLAMCH, ZLANSY 00268 * .. 00269 * .. External Subroutines .. 00270 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTERF, XERBLA, ZDSCAL, 00271 $ ZHETRD, ZSTEMR, ZSTEIN, ZSWAP, ZUNMTR 00272 * .. 00273 * .. Intrinsic Functions .. 00274 INTRINSIC DBLE, MAX, MIN, SQRT 00275 * .. 00276 * .. Executable Statements .. 00277 * 00278 * Test the input parameters. 00279 * 00280 IEEEOK = ILAENV( 10, 'ZHEEVR', 'N', 1, 2, 3, 4 ) 00281 * 00282 LOWER = LSAME( UPLO, 'L' ) 00283 WANTZ = LSAME( JOBZ, 'V' ) 00284 ALLEIG = LSAME( RANGE, 'A' ) 00285 VALEIG = LSAME( RANGE, 'V' ) 00286 INDEIG = LSAME( RANGE, 'I' ) 00287 * 00288 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LRWORK.EQ.-1 ) .OR. 00289 $ ( LIWORK.EQ.-1 ) ) 00290 * 00291 LRWMIN = MAX( 1, 24*N ) 00292 LIWMIN = MAX( 1, 10*N ) 00293 LWMIN = MAX( 1, 2*N ) 00294 * 00295 INFO = 0 00296 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00297 INFO = -1 00298 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00299 INFO = -2 00300 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN 00301 INFO = -3 00302 ELSE IF( N.LT.0 ) THEN 00303 INFO = -4 00304 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00305 INFO = -6 00306 ELSE 00307 IF( VALEIG ) THEN 00308 IF( N.GT.0 .AND. VU.LE.VL ) 00309 $ INFO = -8 00310 ELSE IF( INDEIG ) THEN 00311 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00312 INFO = -9 00313 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00314 INFO = -10 00315 END IF 00316 END IF 00317 END IF 00318 IF( INFO.EQ.0 ) THEN 00319 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00320 INFO = -15 00321 END IF 00322 END IF 00323 * 00324 IF( INFO.EQ.0 ) THEN 00325 NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) 00326 NB = MAX( NB, ILAENV( 1, 'ZUNMTR', UPLO, N, -1, -1, -1 ) ) 00327 LWKOPT = MAX( ( NB+1 )*N, LWMIN ) 00328 WORK( 1 ) = LWKOPT 00329 RWORK( 1 ) = LRWMIN 00330 IWORK( 1 ) = LIWMIN 00331 * 00332 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00333 INFO = -18 00334 ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN 00335 INFO = -20 00336 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00337 INFO = -22 00338 END IF 00339 END IF 00340 * 00341 IF( INFO.NE.0 ) THEN 00342 CALL XERBLA( 'ZHEEVR', -INFO ) 00343 RETURN 00344 ELSE IF( LQUERY ) THEN 00345 RETURN 00346 END IF 00347 * 00348 * Quick return if possible 00349 * 00350 M = 0 00351 IF( N.EQ.0 ) THEN 00352 WORK( 1 ) = 1 00353 RETURN 00354 END IF 00355 * 00356 IF( N.EQ.1 ) THEN 00357 WORK( 1 ) = 2 00358 IF( ALLEIG .OR. INDEIG ) THEN 00359 M = 1 00360 W( 1 ) = DBLE( A( 1, 1 ) ) 00361 ELSE 00362 IF( VL.LT.DBLE( A( 1, 1 ) ) .AND. VU.GE.DBLE( A( 1, 1 ) ) ) 00363 $ THEN 00364 M = 1 00365 W( 1 ) = DBLE( A( 1, 1 ) ) 00366 END IF 00367 END IF 00368 IF( WANTZ ) THEN 00369 Z( 1, 1 ) = ONE 00370 ISUPPZ( 1 ) = 1 00371 ISUPPZ( 2 ) = 1 00372 END IF 00373 RETURN 00374 END IF 00375 * 00376 * Get machine constants. 00377 * 00378 SAFMIN = DLAMCH( 'Safe minimum' ) 00379 EPS = DLAMCH( 'Precision' ) 00380 SMLNUM = SAFMIN / EPS 00381 BIGNUM = ONE / SMLNUM 00382 RMIN = SQRT( SMLNUM ) 00383 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00384 * 00385 * Scale matrix to allowable range, if necessary. 00386 * 00387 ISCALE = 0 00388 ABSTLL = ABSTOL 00389 IF (VALEIG) THEN 00390 VLL = VL 00391 VUU = VU 00392 END IF 00393 ANRM = ZLANSY( 'M', UPLO, N, A, LDA, RWORK ) 00394 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN 00395 ISCALE = 1 00396 SIGMA = RMIN / ANRM 00397 ELSE IF( ANRM.GT.RMAX ) THEN 00398 ISCALE = 1 00399 SIGMA = RMAX / ANRM 00400 END IF 00401 IF( ISCALE.EQ.1 ) THEN 00402 IF( LOWER ) THEN 00403 DO 10 J = 1, N 00404 CALL ZDSCAL( N-J+1, SIGMA, A( J, J ), 1 ) 00405 10 CONTINUE 00406 ELSE 00407 DO 20 J = 1, N 00408 CALL ZDSCAL( J, SIGMA, A( 1, J ), 1 ) 00409 20 CONTINUE 00410 END IF 00411 IF( ABSTOL.GT.0 ) 00412 $ ABSTLL = ABSTOL*SIGMA 00413 IF( VALEIG ) THEN 00414 VLL = VL*SIGMA 00415 VUU = VU*SIGMA 00416 END IF 00417 END IF 00418 00419 * Initialize indices into workspaces. Note: The IWORK indices are 00420 * used only if DSTERF or ZSTEMR fail. 00421 00422 * WORK(INDTAU:INDTAU+N-1) stores the complex scalar factors of the 00423 * elementary reflectors used in ZHETRD. 00424 INDTAU = 1 00425 * INDWK is the starting offset of the remaining complex workspace, 00426 * and LLWORK is the remaining complex workspace size. 00427 INDWK = INDTAU + N 00428 LLWORK = LWORK - INDWK + 1 00429 00430 * RWORK(INDRD:INDRD+N-1) stores the real tridiagonal's diagonal 00431 * entries. 00432 INDRD = 1 00433 * RWORK(INDRE:INDRE+N-1) stores the off-diagonal entries of the 00434 * tridiagonal matrix from ZHETRD. 00435 INDRE = INDRD + N 00436 * RWORK(INDRDD:INDRDD+N-1) is a copy of the diagonal entries over 00437 * -written by ZSTEMR (the DSTERF path copies the diagonal to W). 00438 INDRDD = INDRE + N 00439 * RWORK(INDREE:INDREE+N-1) is a copy of the off-diagonal entries over 00440 * -written while computing the eigenvalues in DSTERF and ZSTEMR. 00441 INDREE = INDRDD + N 00442 * INDRWK is the starting offset of the left-over real workspace, and 00443 * LLRWORK is the remaining workspace size. 00444 INDRWK = INDREE + N 00445 LLRWORK = LRWORK - INDRWK + 1 00446 00447 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and 00448 * stores the block indices of each of the M<=N eigenvalues. 00449 INDIBL = 1 00450 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and 00451 * stores the starting and finishing indices of each block. 00452 INDISP = INDIBL + N 00453 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors 00454 * that corresponding to eigenvectors that fail to converge in 00455 * DSTEIN. This information is discarded; if any fail, the driver 00456 * returns INFO > 0. 00457 INDIFL = INDISP + N 00458 * INDIWO is the offset of the remaining integer workspace. 00459 INDIWO = INDISP + N 00460 00461 * 00462 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form. 00463 * 00464 CALL ZHETRD( UPLO, N, A, LDA, RWORK( INDRD ), RWORK( INDRE ), 00465 $ WORK( INDTAU ), WORK( INDWK ), LLWORK, IINFO ) 00466 * 00467 * If all eigenvalues are desired 00468 * then call DSTERF or ZSTEMR and ZUNMTR. 00469 * 00470 TEST = .FALSE. 00471 IF( INDEIG ) THEN 00472 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00473 TEST = .TRUE. 00474 END IF 00475 END IF 00476 IF( ( ALLEIG.OR.TEST ) .AND. ( IEEEOK.EQ.1 ) ) THEN 00477 IF( .NOT.WANTZ ) THEN 00478 CALL DCOPY( N, RWORK( INDRD ), 1, W, 1 ) 00479 CALL DCOPY( N-1, RWORK( INDRE ), 1, RWORK( INDREE ), 1 ) 00480 CALL DSTERF( N, W, RWORK( INDREE ), INFO ) 00481 ELSE 00482 CALL DCOPY( N-1, RWORK( INDRE ), 1, RWORK( INDREE ), 1 ) 00483 CALL DCOPY( N, RWORK( INDRD ), 1, RWORK( INDRDD ), 1 ) 00484 * 00485 IF (ABSTOL .LE. TWO*N*EPS) THEN 00486 TRYRAC = .TRUE. 00487 ELSE 00488 TRYRAC = .FALSE. 00489 END IF 00490 CALL ZSTEMR( JOBZ, 'A', N, RWORK( INDRDD ), 00491 $ RWORK( INDREE ), VL, VU, IL, IU, M, W, 00492 $ Z, LDZ, N, ISUPPZ, TRYRAC, 00493 $ RWORK( INDRWK ), LLRWORK, 00494 $ IWORK, LIWORK, INFO ) 00495 * 00496 * Apply unitary matrix used in reduction to tridiagonal 00497 * form to eigenvectors returned by ZSTEIN. 00498 * 00499 IF( WANTZ .AND. INFO.EQ.0 ) THEN 00500 INDWKN = INDWK 00501 LLWRKN = LWORK - INDWKN + 1 00502 CALL ZUNMTR( 'L', UPLO, 'N', N, M, A, LDA, 00503 $ WORK( INDTAU ), Z, LDZ, WORK( INDWKN ), 00504 $ LLWRKN, IINFO ) 00505 END IF 00506 END IF 00507 * 00508 * 00509 IF( INFO.EQ.0 ) THEN 00510 M = N 00511 GO TO 30 00512 END IF 00513 INFO = 0 00514 END IF 00515 * 00516 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN. 00517 * Also call DSTEBZ and ZSTEIN if ZSTEMR fails. 00518 * 00519 IF( WANTZ ) THEN 00520 ORDER = 'B' 00521 ELSE 00522 ORDER = 'E' 00523 END IF 00524 00525 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, 00526 $ RWORK( INDRD ), RWORK( INDRE ), M, NSPLIT, W, 00527 $ IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ), 00528 $ IWORK( INDIWO ), INFO ) 00529 * 00530 IF( WANTZ ) THEN 00531 CALL ZSTEIN( N, RWORK( INDRD ), RWORK( INDRE ), M, W, 00532 $ IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ, 00533 $ RWORK( INDRWK ), IWORK( INDIWO ), IWORK( INDIFL ), 00534 $ INFO ) 00535 * 00536 * Apply unitary matrix used in reduction to tridiagonal 00537 * form to eigenvectors returned by ZSTEIN. 00538 * 00539 INDWKN = INDWK 00540 LLWRKN = LWORK - INDWKN + 1 00541 CALL ZUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z, 00542 $ LDZ, WORK( INDWKN ), LLWRKN, IINFO ) 00543 END IF 00544 * 00545 * If matrix was scaled, then rescale eigenvalues appropriately. 00546 * 00547 30 CONTINUE 00548 IF( ISCALE.EQ.1 ) THEN 00549 IF( INFO.EQ.0 ) THEN 00550 IMAX = M 00551 ELSE 00552 IMAX = INFO - 1 00553 END IF 00554 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00555 END IF 00556 * 00557 * If eigenvalues are not in order, then sort them, along with 00558 * eigenvectors. 00559 * 00560 IF( WANTZ ) THEN 00561 DO 50 J = 1, M - 1 00562 I = 0 00563 TMP1 = W( J ) 00564 DO 40 JJ = J + 1, M 00565 IF( W( JJ ).LT.TMP1 ) THEN 00566 I = JJ 00567 TMP1 = W( JJ ) 00568 END IF 00569 40 CONTINUE 00570 * 00571 IF( I.NE.0 ) THEN 00572 ITMP1 = IWORK( INDIBL+I-1 ) 00573 W( I ) = W( J ) 00574 IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 ) 00575 W( J ) = TMP1 00576 IWORK( INDIBL+J-1 ) = ITMP1 00577 CALL ZSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00578 END IF 00579 50 CONTINUE 00580 END IF 00581 * 00582 * Set WORK(1) to optimal workspace size. 00583 * 00584 WORK( 1 ) = LWKOPT 00585 RWORK( 1 ) = LRWMIN 00586 IWORK( 1 ) = LIWMIN 00587 * 00588 RETURN 00589 * 00590 * End of ZHEEVR 00591 * 00592 END
1.7.3 
