 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SSTEVR( 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 REAL ABSTOL, VL, VU 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER ISUPPZ( * ), IWORK( * ) 00017 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * SSTEVR 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, SSTEVR calls SSTEMR to compute the 00029 * eigenspectrum using Relatively Robust Representations. SSTEMR 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 : SSTEVR calls SSTEMR when the full spectrum is requested 00055 * on machines which conform to the ieee-754 floating point standard. 00056 * SSTEVR calls SSTEBZ and SSTEIN on non-ieee machines and 00057 * when partial spectrum requests are made. 00058 * 00059 * Normal execution of SSTEMR 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, SSTEBZ and 00077 ********** SSTEIN are called 00078 * 00079 * N (input) INTEGER 00080 * The order of the matrix. N >= 0. 00081 * 00082 * D (input/output) REAL 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) REAL 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) REAL 00095 * VU (input) REAL 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) REAL 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 * SLAMCH( '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) REAL array, dimension (N) 00140 * The first M elements contain the selected eigenvalues in 00141 * ascending order. 00142 * 00143 * Z (output) REAL 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) REAL 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 >= 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 >= 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 * Jason Riedy, Computer Science Division, University of 00203 * California at Berkeley, USA 00204 * 00205 * ===================================================================== 00206 * 00207 * .. Parameters .. 00208 REAL ZERO, ONE, TWO 00209 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 00210 * .. 00211 * .. Local Scalars .. 00212 LOGICAL ALLEIG, INDEIG, TEST, LQUERY, VALEIG, WANTZ, 00213 $ TRYRAC 00214 CHARACTER ORDER 00215 INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP, 00216 $ INDIWO, ISCALE, J, JJ, LIWMIN, LWMIN, NSPLIT 00217 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, 00218 $ TMP1, TNRM, VLL, VUU 00219 * .. 00220 * .. External Functions .. 00221 LOGICAL LSAME 00222 INTEGER ILAENV 00223 REAL SLAMCH, SLANST 00224 EXTERNAL LSAME, ILAENV, SLAMCH, SLANST 00225 * .. 00226 * .. External Subroutines .. 00227 EXTERNAL SCOPY, SSCAL, SSTEBZ, SSTEMR, SSTEIN, SSTERF, 00228 $ SSWAP, XERBLA 00229 * .. 00230 * .. Intrinsic Functions .. 00231 INTRINSIC MAX, MIN, SQRT 00232 * .. 00233 * .. Executable Statements .. 00234 * 00235 * 00236 * Test the input parameters. 00237 * 00238 IEEEOK = ILAENV( 10, 'SSTEVR', 'N', 1, 2, 3, 4 ) 00239 * 00240 WANTZ = LSAME( JOBZ, 'V' ) 00241 ALLEIG = LSAME( RANGE, 'A' ) 00242 VALEIG = LSAME( RANGE, 'V' ) 00243 INDEIG = LSAME( RANGE, 'I' ) 00244 * 00245 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) ) 00246 LWMIN = MAX( 1, 20*N ) 00247 LIWMIN = MAX(1, 10*N ) 00248 * 00249 * 00250 INFO = 0 00251 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 00252 INFO = -1 00253 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 00254 INFO = -2 00255 ELSE IF( N.LT.0 ) THEN 00256 INFO = -3 00257 ELSE 00258 IF( VALEIG ) THEN 00259 IF( N.GT.0 .AND. VU.LE.VL ) 00260 $ INFO = -7 00261 ELSE IF( INDEIG ) THEN 00262 IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN 00263 INFO = -8 00264 ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN 00265 INFO = -9 00266 END IF 00267 END IF 00268 END IF 00269 IF( INFO.EQ.0 ) THEN 00270 IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 00271 INFO = -14 00272 END IF 00273 END IF 00274 * 00275 IF( INFO.EQ.0 ) THEN 00276 WORK( 1 ) = LWMIN 00277 IWORK( 1 ) = LIWMIN 00278 * 00279 IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 00280 INFO = -17 00281 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 00282 INFO = -19 00283 END IF 00284 END IF 00285 * 00286 IF( INFO.NE.0 ) THEN 00287 CALL XERBLA( 'SSTEVR', -INFO ) 00288 RETURN 00289 ELSE IF( LQUERY ) THEN 00290 RETURN 00291 END IF 00292 * 00293 * Quick return if possible 00294 * 00295 M = 0 00296 IF( N.EQ.0 ) 00297 $ RETURN 00298 * 00299 IF( N.EQ.1 ) THEN 00300 IF( ALLEIG .OR. INDEIG ) THEN 00301 M = 1 00302 W( 1 ) = D( 1 ) 00303 ELSE 00304 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN 00305 M = 1 00306 W( 1 ) = D( 1 ) 00307 END IF 00308 END IF 00309 IF( WANTZ ) 00310 $ Z( 1, 1 ) = ONE 00311 RETURN 00312 END IF 00313 * 00314 * Get machine constants. 00315 * 00316 SAFMIN = SLAMCH( 'Safe minimum' ) 00317 EPS = SLAMCH( 'Precision' ) 00318 SMLNUM = SAFMIN / EPS 00319 BIGNUM = ONE / SMLNUM 00320 RMIN = SQRT( SMLNUM ) 00321 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 00322 * 00323 * 00324 * Scale matrix to allowable range, if necessary. 00325 * 00326 ISCALE = 0 00327 VLL = VL 00328 VUU = VU 00329 * 00330 TNRM = SLANST( 'M', N, D, E ) 00331 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN 00332 ISCALE = 1 00333 SIGMA = RMIN / TNRM 00334 ELSE IF( TNRM.GT.RMAX ) THEN 00335 ISCALE = 1 00336 SIGMA = RMAX / TNRM 00337 END IF 00338 IF( ISCALE.EQ.1 ) THEN 00339 CALL SSCAL( N, SIGMA, D, 1 ) 00340 CALL SSCAL( N-1, SIGMA, E( 1 ), 1 ) 00341 IF( VALEIG ) THEN 00342 VLL = VL*SIGMA 00343 VUU = VU*SIGMA 00344 END IF 00345 END IF 00346 00347 * Initialize indices into workspaces. Note: These indices are used only 00348 * if SSTERF or SSTEMR fail. 00349 00350 * IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and 00351 * stores the block indices of each of the M<=N eigenvalues. 00352 INDIBL = 1 00353 * IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and 00354 * stores the starting and finishing indices of each block. 00355 INDISP = INDIBL + N 00356 * IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors 00357 * that corresponding to eigenvectors that fail to converge in 00358 * SSTEIN. This information is discarded; if any fail, the driver 00359 * returns INFO > 0. 00360 INDIFL = INDISP + N 00361 * INDIWO is the offset of the remaining integer workspace. 00362 INDIWO = INDISP + N 00363 * 00364 * If all eigenvalues are desired, then 00365 * call SSTERF or SSTEMR. If this fails for some eigenvalue, then 00366 * try SSTEBZ. 00367 * 00368 * 00369 TEST = .FALSE. 00370 IF( INDEIG ) THEN 00371 IF( IL.EQ.1 .AND. IU.EQ.N ) THEN 00372 TEST = .TRUE. 00373 END IF 00374 END IF 00375 IF( ( ALLEIG .OR. TEST ) .AND. IEEEOK.EQ.1 ) THEN 00376 CALL SCOPY( N-1, E( 1 ), 1, WORK( 1 ), 1 ) 00377 IF( .NOT.WANTZ ) THEN 00378 CALL SCOPY( N, D, 1, W, 1 ) 00379 CALL SSTERF( N, W, WORK, INFO ) 00380 ELSE 00381 CALL SCOPY( N, D, 1, WORK( N+1 ), 1 ) 00382 IF (ABSTOL .LE. TWO*N*EPS) THEN 00383 TRYRAC = .TRUE. 00384 ELSE 00385 TRYRAC = .FALSE. 00386 END IF 00387 CALL SSTEMR( JOBZ, 'A', N, WORK( N+1 ), WORK, VL, VU, IL, 00388 $ IU, M, W, Z, LDZ, N, ISUPPZ, TRYRAC, 00389 $ WORK( 2*N+1 ), LWORK-2*N, IWORK, LIWORK, INFO ) 00390 * 00391 END IF 00392 IF( INFO.EQ.0 ) THEN 00393 M = N 00394 GO TO 10 00395 END IF 00396 INFO = 0 00397 END IF 00398 * 00399 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN. 00400 * 00401 IF( WANTZ ) THEN 00402 ORDER = 'B' 00403 ELSE 00404 ORDER = 'E' 00405 END IF 00406 00407 CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, M, 00408 $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK, 00409 $ IWORK( INDIWO ), INFO ) 00410 * 00411 IF( WANTZ ) THEN 00412 CALL SSTEIN( N, D, E, M, W, IWORK( INDIBL ), IWORK( INDISP ), 00413 $ Z, LDZ, WORK, IWORK( INDIWO ), IWORK( INDIFL ), 00414 $ INFO ) 00415 END IF 00416 * 00417 * If matrix was scaled, then rescale eigenvalues appropriately. 00418 * 00419 10 CONTINUE 00420 IF( ISCALE.EQ.1 ) THEN 00421 IF( INFO.EQ.0 ) THEN 00422 IMAX = M 00423 ELSE 00424 IMAX = INFO - 1 00425 END IF 00426 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 ) 00427 END IF 00428 * 00429 * If eigenvalues are not in order, then sort them, along with 00430 * eigenvectors. 00431 * 00432 IF( WANTZ ) THEN 00433 DO 30 J = 1, M - 1 00434 I = 0 00435 TMP1 = W( J ) 00436 DO 20 JJ = J + 1, M 00437 IF( W( JJ ).LT.TMP1 ) THEN 00438 I = JJ 00439 TMP1 = W( JJ ) 00440 END IF 00441 20 CONTINUE 00442 * 00443 IF( I.NE.0 ) THEN 00444 W( I ) = W( J ) 00445 W( J ) = TMP1 00446 CALL SSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 00447 END IF 00448 30 CONTINUE 00449 END IF 00450 * 00451 * Causes problems with tests 19 & 20: 00452 * IF (wantz .and. INDEIG ) Z( 1,1) = Z(1,1) / 1.002 + .002 00453 * 00454 * 00455 WORK( 1 ) = LWMIN 00456 IWORK( 1 ) = LIWMIN 00457 RETURN 00458 * 00459 * End of SSTEVR 00460 * 00461 END
1.7.3 
