1 SUBROUTINE CSTEGR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, 2 $ M, W, Z, LDZ, ISUPPZ, WORK, LWORK, IWORK, 3 $ LIWORK, INFO ) 4* 5* -- LAPACK computational routine (version 3.0) -- 6* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 7* Courant Institute, Argonne National Lab, and Rice University 8* October 31, 1999 9* 10* .. Scalar Arguments .. 11 CHARACTER JOBZ, RANGE 12 INTEGER IL, INFO, IU, LDZ, LIWORK, LWORK, M, N 13 REAL ABSTOL, VL, VU 14* .. 15* .. Array Arguments .. 16 INTEGER ISUPPZ( * ), IWORK( * ) 17 REAL D( * ), E( * ), W( * ), WORK( * ) 18 COMPLEX Z( LDZ, * ) 19* .. 20* 21* Purpose 22* ======= 23* 24* CSTEGR computes selected eigenvalues and, optionally, eigenvectors 25* of a real symmetric tridiagonal matrix T. Eigenvalues and 26* eigenvectors can be selected by specifying either a range of values 27* or a range of indices for the desired eigenvalues. The eigenvalues 28* are computed by the dqds algorithm, while orthogonal eigenvectors are 29* computed from various ``good'' L D L^T representations (also known as 30* Relatively Robust Representations). Gram-Schmidt orthogonalization is 31* avoided as far as possible. More specifically, the various steps of 32* the algorithm are as follows. For the i-th unreduced block of T, 33* (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T 34* is a relatively robust representation, 35* (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high 36* relative accuracy by the dqds algorithm, 37* (c) If there is a cluster of close eigenvalues, "choose" sigma_i 38* close to the cluster, and go to step (a), 39* (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, 40* compute the corresponding eigenvector by forming a 41* rank-revealing twisted factorization. 42* The desired accuracy of the output can be specified by the input 43* parameter ABSTOL. 44* 45* For more details, see "A new O(n^2) algorithm for the symmetric 46* tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, 47* Computer Science Division Technical Report No. UCB/CSD-97-971, 48* UC Berkeley, May 1997. 49* 50* Note 1 : Currently CSTEGR is only set up to find ALL the n 51* eigenvalues and eigenvectors of T in O(n^2) time 52* Note 2 : Currently the routine CSTEIN is called when an appropriate 53* sigma_i cannot be chosen in step (c) above. CSTEIN invokes modified 54* Gram-Schmidt when eigenvalues are close. 55* Note 3 : CSTEGR works only on machines which follow ieee-754 56* floating-point standard in their handling of infinities and NaNs. 57* Normal execution of CSTEGR may create NaNs and infinities and hence 58* may abort due to a floating point exception in environments which 59* do not conform to the ieee standard. 60* 61* Arguments 62* ========= 63* 64* JOBZ (input) CHARACTER*1 65* = 'N': Compute eigenvalues only; 66* = 'V': Compute eigenvalues and eigenvectors. 67* 68* RANGE (input) CHARACTER*1 69* = 'A': all eigenvalues will be found. 70* = 'V': all eigenvalues in the half-open interval (VL,VU] 71* will be found. 72* = 'I': the IL-th through IU-th eigenvalues will be found. 73********** Only RANGE = 'A' is currently supported ********************* 74* 75* N (input) INTEGER 76* The order of the matrix. N >= 0. 77* 78* D (input/output) REAL array, dimension (N) 79* On entry, the n diagonal elements of the tridiagonal matrix 80* T. On exit, D is overwritten. 81* 82* E (input/output) REAL array, dimension (N) 83* On entry, the (n-1) subdiagonal elements of the tridiagonal 84* matrix T in elements 1 to N-1 of E; E(N) need not be set. 85* On exit, E is overwritten. 86* 87* VL (input) REAL 88* VU (input) REAL 89* If RANGE='V', the lower and upper bounds of the interval to 90* be searched for eigenvalues. VL < VU. 91* Not referenced if RANGE = 'A' or 'I'. 92* 93* IL (input) INTEGER 94* IU (input) INTEGER 95* If RANGE='I', the indices (in ascending order) of the 96* smallest and largest eigenvalues to be returned. 97* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. 98* Not referenced if RANGE = 'A' or 'V'. 99* 100* ABSTOL (input) REAL 101* The absolute error tolerance for the 102* eigenvalues/eigenvectors. IF JOBZ = 'V', the eigenvalues and 103* eigenvectors output have residual norms bounded by ABSTOL, 104* and the dot products between different eigenvectors are 105* bounded by ABSTOL. If ABSTOL is less than N*EPS*|T|, then 106* N*EPS*|T| will be used in its place, where EPS is the 107* machine precision and |T| is the 1-norm of the tridiagonal 108* matrix. The eigenvalues are computed to an accuracy of 109* EPS*|T| irrespective of ABSTOL. If high relative accuracy 110* is important, set ABSTOL to DLAMCH( 'Safe minimum' ). 111* See Barlow and Demmel "Computing Accurate Eigensystems of 112* Scaled Diagonally Dominant Matrices", LAPACK Working Note #7 113* for a discussion of which matrices define their eigenvalues 114* to high relative accuracy. 115* 116* M (output) INTEGER 117* The total number of eigenvalues found. 0 <= M <= N. 118* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. 119* 120* W (output) REAL array, dimension (N) 121* The first M elements contain the selected eigenvalues in 122* ascending order. 123* 124* Z (output) COMPLEX array, dimension (LDZ, max(1,M) ) 125* If JOBZ = 'V', then if INFO = 0, the first M columns of Z 126* contain the orthonormal eigenvectors of the matrix T 127* corresponding to the selected eigenvalues, with the i-th 128* column of Z holding the eigenvector associated with W(i). 129* If JOBZ = 'N', then Z is not referenced. 130* Note: the user must ensure that at least max(1,M) columns are 131* supplied in the array Z; if RANGE = 'V', the exact value of M 132* is not known in advance and an upper bound must be used. 133* 134* LDZ (input) INTEGER 135* The leading dimension of the array Z. LDZ >= 1, and if 136* JOBZ = 'V', LDZ >= max(1,N). 137* 138* ISUPPZ (output) INTEGER ARRAY, dimension ( 2*max(1,M) ) 139* The support of the eigenvectors in Z, i.e., the indices 140* indicating the nonzero elements in Z. The i-th eigenvector 141* is nonzero only in elements ISUPPZ( 2*i-1 ) through 142* ISUPPZ( 2*i ). 143* 144* WORK (workspace/output) REAL array, dimension (LWORK) 145* On exit, if INFO = 0, WORK(1) returns the optimal 146* (and minimal) LWORK. 147* 148* LWORK (input) INTEGER 149* The dimension of the array WORK. LWORK >= max(1,18*N) 150* 151* If LWORK = -1, then a workspace query is assumed; the routine 152* only calculates the optimal size of the WORK array, returns 153* this value as the first entry of the WORK array, and no error 154* message related to LWORK is issued by XERBLA. 155* 156* IWORK (workspace/output) INTEGER array, dimension (LIWORK) 157* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. 158* 159* LIWORK (input) INTEGER 160* The dimension of the array IWORK. LIWORK >= max(1,10*N) 161* 162* If LIWORK = -1, then a workspace query is assumed; the 163* routine only calculates the optimal size of the IWORK array, 164* returns this value as the first entry of the IWORK array, and 165* no error message related to LIWORK is issued by XERBLA. 166* 167* INFO (output) INTEGER 168* = 0: successful exit 169* < 0: if INFO = -i, the i-th argument had an illegal value 170* > 0: if INFO = 1, internal error in SLARRE, 171* if INFO = 2, internal error in CLARRV. 172* 173* Further Details 174* =============== 175* 176* Based on contributions by 177* Inderjit Dhillon, IBM Almaden, USA 178* Osni Marques, LBNL/NERSC, USA 179* Ken Stanley, Computer Science Division, University of 180* California at Berkeley, USA 181* 182* ===================================================================== 183* 184* .. Parameters .. 185 REAL ZERO, ONE 186 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 187 COMPLEX CZERO 188 PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ) ) 189* .. 190* .. Local Scalars .. 191 LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ 192 INTEGER I, IBEGIN, IEND, IINDBL, IINDWK, IINFO, IINSPL, 193 $ INDGRS, INDWOF, INDWRK, ITMP, J, JJ, LIWMIN, 194 $ LWMIN, NSPLIT 195 REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SCALE, SMLNUM, 196 $ THRESH, TMP, TNRM, TOL 197* .. 198* .. External Functions .. 199 LOGICAL LSAME 200 REAL SLAMCH, SLANST 201 EXTERNAL LSAME, SLAMCH, SLANST 202* .. 203* .. External Subroutines .. 204 EXTERNAL CLARRV, CLASET, CSWAP, SLARRE, SSCAL, XERBLA 205* .. 206* .. Intrinsic Functions .. 207 INTRINSIC MAX, MIN, REAL, SQRT 208* .. 209* .. Executable Statements .. 210* 211* Test the input parameters. 212* 213 WANTZ = LSAME( JOBZ, 'V' ) 214 ALLEIG = LSAME( RANGE, 'A' ) 215 VALEIG = LSAME( RANGE, 'V' ) 216 INDEIG = LSAME( RANGE, 'I' ) 217* 218 LQUERY = ( ( LWORK.EQ.-1 ) .OR. ( LIWORK.EQ.-1 ) ) 219 LWMIN = 18*N 220 LIWMIN = 10*N 221* 222 INFO = 0 223 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN 224 INFO = -1 225 ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN 226 INFO = -2 227* 228* The following two lines need to be removed once the 229* RANGE = 'V' and RANGE = 'I' options are provided. 230* 231 ELSE IF( VALEIG .OR. INDEIG ) THEN 232 INFO = -2 233 ELSE IF( N.LT.0 ) THEN 234 INFO = -3 235 ELSE IF( VALEIG .AND. N.GT.0 .AND. VU.LE.VL ) THEN 236 INFO = -7 237 ELSE IF( INDEIG .AND. IL.LT.1 ) THEN 238 INFO = -8 239* The following change should be made in DSTEVX also, otherwise 240* IL can be specified as N+1 and IU as N. 241* ELSE IF( INDEIG .AND. ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN 242 ELSE IF( INDEIG .AND. ( IU.LT.IL .OR. IU.GT.N ) ) THEN 243 INFO = -9 244 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN 245 INFO = -14 246 ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN 247 INFO = -17 248 ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN 249 INFO = -19 250 END IF 251 IF( INFO.EQ.0 ) THEN 252 WORK( 1 ) = LWMIN 253 IWORK( 1 ) = LIWMIN 254 END IF 255* 256 IF( INFO.NE.0 ) THEN 257 CALL XERBLA( 'CSTEGR', -INFO ) 258 RETURN 259 ELSE IF( LQUERY ) THEN 260 RETURN 261 END IF 262* 263* Quick return if possible 264* 265 M = 0 266 IF( N.EQ.0 ) 267 $ RETURN 268* 269 IF( N.EQ.1 ) THEN 270 IF( ALLEIG .OR. INDEIG ) THEN 271 M = 1 272 W( 1 ) = D( 1 ) 273 ELSE 274 IF( VL.LT.D( 1 ) .AND. VU.GE.D( 1 ) ) THEN 275 M = 1 276 W( 1 ) = D( 1 ) 277 END IF 278 END IF 279 IF( WANTZ ) 280 $ Z( 1, 1 ) = ONE 281 RETURN 282 END IF 283* 284* Get machine constants. 285* 286 SAFMIN = SLAMCH( 'Safe minimum' ) 287 EPS = SLAMCH( 'Precision' ) 288 SMLNUM = SAFMIN / EPS 289 BIGNUM = ONE / SMLNUM 290 RMIN = SQRT( SMLNUM ) 291 RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) ) 292* 293* Scale matrix to allowable range, if necessary. 294* 295 SCALE = ONE 296 TNRM = SLANST( 'M', N, D, E ) 297 IF( TNRM.GT.ZERO .AND. TNRM.LT.RMIN ) THEN 298 SCALE = RMIN / TNRM 299 ELSE IF( TNRM.GT.RMAX ) THEN 300 SCALE = RMAX / TNRM 301 END IF 302 IF( SCALE.NE.ONE ) THEN 303 CALL SSCAL( N, SCALE, D, 1 ) 304 CALL SSCAL( N-1, SCALE, E, 1 ) 305 TNRM = TNRM*SCALE 306 END IF 307 INDGRS = 1 308 INDWOF = 2*N + 1 309 INDWRK = 3*N + 1 310* 311 IINSPL = 1 312 IINDBL = N + 1 313 IINDWK = 2*N + 1 314* 315 CALL CLASET( 'Full', N, N, CZERO, CZERO, Z, LDZ ) 316* 317* Compute the desired eigenvalues of the tridiagonal after splitting 318* into smaller subblocks if the corresponding of-diagonal elements 319* are small 320* 321 THRESH = EPS*TNRM 322 CALL SLARRE( N, D, E, THRESH, NSPLIT, IWORK( IINSPL ), M, W, 323 $ WORK( INDWOF ), WORK( INDGRS ), WORK( INDWRK ), 324 $ IINFO ) 325 IF( IINFO.NE.0 ) THEN 326 INFO = 1 327 RETURN 328 END IF 329* 330 IF( WANTZ ) THEN 331* 332* Compute the desired eigenvectors corresponding to the computed 333* eigenvalues 334* 335 TOL = MAX( ABSTOL, REAL( N )*THRESH ) 336 IBEGIN = 1 337 DO 20 I = 1, NSPLIT 338 IEND = IWORK( IINSPL+I-1 ) 339 DO 10 J = IBEGIN, IEND 340 IWORK( IINDBL+J-1 ) = I 341 10 CONTINUE 342 IBEGIN = IEND + 1 343 20 CONTINUE 344* 345 CALL CLARRV( N, D, E, IWORK( IINSPL ), M, W, IWORK( IINDBL ), 346 $ WORK( INDGRS ), TOL, Z, LDZ, ISUPPZ, 347 $ WORK( INDWRK ), IWORK( IINDWK ), IINFO ) 348 IF( IINFO.NE.0 ) THEN 349 INFO = 2 350 RETURN 351 END IF 352* 353 END IF 354* 355 IBEGIN = 1 356 DO 40 I = 1, NSPLIT 357 IEND = IWORK( IINSPL+I-1 ) 358 DO 30 J = IBEGIN, IEND 359 W( J ) = W( J ) + WORK( INDWOF+I-1 ) 360 30 CONTINUE 361 IBEGIN = IEND + 1 362 40 CONTINUE 363* 364* If matrix was scaled, then rescale eigenvalues appropriately. 365* 366 IF( SCALE.NE.ONE ) THEN 367 CALL SSCAL( M, ONE / SCALE, W, 1 ) 368 END IF 369* 370* If eigenvalues are not in order, then sort them, along with 371* eigenvectors. 372* 373 IF( NSPLIT.GT.1 ) THEN 374 DO 60 J = 1, M - 1 375 I = 0 376 TMP = W( J ) 377 DO 50 JJ = J + 1, M 378 IF( W( JJ ).LT.TMP ) THEN 379 I = JJ 380 TMP = W( JJ ) 381 END IF 382 50 CONTINUE 383 IF( I.NE.0 ) THEN 384 W( I ) = W( J ) 385 W( J ) = TMP 386 IF( WANTZ ) THEN 387 CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 ) 388 ITMP = ISUPPZ( 2*I-1 ) 389 ISUPPZ( 2*I-1 ) = ISUPPZ( 2*J-1 ) 390 ISUPPZ( 2*J-1 ) = ITMP 391 ITMP = ISUPPZ( 2*I ) 392 ISUPPZ( 2*I ) = ISUPPZ( 2*J ) 393 ISUPPZ( 2*J ) = ITMP 394 END IF 395 END IF 396 60 CONTINUE 397 END IF 398* 399 WORK( 1 ) = LWMIN 400 IWORK( 1 ) = LIWMIN 401 RETURN 402* 403* End of CSTEGR 404* 405 END 406