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