tql2.f
170 lines
| 4.8 KiB
| text/x-fortran
|
FortranFixedLexer
| r1601 | subroutine tql2(nm,n,d,e,z,ierr) | |||
| c | ||||
| integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr | ||||
| real d(n),e(n),z(nm,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 tql2, | ||||
| 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 and eigenvectors | ||||
| c of a symmetric tridiagonal matrix by the ql method. | ||||
| c the eigenvectors of a full symmetric matrix can also | ||||
| c be found if tred2 has been used to reduce this | ||||
| c full matrix to tridiagonal form. | ||||
| c | ||||
| c on input | ||||
| c | ||||
| c nm must be set to the row dimension of two-dimensional | ||||
| c array parameters as declared in the calling program | ||||
| c dimension statement. | ||||
| 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 z contains the transformation matrix produced in the | ||||
| c reduction by tred2, if performed. if the eigenvectors | ||||
| c of the tridiagonal matrix are desired, z must contain | ||||
| c the identity matrix. | ||||
| c | ||||
| c on output | ||||
| c | ||||
| c d contains the eigenvalues in ascending order. if an | ||||
| c error exit is made, the eigenvalues are correct but | ||||
| c unordered for indices 1,2,...,ierr-1. | ||||
| c | ||||
| c e has been destroyed. | ||||
| c | ||||
| c z contains orthonormal eigenvectors of the symmetric | ||||
| c tridiagonal (or full) matrix. if an error exit is made, | ||||
| c z contains the eigenvectors associated with the stored | ||||
| c eigenvalues. | ||||
| 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 240 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 220 | ||||
| 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)) | ||||
| c .......... form vector .......... | ||||
| do 180 k = 1, n | ||||
| h = z(k,i+1) | ||||
| z(k,i+1) = s * z(k,i) + c * h | ||||
| z(k,i) = c * z(k,i) - s * h | ||||
| 180 continue | ||||
| c | ||||
| 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 | ||||
| 220 d(l) = d(l) + f | ||||
| 240 continue | ||||
| c .......... order eigenvalues and eigenvectors .......... | ||||
| do 300 ii = 2, n | ||||
| i = ii - 1 | ||||
| k = i | ||||
| p = d(i) | ||||
| c | ||||
| do 260 j = ii, n | ||||
| if (d(j) .ge. p) go to 260 | ||||
| k = j | ||||
| p = d(j) | ||||
| 260 continue | ||||
| c | ||||
| if (k .eq. i) go to 300 | ||||
| d(k) = d(i) | ||||
| d(i) = p | ||||
| c | ||||
| do 280 j = 1, n | ||||
| p = z(j,i) | ||||
| z(j,i) = z(j,k) | ||||
| z(j,k) = p | ||||
| 280 continue | ||||
| c | ||||
| 300 continue | ||||
| c | ||||
| go to 1001 | ||||
| c .......... set error -- no convergence to an | ||||
| c eigenvalue after 30 iterations .......... | ||||
| 1000 ierr = l | ||||
| 1001 return | ||||
| end | ||||
