1      subroutine tridib(n,eps1,d,e,e2,lb,ub,m11,m,w,ind,ierr,rv4,rv5)
2c
3      integer i,j,k,l,m,n,p,q,r,s,ii,m1,m2,m11,m22,tag,ierr,isturm
4      double precision d(n),e(n),e2(n),w(m),rv4(n),rv5(n)
5      double precision u,v,lb,t1,t2,ub,xu,x0,x1,eps1,tst1,tst2,epslon
6      integer ind(m)
7c
8c     this subroutine is a translation of the algol procedure bisect,
9c     num. math. 9, 386-393(1967) by barth, martin, and wilkinson.
10c     handbook for auto. comp., vol.ii-linear algebra, 249-256(1971).
11c
12c     this subroutine finds those eigenvalues of a tridiagonal
13c     symmetric matrix between specified boundary indices,
14c     using bisection.
15c
16c     on input
17c
18c        n is the order of the matrix.
19c
20c        eps1 is an absolute error tolerance for the computed
21c          eigenvalues.  if the input eps1 is non-positive,
22c          it is reset for each submatrix to a default value,
23c          namely, minus the product of the relative machine
24c          precision and the 1-norm of the submatrix.
25c
26c        d contains the diagonal elements of the input matrix.
27c
28c        e contains the subdiagonal elements of the input matrix
29c          in its last n-1 positions.  e(1) is arbitrary.
30c
31c        e2 contains the squares of the corresponding elements of e.
32c          e2(1) is arbitrary.
33c
34c        m11 specifies the lower boundary index for the desired
35c          eigenvalues.
36c
37c        m specifies the number of eigenvalues desired.  the upper
38c          boundary index m22 is then obtained as m22=m11+m-1.
39c
40c     on output
41c
42c        eps1 is unaltered unless it has been reset to its
43c          (last) default value.
44c
45c        d and e are unaltered.
46c
47c        elements of e2, corresponding to elements of e regarded
48c          as negligible, have been replaced by zero causing the
49c          matrix to split into a direct sum of submatrices.
50c          e2(1) is also set to zero.
51c
52c        lb and ub define an interval containing exactly the desired
53c          eigenvalues.
54c
55c        w contains, in its first m positions, the eigenvalues
56c          between indices m11 and m22 in ascending order.
57c
58c        ind contains in its first m positions the submatrix indices
59c          associated with the corresponding eigenvalues in w --
60c          1 for eigenvalues belonging to the first submatrix from
61c          the top, 2 for those belonging to the second submatrix, etc..
62c
63c        ierr is set to
64c          zero       for normal return,
65c          3*n+1      if multiple eigenvalues at index m11 make
66c                     unique selection impossible,
67c          3*n+2      if multiple eigenvalues at index m22 make
68c                     unique selection impossible.
69c
70c        rv4 and rv5 are temporary storage arrays.
71c
72c     note that subroutine tql1, imtql1, or tqlrat is generally faster
73c     than tridib, if more than n/4 eigenvalues are to be found.
74c
75c     questions and comments should be directed to burton s. garbow,
76c     mathematics and computer science div, argonne national laboratory
77c
78c     this version dated august 1983.
79c
80c     ------------------------------------------------------------------
81c
82      ierr = 0
83      tag = 0
84      xu = d(1)
85      x0 = d(1)
86      u = 0.0d0
87c     .......... look for small sub-diagonal entries and determine an
88c                interval containing all the eigenvalues ..........
89      do 40 i = 1, n
90         x1 = u
91         u = 0.0d0
92         if (i .ne. n) u = dabs(e(i+1))
93         xu = dmin1(d(i)-(x1+u),xu)
94         x0 = dmax1(d(i)+(x1+u),x0)
95         if (i .eq. 1) go to 20
96         tst1 = dabs(d(i)) + dabs(d(i-1))
97         tst2 = tst1 + dabs(e(i))
98         if (tst2 .gt. tst1) go to 40
99   20    e2(i) = 0.0d0
100   40 continue
101c
102      x1 = n
103      x1 = x1 * epslon(dmax1(dabs(xu),dabs(x0)))
104      xu = xu - x1
105      t1 = xu
106      x0 = x0 + x1
107      t2 = x0
108c     .......... determine an interval containing exactly
109c                the desired eigenvalues ..........
110      p = 1
111      q = n
112      m1 = m11 - 1
113      if (m1 .eq. 0) go to 75
114      isturm = 1
115   50 v = x1
116      x1 = xu + (x0 - xu) * 0.5d0
117      if (x1 .eq. v) go to 980
118      go to 320
119   60 if (s - m1) 65, 73, 70
120   65 xu = x1
121      go to 50
122   70 x0 = x1
123      go to 50
124   73 xu = x1
125      t1 = x1
126   75 m22 = m1 + m
127      if (m22 .eq. n) go to 90
128      x0 = t2
129      isturm = 2
130      go to 50
131   80 if (s - m22) 65, 85, 70
132   85 t2 = x1
133   90 q = 0
134      r = 0
135c     .......... establish and process next submatrix, refining
136c                interval by the gerschgorin bounds ..........
137  100 if (r .eq. m) go to 1001
138      tag = tag + 1
139      p = q + 1
140      xu = d(p)
141      x0 = d(p)
142      u = 0.0d0
143c
144      do 120 q = p, n
145         x1 = u
146         u = 0.0d0
147         v = 0.0d0
148         if (q .eq. n) go to 110
149         u = dabs(e(q+1))
150         v = e2(q+1)
151  110    xu = dmin1(d(q)-(x1+u),xu)
152         x0 = dmax1(d(q)+(x1+u),x0)
153         if (v .eq. 0.0d0) go to 140
154  120 continue
155c
156  140 x1 = epslon(dmax1(dabs(xu),dabs(x0)))
157      if (eps1 .le. 0.0d0) eps1 = -x1
158      if (p .ne. q) go to 180
159c     .......... check for isolated root within interval ..........
160      if (t1 .gt. d(p) .or. d(p) .ge. t2) go to 940
161      m1 = p
162      m2 = p
163      rv5(p) = d(p)
164      go to 900
165  180 x1 = x1 * (q - p + 1)
166      lb = dmax1(t1,xu-x1)
167      ub = dmin1(t2,x0+x1)
168      x1 = lb
169      isturm = 3
170      go to 320
171  200 m1 = s + 1
172      x1 = ub
173      isturm = 4
174      go to 320
175  220 m2 = s
176      if (m1 .gt. m2) go to 940
177c     .......... find roots by bisection ..........
178      x0 = ub
179      isturm = 5
180c
181      do 240 i = m1, m2
182         rv5(i) = ub
183         rv4(i) = lb
184  240 continue
185c     .......... loop for k-th eigenvalue
186c                for k=m2 step -1 until m1 do --
187c                (-do- not used to legalize -computed go to-) ..........
188      k = m2
189  250    xu = lb
190c     .......... for i=k step -1 until m1 do -- ..........
191         do 260 ii = m1, k
192            i = m1 + k - ii
193            if (xu .ge. rv4(i)) go to 260
194            xu = rv4(i)
195            go to 280
196  260    continue
197c
198  280    if (x0 .gt. rv5(k)) x0 = rv5(k)
199c     .......... next bisection step ..........
200  300    x1 = (xu + x0) * 0.5d0
201         if ((x0 - xu) .le. dabs(eps1)) go to 420
202         tst1 = 2.0d0 * (dabs(xu) + dabs(x0))
203         tst2 = tst1 + (x0 - xu)
204         if (tst2 .eq. tst1) go to 420
205c     .......... in-line procedure for sturm sequence ..........
206  320    s = p - 1
207         u = 1.0d0
208c
209         do 340 i = p, q
210            if (u .ne. 0.0d0) go to 325
211            v = dabs(e(i)) / epslon(1.0d0)
212            if (e2(i) .eq. 0.0d0) v = 0.0d0
213            go to 330
214  325       v = e2(i) / u
215  330       u = d(i) - x1 - v
216            if (u .lt. 0.0d0) s = s + 1
217  340    continue
218c
219         go to (60,80,200,220,360), isturm
220c     .......... refine intervals ..........
221  360    if (s .ge. k) go to 400
222         xu = x1
223         if (s .ge. m1) go to 380
224         rv4(m1) = x1
225         go to 300
226  380    rv4(s+1) = x1
227         if (rv5(s) .gt. x1) rv5(s) = x1
228         go to 300
229  400    x0 = x1
230         go to 300
231c     .......... k-th eigenvalue found ..........
232  420    rv5(k) = x1
233      k = k - 1
234      if (k .ge. m1) go to 250
235c     .......... order eigenvalues tagged with their
236c                submatrix associations ..........
237  900 s = r
238      r = r + m2 - m1 + 1
239      j = 1
240      k = m1
241c
242      do 920 l = 1, r
243         if (j .gt. s) go to 910
244         if (k .gt. m2) go to 940
245         if (rv5(k) .ge. w(l)) go to 915
246c
247         do 905 ii = j, s
248            i = l + s - ii
249            w(i+1) = w(i)
250            ind(i+1) = ind(i)
251  905    continue
252c
253  910    w(l) = rv5(k)
254         ind(l) = tag
255         k = k + 1
256         go to 920
257  915    j = j + 1
258  920 continue
259c
260  940 if (q .lt. n) go to 100
261      go to 1001
262c     .......... set error -- interval cannot be found containing
263c                exactly the desired eigenvalues ..........
264  980 ierr = 3 * n + isturm
265 1001 lb = t1
266      ub = t2
267      return
268      end
269