1*DECK TSTURM 2 SUBROUTINE TSTURM (NM, N, EPS1, D, E, E2, LB, UB, MM, M, W, Z, 3 + IERR, RV1, RV2, RV3, RV4, RV5, RV6) 4C***BEGIN PROLOGUE TSTURM 5C***PURPOSE Find those eigenvalues of a symmetric tridiagonal matrix 6C in a given interval and their associated eigenvectors by 7C Sturm sequencing. 8C***LIBRARY SLATEC (EISPACK) 9C***CATEGORY D4A5, D4C2A 10C***TYPE SINGLE PRECISION (TSTURM-S) 11C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK 12C***AUTHOR Smith, B. T., et al. 13C***DESCRIPTION 14C 15C This subroutine finds those eigenvalues of a TRIDIAGONAL 16C SYMMETRIC matrix which lie in a specified interval and their 17C associated eigenvectors, using bisection and inverse iteration. 18C 19C On Input 20C 21C NM must be set to the row dimension of the two-dimensional 22C array parameter, Z, as declared in the calling program 23C dimension statement. NM is an INTEGER variable. 24C 25C N is the order of the matrix. N is an INTEGER variable. 26C N must be less than or equal to NM. 27C 28C EPS1 is an absolute error tolerance for the computed eigen- 29C values. It should be chosen so that the accuracy of these 30C eigenvalues is commensurate with relative perturbations of 31C the order of the relative machine precision in the matrix 32C elements. If the input EPS1 is non-positive, it is reset 33C for each submatrix to a default value, namely, minus the 34C product of the relative machine precision and the 1-norm of 35C the submatrix. EPS1 is a REAL variable. 36C 37C D contains the diagonal elements of the symmetric tridiagonal 38C matrix. D is a one-dimensional REAL array, dimensioned D(N). 39C 40C E contains the subdiagonal elements of the symmetric 41C tridiagonal matrix in its last N-1 positions. E(1) is 42C arbitrary. E is a one-dimensional REAL array, dimensioned 43C E(N). 44C 45C E2 contains the squares of the corresponding elements of E. 46C E2(1) is arbitrary. E2 is a one-dimensional REAL array, 47C dimensioned E2(N). 48C 49C LB and UB define the interval to be searched for eigenvalues. 50C If LB is not less than UB, no eigenvalues will be found. 51C LB and UB are REAL variables. 52C 53C MM should be set to an upper bound for the number of 54C eigenvalues in the interval. MM is an INTEGER variable. 55C WARNING - If more than MM eigenvalues are determined to lie 56C in the interval, an error return is made with no values or 57C vectors found. 58C 59C On Output 60C 61C EPS1 is unaltered unless it has been reset to its 62C (last) default value. 63C 64C D and E are unaltered. 65C 66C Elements of E2, corresponding to elements of E regarded as 67C negligible, have been replaced by zero causing the matrix to 68C split into a direct sum of submatrices. E2(1) is also set 69C to zero. 70C 71C M is the number of eigenvalues determined to lie in (LB,UB). 72C M is an INTEGER variable. 73C 74C W contains the M eigenvalues in ascending order if the matrix 75C does not split. If the matrix splits, the eigenvalues are 76C in ascending order for each submatrix. If a vector error 77C exit is made, W contains those values already found. W is a 78C one-dimensional REAL array, dimensioned W(MM). 79C 80C Z contains the associated set of orthonormal eigenvectors. 81C If an error exit is made, Z contains those vectors already 82C found. Z is a one-dimensional REAL array, dimensioned 83C Z(NM,MM). 84C 85C IERR is an INTEGER flag set to 86C Zero for normal return, 87C 3*N+1 if M exceeds MM no eigenvalues or eigenvectors 88C are computed, 89C 4*N+J if the eigenvector corresponding to the J-th 90C eigenvalue fails to converge in 5 iterations, then 91C the eigenvalues and eigenvectors in W and Z should 92C be correct for indices 1, 2, ..., J-1. 93C 94C RV1, RV2, RV3, RV4, RV5, and RV6 are temporary storage arrays, 95C dimensioned RV1(N), RV2(N), RV3(N), RV4(N), RV5(N), and 96C RV6(N). 97C 98C The ALGOL procedure STURMCNT contained in TRISTURM 99C appears in TSTURM in-line. 100C 101C Questions and comments should be directed to B. S. Garbow, 102C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY 103C ------------------------------------------------------------------ 104C 105C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, 106C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- 107C system Routines - EISPACK Guide, Springer-Verlag, 108C 1976. 109C***ROUTINES CALLED R1MACH 110C***REVISION HISTORY (YYMMDD) 111C 760101 DATE WRITTEN 112C 890531 Changed all specific intrinsics to generic. (WRB) 113C 890531 REVISION DATE from Version 3.2 114C 891214 Prologue converted to Version 4.0 format. (BAB) 115C 920501 Reformatted the REFERENCES section. (WRB) 116C***END PROLOGUE TSTURM 117C 118 INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS 119 INTEGER IERR,GROUP,ISTURM 120 REAL D(*),E(*),E2(*),W(*),Z(NM,*) 121 REAL RV1(*),RV2(*),RV3(*),RV4(*),RV5(*),RV6(*) 122 REAL U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4 123 REAL NORM,MACHEP,S1,S2 124 LOGICAL FIRST 125C 126 SAVE FIRST, MACHEP 127 DATA FIRST /.TRUE./ 128C***FIRST EXECUTABLE STATEMENT TSTURM 129 IF (FIRST) THEN 130 MACHEP = R1MACH(4) 131 ENDIF 132 FIRST = .FALSE. 133C 134 IERR = 0 135 T1 = LB 136 T2 = UB 137C .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES .......... 138 DO 40 I = 1, N 139 IF (I .EQ. 1) GO TO 20 140 S1 = ABS(D(I)) + ABS(D(I-1)) 141 S2 = S1 + ABS(E(I)) 142 IF (S2 .GT. S1) GO TO 40 143 20 E2(I) = 0.0E0 144 40 CONTINUE 145C .......... DETERMINE THE NUMBER OF EIGENVALUES 146C IN THE INTERVAL .......... 147 P = 1 148 Q = N 149 X1 = UB 150 ISTURM = 1 151 GO TO 320 152 60 M = S 153 X1 = LB 154 ISTURM = 2 155 GO TO 320 156 80 M = M - S 157 IF (M .GT. MM) GO TO 980 158 Q = 0 159 R = 0 160C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING 161C INTERVAL BY THE GERSCHGORIN BOUNDS .......... 162 100 IF (R .EQ. M) GO TO 1001 163 P = Q + 1 164 XU = D(P) 165 X0 = D(P) 166 U = 0.0E0 167C 168 DO 120 Q = P, N 169 X1 = U 170 U = 0.0E0 171 V = 0.0E0 172 IF (Q .EQ. N) GO TO 110 173 U = ABS(E(Q+1)) 174 V = E2(Q+1) 175 110 XU = MIN(D(Q)-(X1+U),XU) 176 X0 = MAX(D(Q)+(X1+U),X0) 177 IF (V .EQ. 0.0E0) GO TO 140 178 120 CONTINUE 179C 180 140 X1 = MAX(ABS(XU),ABS(X0)) * MACHEP 181 IF (EPS1 .LE. 0.0E0) EPS1 = -X1 182 IF (P .NE. Q) GO TO 180 183C .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... 184 IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940 185 R = R + 1 186C 187 DO 160 I = 1, N 188 160 Z(I,R) = 0.0E0 189C 190 W(R) = D(P) 191 Z(P,R) = 1.0E0 192 GO TO 940 193 180 X1 = X1 * (Q-P+1) 194 LB = MAX(T1,XU-X1) 195 UB = MIN(T2,X0+X1) 196 X1 = LB 197 ISTURM = 3 198 GO TO 320 199 200 M1 = S + 1 200 X1 = UB 201 ISTURM = 4 202 GO TO 320 203 220 M2 = S 204 IF (M1 .GT. M2) GO TO 940 205C .......... FIND ROOTS BY BISECTION .......... 206 X0 = UB 207 ISTURM = 5 208C 209 DO 240 I = M1, M2 210 RV5(I) = UB 211 RV4(I) = LB 212 240 CONTINUE 213C .......... LOOP FOR K-TH EIGENVALUE 214C FOR K=M2 STEP -1 UNTIL M1 DO -- 215C (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) .......... 216 K = M2 217 250 XU = LB 218C .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... 219 DO 260 II = M1, K 220 I = M1 + K - II 221 IF (XU .GE. RV4(I)) GO TO 260 222 XU = RV4(I) 223 GO TO 280 224 260 CONTINUE 225C 226 280 IF (X0 .GT. RV5(K)) X0 = RV5(K) 227C .......... NEXT BISECTION STEP .......... 228 300 X1 = (XU + X0) * 0.5E0 229 S1 = 2.0E0*(ABS(XU) + ABS(X0) + ABS(EPS1)) 230 S2 = S1 + ABS(X0 - XU) 231 IF (S2 .EQ. S1) GO TO 420 232C .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... 233 320 S = P - 1 234 U = 1.0E0 235C 236 DO 340 I = P, Q 237 IF (U .NE. 0.0E0) GO TO 325 238 V = ABS(E(I)) / MACHEP 239 IF (E2(I) .EQ. 0.0E0) V = 0.0E0 240 GO TO 330 241 325 V = E2(I) / U 242 330 U = D(I) - X1 - V 243 IF (U .LT. 0.0E0) S = S + 1 244 340 CONTINUE 245C 246 GO TO (60,80,200,220,360), ISTURM 247C .......... REFINE INTERVALS .......... 248 360 IF (S .GE. K) GO TO 400 249 XU = X1 250 IF (S .GE. M1) GO TO 380 251 RV4(M1) = X1 252 GO TO 300 253 380 RV4(S+1) = X1 254 IF (RV5(S) .GT. X1) RV5(S) = X1 255 GO TO 300 256 400 X0 = X1 257 GO TO 300 258C .......... K-TH EIGENVALUE FOUND .......... 259 420 RV5(K) = X1 260 K = K - 1 261 IF (K .GE. M1) GO TO 250 262C .......... FIND VECTORS BY INVERSE ITERATION .......... 263 NORM = ABS(D(P)) 264 IP = P + 1 265C 266 DO 500 I = IP, Q 267 500 NORM = MAX(NORM, ABS(D(I)) + ABS(E(I))) 268C .......... EPS2 IS THE CRITERION FOR GROUPING, 269C EPS3 REPLACES ZERO PIVOTS AND EQUAL 270C ROOTS ARE MODIFIED BY EPS3, 271C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .......... 272 EPS2 = 1.0E-3 * NORM 273 UK = SQRT(REAL(Q-P+5)) 274 EPS3 = UK * MACHEP * NORM 275 EPS4 = UK * EPS3 276 UK = EPS4 / SQRT(UK) 277 GROUP = 0 278 S = P 279C 280 DO 920 K = M1, M2 281 R = R + 1 282 ITS = 1 283 W(R) = RV5(K) 284 X1 = RV5(K) 285C .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... 286 IF (K .EQ. M1) GO TO 520 287 IF (X1 - X0 .GE. EPS2) GROUP = -1 288 GROUP = GROUP + 1 289 IF (X1 .LE. X0) X1 = X0 + EPS3 290C .......... ELIMINATION WITH INTERCHANGES AND 291C INITIALIZATION OF VECTOR .......... 292 520 V = 0.0E0 293C 294 DO 580 I = P, Q 295 RV6(I) = UK 296 IF (I .EQ. P) GO TO 560 297 IF (ABS(E(I)) .LT. ABS(U)) GO TO 540 298 XU = U / E(I) 299 RV4(I) = XU 300 RV1(I-1) = E(I) 301 RV2(I-1) = D(I) - X1 302 RV3(I-1) = 0.0E0 303 IF (I .NE. Q) RV3(I-1) = E(I+1) 304 U = V - XU * RV2(I-1) 305 V = -XU * RV3(I-1) 306 GO TO 580 307 540 XU = E(I) / U 308 RV4(I) = XU 309 RV1(I-1) = U 310 RV2(I-1) = V 311 RV3(I-1) = 0.0E0 312 560 U = D(I) - X1 - XU * V 313 IF (I .NE. Q) V = E(I+1) 314 580 CONTINUE 315C 316 IF (U .EQ. 0.0E0) U = EPS3 317 RV1(Q) = U 318 RV2(Q) = 0.0E0 319 RV3(Q) = 0.0E0 320C .......... BACK SUBSTITUTION 321C FOR I=Q STEP -1 UNTIL P DO -- .......... 322 600 DO 620 II = P, Q 323 I = P + Q - II 324 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I) 325 V = U 326 U = RV6(I) 327 620 CONTINUE 328C .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS 329C MEMBERS OF GROUP .......... 330 IF (GROUP .EQ. 0) GO TO 700 331C 332 DO 680 JJ = 1, GROUP 333 J = R - GROUP - 1 + JJ 334 XU = 0.0E0 335C 336 DO 640 I = P, Q 337 640 XU = XU + RV6(I) * Z(I,J) 338C 339 DO 660 I = P, Q 340 660 RV6(I) = RV6(I) - XU * Z(I,J) 341C 342 680 CONTINUE 343C 344 700 NORM = 0.0E0 345C 346 DO 720 I = P, Q 347 720 NORM = NORM + ABS(RV6(I)) 348C 349 IF (NORM .GE. 1.0E0) GO TO 840 350C .......... FORWARD SUBSTITUTION .......... 351 IF (ITS .EQ. 5) GO TO 960 352 IF (NORM .NE. 0.0E0) GO TO 740 353 RV6(S) = EPS4 354 S = S + 1 355 IF (S .GT. Q) S = P 356 GO TO 780 357 740 XU = EPS4 / NORM 358C 359 DO 760 I = P, Q 360 760 RV6(I) = RV6(I) * XU 361C .......... ELIMINATION OPERATIONS ON NEXT VECTOR 362C ITERATE .......... 363 780 DO 820 I = IP, Q 364 U = RV6(I) 365C .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE 366C WAS PERFORMED EARLIER IN THE 367C TRIANGULARIZATION PROCESS .......... 368 IF (RV1(I-1) .NE. E(I)) GO TO 800 369 U = RV6(I-1) 370 RV6(I-1) = RV6(I) 371 800 RV6(I) = U - RV4(I) * RV6(I-1) 372 820 CONTINUE 373C 374 ITS = ITS + 1 375 GO TO 600 376C .......... NORMALIZE SO THAT SUM OF SQUARES IS 377C 1 AND EXPAND TO FULL ORDER .......... 378 840 U = 0.0E0 379C 380 DO 860 I = P, Q 381 860 U = U + RV6(I)**2 382C 383 XU = 1.0E0 / SQRT(U) 384C 385 DO 880 I = 1, N 386 880 Z(I,R) = 0.0E0 387C 388 DO 900 I = P, Q 389 900 Z(I,R) = RV6(I) * XU 390C 391 X0 = X1 392 920 CONTINUE 393C 394 940 IF (Q .LT. N) GO TO 100 395 GO TO 1001 396C .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... 397 960 IERR = 4 * N + R 398 GO TO 1001 399C .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF 400C EIGENVALUES IN INTERVAL .......... 401 980 IERR = 3 * N + 1 402 1001 LB = T1 403 UB = T2 404 RETURN 405 END 406