subroutine imtql1(n,d,e,ierr) c integer i,j,l,m,n,ii,mml,ierr double precision d(n),e(n) double precision b,c,f,g,p,r,s,tst1,tst2,pythag c c this subroutine is a translation of the algol procedure imtql1, c num. math. 12, 377-383(1968) by martin and wilkinson, c as modified in num. math. 15, 450(1970) by dubrulle. c handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). c c this subroutine finds the eigenvalues of a symmetric c tridiagonal matrix by the implicit ql method. c c on input c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues. c c e has been destroyed. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e(i-1) = e(i) c e(n) = 0.0d0 c do 290 l = 1, n j = 0 c .......... look for small sub-diagonal element .......... 105 do 110 m = l, n if (m .eq. n) go to 120 tst1 = dabs(d(m)) + dabs(d(m+1)) tst2 = tst1 + dabs(e(m)) if (tst2 .eq. tst1) go to 120 110 continue c 120 p = d(l) if (m .eq. l) go to 215 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... g = (d(l+1) - p) / (2.0d0 * e(l)) r = pythag(g,1.0d0) g = d(m) - p + e(l) / (g + dsign(r,g)) s = 1.0d0 c = 1.0d0 p = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml i = m - ii f = s * e(i) b = c * e(i) r = pythag(f,g) e(i+1) = r if (r .eq. 0.0d0) go to 210 s = f / r c = g / r g = d(i+1) - p r = (d(i) - g) * s + 2.0d0 * c * b p = s * r d(i+1) = g + p g = c * r - b 200 continue c d(l) = d(l) - p e(l) = g e(m) = 0.0d0 go to 105 c .......... recover from underflow .......... 210 d(i+1) = d(i+1) - p e(m) = 0.0d0 go to 105 c .......... order eigenvalues .......... 215 if (l .eq. 1) go to 250 c .......... for i=l step -1 until 2 do -- .......... do 230 ii = 2, l i = l + 2 - ii if (p .ge. d(i-1)) go to 270 d(i) = d(i-1) 230 continue c 250 i = 1 270 d(i) = p 290 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end