1      subroutine imtql1(n,d,e,ierr)
2c
3      integer i,j,l,m,n,ii,mml,ierr
4      double precision d(n),e(n)
5      double precision b,c,f,g,p,r,s,tst1,tst2,pythag
6c
7c     this subroutine is a translation of the algol procedure imtql1,
8c     num. math. 12, 377-383(1968) by martin and wilkinson,
9c     as modified in num. math. 15, 450(1970) by dubrulle.
10c     handbook for auto. comp., vol.ii-linear algebra, 241-248(1971).
11c
12c     this subroutine finds the eigenvalues of a symmetric
13c     tridiagonal matrix by the implicit ql method.
14c
15c     on input
16c
17c        n is the order of the matrix.
18c
19c        d contains the diagonal elements of the input matrix.
20c
21c        e contains the subdiagonal elements of the input matrix
22c          in its last n-1 positions.  e(1) is arbitrary.
23c
24c      on output
25c
26c        d contains the eigenvalues in ascending order.  if an
27c          error exit is made, the eigenvalues are correct and
28c          ordered for indices 1,2,...ierr-1, but may not be
29c          the smallest eigenvalues.
30c
31c        e has been destroyed.
32c
33c        ierr is set to
34c          zero       for normal return,
35c          j          if the j-th eigenvalue has not been
36c                     determined after 30 iterations.
37c
38c     calls pythag for  dsqrt(a*a + b*b) .
39c
40c     questions and comments should be directed to burton s. garbow,
41c     mathematics and computer science div, argonne national laboratory
42c
43c     this version dated august 1983.
44c
45c     ------------------------------------------------------------------
46c
47      ierr = 0
48      if (n .eq. 1) go to 1001
49c
50      do 100 i = 2, n
51  100 e(i-1) = e(i)
52c
53      e(n) = 0.0d0
54c
55      do 290 l = 1, n
56         j = 0
57c     .......... look for small sub-diagonal element ..........
58  105    do 110 m = l, n
59            if (m .eq. n) go to 120
60            tst1 = dabs(d(m)) + dabs(d(m+1))
61            tst2 = tst1 + dabs(e(m))
62            if (tst2 .eq. tst1) go to 120
63  110    continue
64c
65  120    p = d(l)
66         if (m .eq. l) go to 215
67         if (j .eq. 30) go to 1000
68         j = j + 1
69c     .......... form shift ..........
70         g = (d(l+1) - p) / (2.0d0 * e(l))
71         r = pythag(g,1.0d0)
72         g = d(m) - p + e(l) / (g + dsign(r,g))
73         s = 1.0d0
74         c = 1.0d0
75         p = 0.0d0
76         mml = m - l
77c     .......... for i=m-1 step -1 until l do -- ..........
78         do 200 ii = 1, mml
79            i = m - ii
80            f = s * e(i)
81            b = c * e(i)
82            r = pythag(f,g)
83            e(i+1) = r
84            if (r .eq. 0.0d0) go to 210
85            s = f / r
86            c = g / r
87            g = d(i+1) - p
88            r = (d(i) - g) * s + 2.0d0 * c * b
89            p = s * r
90            d(i+1) = g + p
91            g = c * r - b
92  200    continue
93c
94         d(l) = d(l) - p
95         e(l) = g
96         e(m) = 0.0d0
97         go to 105
98c     .......... recover from underflow ..........
99  210    d(i+1) = d(i+1) - p
100         e(m) = 0.0d0
101         go to 105
102c     .......... order eigenvalues ..........
103  215    if (l .eq. 1) go to 250
104c     .......... for i=l step -1 until 2 do -- ..........
105         do 230 ii = 2, l
106            i = l + 2 - ii
107            if (p .ge. d(i-1)) go to 270
108            d(i) = d(i-1)
109  230    continue
110c
111  250    i = 1
112  270    d(i) = p
113  290 continue
114c
115      go to 1001
116c     .......... set error -- no convergence to an
117c                eigenvalue after 30 iterations ..........
118 1000 ierr = l
119 1001 return
120      end
121