X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=website%2Farchive%2Fbinaries%2Fmac%2Fsrc%2Fdisembl%2FTisean_3.0.1%2Fsource_f%2Fslatec%2Ftql2.f;fp=website%2Farchive%2Fbinaries%2Fmac%2Fsrc%2Fdisembl%2FTisean_3.0.1%2Fsource_f%2Fslatec%2Ftql2.f;h=40d99385ffa473eab8024a34b7b5fa4862431f91;hb=dbde3fb6f00b9bb770343631a517c0e599db8528;hp=0000000000000000000000000000000000000000;hpb=85f830bbd51a7277994bd4233141016304e210c9;p=jabaws.git diff --git a/website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tql2.f b/website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tql2.f new file mode 100644 index 0000000..40d9938 --- /dev/null +++ b/website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tql2.f @@ -0,0 +1,203 @@ +*DECK TQL2 + SUBROUTINE TQL2 (NM, N, D, E, Z, IERR) +C***BEGIN PROLOGUE TQL2 +C***PURPOSE Compute the eigenvalues and eigenvectors of symmetric +C tridiagonal matrix. +C***LIBRARY SLATEC (EISPACK) +C***CATEGORY D4A5, D4C2A +C***TYPE SINGLE PRECISION (TQL2-S) +C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK +C***AUTHOR Smith, B. T., et al. +C***DESCRIPTION +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 the two-dimensional +C array parameter, Z, as declared in the calling program +C dimension statement. NM is an INTEGER variable. +C +C N is the order of the matrix. N is an INTEGER variable. +C N must be less than or equal to NM. +C +C D contains the diagonal elements of the symmetric tridiagonal +C matrix. D is a one-dimensional REAL array, dimensioned D(N). +C +C E contains the subdiagonal elements of the symmetric +C tridiagonal matrix in its last N-1 positions. E(1) is +C arbitrary. E is a one-dimensional REAL array, dimensioned +C E(N). +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. Z is a two-dimensional REAL array, +C dimensioned Z(NM,N). +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 an INTEGER flag 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(A,B) for sqrt(A**2 + B**2). +C +C Questions and comments should be directed to B. S. Garbow, +C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY +C ------------------------------------------------------------------ +C +C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, +C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- +C system Routines - EISPACK Guide, Springer-Verlag, +C 1976. +C***ROUTINES CALLED PYTHAG +C***REVISION HISTORY (YYMMDD) +C 760101 DATE WRITTEN +C 890831 Modified array declarations. (WRB) +C 890831 REVISION DATE from Version 3.2 +C 891214 Prologue converted to Version 4.0 format. (BAB) +C 920501 Reformatted the REFERENCES section. (WRB) +C***END PROLOGUE TQL2 +C + INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR + REAL D(*),E(*),Z(NM,*) + REAL B,C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2 + REAL PYTHAG +C +C***FIRST EXECUTABLE STATEMENT TQL2 + 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 + B = 0.0E0 + E(N) = 0.0E0 +C + DO 240 L = 1, N + J = 0 + H = ABS(D(L)) + ABS(E(L)) + IF (B .LT. H) B = H +C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... + DO 110 M = L, N + IF (B + ABS(E(M)) .EQ. B) 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 + IF (ABS(P) .LT. ABS(E(I))) GO TO 150 + C = E(I) / P + R = SQRT(C*C+1.0E0) + E(I+1) = S * P * R + S = C / R + C = 1.0E0 / R + GO TO 160 + 150 C = P / E(I) + R = SQRT(C*C+1.0E0) + E(I+1) = S * E(I) * R + S = 1.0E0 / R + C = C * S + 160 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 + IF (B + ABS(E(L)) .GT. B) 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