00001 SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, 00002 $ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, 00003 $ LIWORK, INFO ) 00004 * 00005 * -- LAPACK driver routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER JOBZ, RANGE 00012 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N 00013 DOUBLE PRECISION ABSTOL, VL, VU 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER ISUPPZ( * ), IWORK( * ) 00017 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DSTEVR computes selected eigenvalues and, optionally, eigenvectors 00024 * of a real symmetric tridiagonal matrix T. Eigenvalues and 00025 * eigenvectors can be selected by specifying either a range of values 00026 * or a range of indices for the desired eigenvalues. 00027 * 00028 * Whenever possible, DSTEVR calls DSTEMR to compute the 00029 * eigenspectrum using Relatively Robust Representations. DSTEMR 00030 * computes eigenvalues by the dqds algorithm, while orthogonal 00031 * eigenvectors are computed from various "good" L D L^T representations 00032 * (also known as Relatively Robust Representations). Gram-Schmidt 00033 * orthogonalization is avoided as far as possible. More specifically, 00034 * the various steps of the algorithm are as follows. For the i-th 00035 * unreduced block of T, 00036 * (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T 00037 * is a relatively robust representation, 00038 * (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high 00039 * relative accuracy by the dqds algorithm, 00040 * (c) If there is a cluster of close eigenvalues, "choose" sigma_i 00041 * close to the cluster, and go to step (a), 00042 * (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, 00043 * compute the corresponding eigenvector by forming a 00044 * rank-revealing twisted factorization. 00045 * The desired accuracy of the output can be specified by the input 00046 * parameter ABSTOL. 00047 * 00048 * For more details, see "A new O(n^2) algorithm for the symmetric 00049 * tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, 00050 * Computer Science Division Technical Report No. UCB//CSD-97-971, 00051 * UC Berkeley, May 1997. 00052 * 00053 * 00054 * Note 1 : DSTEVR calls DSTEMR when the full spectrum is requested 00055 * on machines which conform to the ieee-754 floating point standard. 00056 * DSTEVR calls DSTEBZ and DSTEIN on non-ieee machines and 00057 * when partial spectrum requests are made. 00058 * 00059 * Normal execution of DSTEMR may create NaNs and infinities and 00060 * hence may abort due to a floating point exception in environments 00061 * which do not handle NaNs and infinities in the ieee standard default 00062 * manner. 00063 * 00064 * Arguments 00065 * ========= 00066 * 00067 * JOBZ (input) CHARACTER*1 00068 * = 'N': Compute eigenvalues only; 00069 * = 'V': Compute eigenvalues and eigenvectors. 00070 * 00071 * RANGE (input) CHARACTER*1 00072 * = 'A': all eigenvalues will be found. 00073 * = 'V': all eigenvalues in the half-open interval (VL,VU] 00074 * will be found. 00075 * = 'I': the IL-th through IU-th eigenvalues will be found. 00076 ********** For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and 00077 ********** DSTEIN are called 00078 * 00079 * N (input) INTEGER 00080 * The order of the matrix. N >= 0. 00081 * 00082 * D (input/output) DOUBLE PRECISION array, dimension (N) 00083 * On entry, the n diagonal elements of the tridiagonal matrix 00084 * A. 00085 * On exit, D may be multiplied by a constant factor chosen 00086 * to avoid over/underflow in computing the eigenvalues. 00087 * 00088 * E (input/output) DOUBLE PRECISION array, dimension (max(1,N-1)) 00089 * On entry, the (n-1) subdiagonal elements of the tridiagonal 00090 * matrix A in elements 1 to N-1 of E. 00091 * On exit, E may be multiplied by a constant factor chosen 00092 * to avoid over/underflow in computing the eigenvalues. 00093 * 00094 * VL (input) DOUBLE PRECISION 00095 * VU (input) DOUBLE PRECISION 00096 * If RANGE='V', the lower and upper bounds of the interval to 00097 * be searched for eigenvalues. VL < VU. 00098 * Not referenced if RANGE = 'A' or 'I'. 00099 * 00100 * IL (input) INTEGER 00101 * IU (input) INTEGER 00102 * If RANGE='I', the indices (in ascending order) of the 00103 * smallest and largest eigenvalues to be returned. 00104 * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 00105 * Not referenced if RANGE = 'A' or 'V'. 00106 * 00107 * ABSTOL (input) DOUBLE PRECISION 00108 * The absolute error tolerance for the eigenvalues. 00109 * An approximate eigenvalue is accepted as converged 00110 * when it is determined to lie in an interval [a,b] 00111 * of width less than or equal to 00112 * 00113 * ABSTOL + EPS * max( |a|,|b| ) , 00114 * 00115 * where EPS is the machine precision. If ABSTOL is less than 00116 * or equal to zero, then EPS*|T| will be used in its place, 00117 * where |T| is the 1-norm of the tridiagonal matrix obtained 00118 * by reducing A to tridiagonal form. 00119 * 00120 * See "Computing Small Singular Values of Bidiagonal Matrices 00121 * with Guaranteed High Relative Accuracy," by Demmel and 00122 * Kahan, LAPACK Working Note #3. 00123 * 00124 * If high relative accuracy is important, set ABSTOL to 00125 * DLAMCH( 'Safe minimum' ). Doing so will guarantee that 00126 * eigenvalues are computed to high relative accuracy when 00127 * possible in future releases. The current code does not 00128 * make any guarantees about high relative accuracy, but 00129 * future releases will. See J. Barlow and J. Demmel, 00130 * "Computing Accurate Eigensystems of Scaled Diagonally 00131 * Dominant Matrices", LAPACK Working Note #7, for a discussion 00132 * of which matrices define their eigenvalues to high relative 00133 * accuracy. 00134 * 00135 * M (output) INTEGER 00136 * The total number of eigenvalues found. 0 <= M <= N. 00137 * If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 00138 * 00139 * W (output) DOUBLE PRECISION array, dimension (N) 00140 * The first M elements contain the selected eigenvalues in 00141 * ascending order. 00142 * 00143 * Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) 00144 * If JOBZ = 'V', then if INFO = 0, the first M columns of Z 00145 * contain the orthonormal eigenvectors of the matrix A 00146 * corresponding to the selected eigenvalues, with the i-th 00147 * column of Z holding the eigenvector associated with W(i). 00148 * Note: the user must ensure that at least max(1,M) columns are 00149 * supplied in the array Z; if RANGE = 'V', the exact value of M 00150 * is not known in advance and an upper bound must be used. 00151 * 00152 * LDZ (input) INTEGER 00153 * The leading dimension of the array Z. LDZ >= 1, and if 00154 * JOBZ = 'V', LDZ >= max(1,N). 00155 * 00156 * ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) 00157 * The support of the eigenvectors in Z, i.e., the indices 00158 * indicating the nonzero elements in Z. The i-th eigenvector 00159 * is nonzero only in elements ISUPPZ( 2*i-1 ) through 00160 * ISUPPZ( 2*i ). 00161 ********** Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 00162 * 00163 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00164 * On exit, if INFO = 0, WORK(1) returns the optimal (and 00165 * minimal) LWORK. 00166 * 00167 * LWORK (input) INTEGER 00168 * The dimension of the array WORK. LWORK >= max(1,20*N). 00169 * 00170 * If LWORK = -1, then a workspace query is assumed; the routine 00171 * only calculates the optimal sizes of the WORK and IWORK 00172 * arrays, returns these values as the first entries of the WORK 00173 * and IWORK arrays, and no error message related to LWORK or 00174 * LIWORK is issued by XERBLA. 00175 * 00176 * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK)) 00177 * On exit, if INFO = 0, IWORK(1) returns the optimal (and 00178 * minimal) LIWORK. 00179 * 00180 * LIWORK (input) INTEGER 00181 * The dimension of the array IWORK. LIWORK >= max(1,10*N). 00182 * 00183 * If LIWORK = -1, then a workspace query is assumed; the 00184 * routine only calculates the optimal sizes of the WORK and 00185 * IWORK arrays, returns these values as the first entries of 00186 * the WORK and IWORK arrays, and no error message related to 00187 * LWORK or LIWORK is issued by XERBLA. 00188 * 00189 * INFO (output) INTEGER 00190 * = 0: successful exit 00191 * < 0: if INFO = -i, the i-th argument had an illegal value 00192 * > 0: Internal error 00193 * 00194 * Further Details 00195 * =============== 00196 * 00197 * Based on contributions by 00198 * Inderjit Dhillon, IBM Almaden, USA 00199 * Osni Marques, LBNL/NERSC, USA 00200 * Ken Stanley, Computer Science Division, University of 00201 * California at Berkeley, USA 00202 * 00203 * ===================================================================== 00204 * 00205 * .. Parameters .. 00206 DOUBLE PRECISION ZERO, ONE, TWO 00207 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00208 * .. 00209 * .. Local Scalars .. 00210 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ, 00211 $ TRYRAC 00212 CHARACTER ORDER 00213 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP, 00214 $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN, 00215 $ NSPLIT 00216 DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, 00217 $ TMP1, TNRM, VLL, VUU 00218 * .. 00219 * .. External Functions .. 00220 LOGICAL LSAME 00221 INTEGER ILAENV 00222 DOUBLE PRECISION DLAMCH, DLANST 00223 EXTERNAL LSAME, ILAENV, DLAMCH, DLANST 00224 * .. 00225 * .. External Subroutines .. 00226 EXTERNAL DCOPY, DSCAL, DSTEBZ, DSTEMR, DSTEIN, DSTERF, 00227 $ DSWAP, XERBLA 00228 * .. 00229 * .. Intrinsic Functions .. 00230 INTRINSIC MAX, MIN, SQRT 00231 * .. 00232 * .. Executable Statements .. 00233 * 00234 * 00235 * Test the input parameters. 00236 * 00237 IEEEOK = ILAENV( 10, 'DSTEVR', 'N', 1, 2, 3, 4 ) 00238 * 00239 WANTZ = LSAME( JOBZ, 'V' ) 00240 ALLEIG = LSAME( RANGE, 'A' ) 00241 VALEIG = LSAME( RANGE, 'V' ) 00242 INDEIG = LSAME( RANGE, 'I' ) 00243 * 00244 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) ) 00245 LWMIN = MAX( 1, 20*N ) 00246 LIWMIN = MAX( 1, 10*N ) 00247 * 00248 * 00249 INFO = 0 00250 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00251 INFO = -1 00252 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00253 INFO = -2 00254 ELSE IF( N.LT.0 ) THEN 00255 INFO = -3 00256 ELSE 00257 IF( VALEIG ) THEN 00258 IF( N.GT.0 .AND. VU.LE.VL ) 00259 $ INFO = -7 00260 ELSE IF( INDEIG ) THEN 00261 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00262 INFO = -8 00263 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00264 INFO = -9 00265 END IF 00266 END IF 00267 END IF 00268 IF( INFO.EQ.0 ) THEN 00269 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00270 INFO = -14 00271 END IF 00272 END IF 00273 * 00274 IF( INFO.EQ.0 ) THEN 00275 WORK( 1 ) = LWMIN 00276 IWORK( 1 ) = LIWMIN 00277 * 00278 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00279 INFO = -17 00280 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00281 INFO = -19 00282 END IF 00283 END IF 00284 * 00285 IF( INFO.NE.0 ) THEN 00286 CALL XERBLA( 'DSTEVR', -INFO ) 00287 RETURN 00288 ELSE IF( LQUERY ) THEN 00289 RETURN 00290 END IF 00291 * 00292 * Quick return if possible 00293 * 00294 M = 0 00295 IF( N.EQ.0 ) 00296 $ RETURN 00297 * 00298 IF( N.EQ.1 ) THEN 00299 IF( ALLEIG .OR. INDEIG ) THEN 00300 M = 1 00301 W( 1 ) = D( 1 ) 00302 ELSE 00303 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN 00304 M = 1 00305 W( 1 ) = D( 1 ) 00306 END IF 00307 END IF 00308 IF( WANTZ ) 00309 $ Z( 1, 1 ) = ONE 00310 RETURN 00311 END IF 00312 * 00313 * Get machine constants. 00314 * 00315 SAFMIN = DLAMCH( 'Safe minimum' ) 00316 EPS = DLAMCH( 'Precision' ) 00317 SMLNUM = SAFMIN / EPS 00318 BIGNUM = ONE / SMLNUM 00319 RMIN = SQRT( SMLNUM ) 00320 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00321 * 00322 * 00323 * Scale matrix to allowable range, if necessary. 00324 * 00325 ISCALE = 0 00326 VLL = VL 00327 VUU = VU 00328 * 00329 TNRM = DLANST( 'M', N, D, E ) 00330 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN 00331 ISCALE = 1 00332 SIGMA = RMIN / TNRM 00333 ELSE IF( TNRM.GT.RMAX ) THEN 00334 ISCALE = 1 00335 SIGMA = RMAX / TNRM 00336 END IF 00337 IF( ISCALE.EQ.1 ) THEN 00338 CALL DSCAL( N, SIGMA, D, 1 ) 00339 CALL DSCAL( N-1, SIGMA, E( 1 ), 1 ) 00340 IF( VALEIG ) THEN 00341 VLL = VL*SIGMA 00342 VUU = VU*SIGMA 00343 END IF 00344 END IF 00345 00346 * Initialize indices into workspaces. Note: These indices are used only 00347 * if DSTERF or DSTEMR fail. 00348 00349 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in DSTEBZ and 00350 * stores the block indices of each of the M<=N eigenvalues. 00351 INDIBL = 1 00352 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in DSTEBZ and 00353 * stores the starting and finishing indices of each block. 00354 INDISP = INDIBL + N 00355 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors 00356 * that corresponding to eigenvectors that fail to converge in 00357 * DSTEIN. This information is discarded; if any fail, the driver 00358 * returns INFO > 0. 00359 INDIFL = INDISP + N 00360 * INDIWO is the offset of the remaining integer workspace. 00361 INDIWO = INDISP + N 00362 * 00363 * If all eigenvalues are desired, then 00364 * call DSTERF or DSTEMR. If this fails for some eigenvalue, then 00365 * try DSTEBZ. 00366 * 00367 * 00368 TEST = .FALSE. 00369 IF( INDEIG ) THEN 00370 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00371 TEST = .TRUE. 00372 END IF 00373 END IF 00374 IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN 00375 CALL DCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 ) 00376 IF( .NOT.WANTZ ) THEN 00377 CALL DCOPY( N, D, 1, W, 1 ) 00378 CALL DSTERF( N, W, WORK, INFO ) 00379 ELSE 00380 CALL DCOPY( N, D, 1, WORK( N+1 ), 1 ) 00381 IF (ABSTOL .LE. TWO*N*EPS) THEN 00382 TRYRAC = .TRUE. 00383 ELSE 00384 TRYRAC = .FALSE. 00385 END IF 00386 CALL DSTEMR( JOBZ, 'A', N, WORK( N+1 ), WORK, VL, VU, IL, 00387 $ IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC, 00388 $ WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO ) 00389 * 00390 END IF 00391 IF( INFO.EQ.0 ) THEN 00392 M = N 00393 GO TO 10 00394 END IF 00395 INFO = 0 00396 END IF 00397 * 00398 * Otherwise, call DSTEBZ and, if eigenvectors are desired, DSTEIN. 00399 * 00400 IF( WANTZ ) THEN 00401 ORDER = 'B' 00402 ELSE 00403 ORDER = 'E' 00404 END IF 00405 00406 CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M, 00407 $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK, 00408 $ IWORK( INDIWO ), INFO ) 00409 * 00410 IF( WANTZ ) THEN 00411 CALL DSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ), 00412 $ Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ), 00413 $ INFO ) 00414 END IF 00415 * 00416 * If matrix was scaled, then rescale eigenvalues appropriately. 00417 * 00418 10 CONTINUE 00419 IF( ISCALE.EQ.1 ) THEN 00420 IF( INFO.EQ.0 ) THEN 00421 IMAX = M 00422 ELSE 00423 IMAX = INFO - 1 00424 END IF 00425 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 ) 00426 END IF 00427 * 00428 * If eigenvalues are not in order, then sort them, along with 00429 * eigenvectors. 00430 * 00431 IF( WANTZ ) THEN 00432 DO 30 J = 1, M - 1 00433 I = 0 00434 TMP1 = W( J ) 00435 DO 20 JJ = J + 1, M 00436 IF( W( JJ ).LT.TMP1 ) THEN 00437 I = JJ 00438 TMP1 = W( JJ ) 00439 END IF 00440 20 CONTINUE 00441 * 00442 IF( I.NE.0 ) THEN 00443 ITMP1 = IWORK( I ) 00444 W( I ) = W( J ) 00445 IWORK( I ) = IWORK( J ) 00446 W( J ) = TMP1 00447 IWORK( J ) = ITMP1 00448 CALL DSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00449 END IF 00450 30 CONTINUE 00451 END IF 00452 * 00453 * Causes problems with tests 19 & 20: 00454 * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002 00455 * 00456 * 00457 WORK( 1 ) = LWMIN 00458 IWORK( 1 ) = LIWMIN 00459 RETURN 00460 * 00461 * End of DSTEVR 00462 * 00463 END
1.7.2