Mac binaries

[jabaws.git] / website / archive / binaries / mac / src / disembl / Tisean_3.0.1 / source_f / slatec / tqlrat.f

diff --git a/website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tqlrat.f b/website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tqlrat.f
new file mode 100644 (file)
index 0000000..8cb7b9c
--- /dev/null
+++ b/website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tqlrat.f
@@ -0,0 +1,165 @@
+*DECK TQLRAT
+      SUBROUTINE TQLRAT (N, D, E2, IERR)
+C***BEGIN PROLOGUE  TQLRAT
+C***PURPOSE  Compute the eigenvalues of symmetric tridiagonal matrix
+C            using a rational variant of the QL method.
+C***LIBRARY   SLATEC (EISPACK)
+C***CATEGORY  D4A5, D4C2A
+C***TYPE      SINGLE PRECISION (TQLRAT-S)
+C***KEYWORDS  EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX, EISPACK,
+C             QL METHOD
+C***AUTHOR  Smith, B. T., et al.
+C***DESCRIPTION
+C
+C     This subroutine is a translation of the ALGOL procedure TQLRAT.
+C
+C     This subroutine finds the eigenvalues of a SYMMETRIC
+C     TRIDIAGONAL matrix by the rational QL method.
+C
+C     On Input
+C
+C        N is the order of the matrix.  N is an INTEGER variable.
+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        E2 contains the squares of the subdiagonal elements of the
+C          symmetric tridiagonal matrix in its last N-1 positions.
+C          E2(1) is arbitrary.  E2 is a one-dimensional REAL array,
+C          dimensioned E2(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 and
+C          ordered for indices 1, 2, ..., IERR-1, but may not be
+C          the smallest eigenvalues.
+C
+C        E2 has been destroyed.
+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               C. H. Reinsch, Eigenvalues of a real, symmetric, tri-
+C                 diagonal matrix, Algorithm 464, Communications of the
+C                 ACM 16, 11 (November 1973), pp. 689.
+C***ROUTINES CALLED  PYTHAG, R1MACH
+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  TQLRAT
+C
+      INTEGER I,J,L,M,N,II,L1,MML,IERR
+      REAL D(*),E2(*)
+      REAL B,C,F,G,H,P,R,S,MACHEP
+      REAL PYTHAG
+      LOGICAL FIRST
+C
+      SAVE FIRST, MACHEP
+      DATA FIRST /.TRUE./
+C***FIRST EXECUTABLE STATEMENT  TQLRAT
+      IF (FIRST) THEN
+         MACHEP = R1MACH(4)
+      ENDIF
+      FIRST = .FALSE.
+C
+      IERR = 0
+      IF (N .EQ. 1) GO TO 1001
+C
+      DO 100 I = 2, N
+  100 E2(I-1) = E2(I)
+C
+      F = 0.0E0
+      B = 0.0E0
+      E2(N) = 0.0E0
+C
+      DO 290 L = 1, N
+         J = 0
+         H = MACHEP * (ABS(D(L)) + SQRT(E2(L)))
+         IF (B .GT. H) GO TO 105
+         B = H
+         C = B * B
+C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
+  105    DO 110 M = L, N
+            IF (E2(M) .LE. C) GO TO 120
+C     .......... E2(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
+         S = SQRT(E2(L))
+         G = D(L)
+         P = (D(L1) - G) / (2.0E0 * S)
+         R = PYTHAG(P,1.0E0)
+         D(L) = S / (P + SIGN(R,P))
+         H = G - D(L)
+C
+         DO 140 I = L1, N
+  140    D(I) = D(I) - H
+C
+         F = F + H
+C     .......... RATIONAL QL TRANSFORMATION ..........
+         G = D(M)
+         IF (G .EQ. 0.0E0) G = B
+         H = G
+         S = 0.0E0
+         MML = M - L
+C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
+         DO 200 II = 1, MML
+            I = M - II
+            P = G * H
+            R = P + E2(I)
+            E2(I+1) = S * R
+            S = E2(I) / R
+            D(I+1) = H + S * (H + D(I))
+            G = D(I) - E2(I) / G
+            IF (G .EQ. 0.0E0) G = B
+            H = G * P / R
+  200    CONTINUE
+C
+         E2(L) = S * G
+         D(L) = H
+C     .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
+         IF (H .EQ. 0.0E0) GO TO 210
+         IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
+         E2(L) = H * E2(L)
+         IF (E2(L) .NE. 0.0E0) 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