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