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