tql1.f
135 lines
| 3.7 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | subroutine tql1(n,d,e,ierr) | |||
| c | ||||
| integer i,j,l,m,n,ii,l1,l2,mml,ierr | ||||
| real d(n),e(n) | ||||
| real c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag | ||||
| c | ||||
| c this subroutine is a translation of the algol procedure tql1, | ||||
| c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and | ||||
| c wilkinson. | ||||
| c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). | ||||
| c | ||||
| c this subroutine finds the eigenvalues of a symmetric | ||||
| c tridiagonal matrix by the 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 sqrt(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 | ||||
| f = 0.0e0 | ||||
| tst1 = 0.0e0 | ||||
| e(n) = 0.0e0 | ||||
| c | ||||
| do 290 l = 1, n | ||||
| j = 0 | ||||
| h = abs(d(l)) + abs(e(l)) | ||||
| if (tst1 .lt. h) tst1 = h | ||||
| c .......... look for small sub-diagonal element .......... | ||||
| do 110 m = l, n | ||||
| tst2 = tst1 + abs(e(m)) | ||||
| if (tst2 .eq. tst1) go to 120 | ||||
| c .......... e(n) is always zero, so there is no exit | ||||
| c through the bottom of the loop .......... | ||||
| 110 continue | ||||
| c | ||||
| 120 if (m .eq. l) go to 210 | ||||
| 130 if (j .eq. 30) go to 1000 | ||||
| j = j + 1 | ||||
| c .......... form shift .......... | ||||
| l1 = l + 1 | ||||
| l2 = l1 + 1 | ||||
| g = d(l) | ||||
| p = (d(l1) - g) / (2.0e0 * e(l)) | ||||
| r = pythag(p,1.0e0) | ||||
| d(l) = e(l) / (p + sign(r,p)) | ||||
| d(l1) = e(l) * (p + sign(r,p)) | ||||
| dl1 = d(l1) | ||||
| h = g - d(l) | ||||
| if (l2 .gt. n) go to 145 | ||||
| c | ||||
| do 140 i = l2, n | ||||
| 140 d(i) = d(i) - h | ||||
| c | ||||
| 145 f = f + h | ||||
| c .......... ql transformation .......... | ||||
| p = d(m) | ||||
| c = 1.0e0 | ||||
| c2 = c | ||||
| el1 = e(l1) | ||||
| s = 0.0e0 | ||||
| mml = m - l | ||||
| c .......... for i=m-1 step -1 until l do -- .......... | ||||
| do 200 ii = 1, mml | ||||
| c3 = c2 | ||||
| c2 = c | ||||
| s2 = s | ||||
| i = m - ii | ||||
| g = c * e(i) | ||||
| h = c * p | ||||
| r = pythag(p,e(i)) | ||||
| e(i+1) = s * r | ||||
| s = e(i) / r | ||||
| c = p / r | ||||
| p = c * d(i) - s * g | ||||
| d(i+1) = h + s * (c * g + s * d(i)) | ||||
| 200 continue | ||||
| c | ||||
| p = -s * s2 * c3 * el1 * e(l) / dl1 | ||||
| e(l) = s * p | ||||
| d(l) = c * p | ||||
| tst2 = tst1 + abs(e(l)) | ||||
| if (tst2 .gt. tst1) go to 130 | ||||
| 210 p = d(l) + f | ||||
| c .......... order eigenvalues .......... | ||||
| 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 | ||||
