website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/slatec/tqlrat.f

   1 *DECK TQLRAT
   2       SUBROUTINE TQLRAT (N, D, E2, IERR)
   3 C***BEGIN PROLOGUE  TQLRAT
   4 C***PURPOSE  Compute the eigenvalues of symmetric tridiagonal matrix
   5 C            using a rational variant of the QL method.
   6 C***LIBRARY   SLATEC (EISPACK)
   7 C***CATEGORY  D4A5, D4C2A
   8 C***TYPE      SINGLE PRECISION (TQLRAT-S)
   9 C***KEYWORDS  EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX, EISPACK,
  10 C             QL METHOD
  11 C***AUTHOR  Smith, B. T., et al.
  12 C***DESCRIPTION
  13 C
  14 C     This subroutine is a translation of the ALGOL procedure TQLRAT.
  15 C
  16 C     This subroutine finds the eigenvalues of a SYMMETRIC
  17 C     TRIDIAGONAL matrix by the rational QL method.
  18 C
  19 C     On Input
  20 C
  21 C        N is the order of the matrix.  N is an INTEGER variable.
  22 C
  23 C        D contains the diagonal elements of the symmetric tridiagonal
  24 C          matrix.  D is a one-dimensional REAL array, dimensioned D(N).
  25 C
  26 C        E2 contains the squares of the subdiagonal elements of the
  27 C          symmetric tridiagonal matrix in its last N-1 positions.
  28 C          E2(1) is arbitrary.  E2 is a one-dimensional REAL array,
  29 C          dimensioned E2(N).
  30 C
  31 C      On Output
  32 C
  33 C        D contains the eigenvalues in ascending order.  If an
  34 C          error exit is made, the eigenvalues are correct and
  35 C          ordered for indices 1, 2, ..., IERR-1, but may not be
  36 C          the smallest eigenvalues.
  37 C
  38 C        E2 has been destroyed.
  39 C
  40 C        IERR is an INTEGER flag set to
  41 C          Zero       for normal return,
  42 C          J          if the J-th eigenvalue has not been
  43 C                     determined after 30 iterations.
  44 C
  45 C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
  46 C
  47 C     Questions and comments should be directed to B. S. Garbow,
  48 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  49 C     ------------------------------------------------------------------
  50 C
  51 C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  52 C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  53 C                 system Routines - EISPACK Guide, Springer-Verlag,
  54 C                 1976.
  55 C               C. H. Reinsch, Eigenvalues of a real, symmetric, tri-
  56 C                 diagonal matrix, Algorithm 464, Communications of the
  57 C                 ACM 16, 11 (November 1973), pp. 689.
  58 C***ROUTINES CALLED  PYTHAG, R1MACH
  59 C***REVISION HISTORY  (YYMMDD)
  60 C   760101  DATE WRITTEN
  61 C   890831  Modified array declarations.  (WRB)
  62 C   890831  REVISION DATE from Version 3.2
  63 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  64 C   920501  Reformatted the REFERENCES section.  (WRB)
  65 C***END PROLOGUE  TQLRAT
  66 C
  67       INTEGER I,J,L,M,N,II,L1,MML,IERR
  68       REAL D(*),E2(*)
  69       REAL B,C,F,G,H,P,R,S,MACHEP
  70       REAL PYTHAG
  71       LOGICAL FIRST
  72 C
  73       SAVE FIRST, MACHEP
  74       DATA FIRST /.TRUE./
  75 C***FIRST EXECUTABLE STATEMENT  TQLRAT
  76       IF (FIRST) THEN
  77          MACHEP = R1MACH(4)
  78       ENDIF
  79       FIRST = .FALSE.
  80 C
  81       IERR = 0
  82       IF (N .EQ. 1) GO TO 1001
  83 C
  84       DO 100 I = 2, N
  85   100 E2(I-1) = E2(I)
  86 C
  87       F = 0.0E0
  88       B = 0.0E0
  89       E2(N) = 0.0E0
  90 C
  91       DO 290 L = 1, N
  92          J = 0
  93          H = MACHEP * (ABS(D(L)) + SQRT(E2(L)))
  94          IF (B .GT. H) GO TO 105
  95          B = H
  96          C = B * B
  97 C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
  98   105    DO 110 M = L, N
  99             IF (E2(M) .LE. C) GO TO 120
 100 C     .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
 101 C                THROUGH THE BOTTOM OF THE LOOP ..........
 102   110    CONTINUE
 103 C
 104   120    IF (M .EQ. L) GO TO 210
 105   130    IF (J .EQ. 30) GO TO 1000
 106          J = J + 1
 107 C     .......... FORM SHIFT ..........
 108          L1 = L + 1
 109          S = SQRT(E2(L))
 110          G = D(L)
 111          P = (D(L1) - G) / (2.0E0 * S)
 112          R = PYTHAG(P,1.0E0)
 113          D(L) = S / (P + SIGN(R,P))
 114          H = G - D(L)
 115 C
 116          DO 140 I = L1, N
 117   140    D(I) = D(I) - H
 118 C
 119          F = F + H
 120 C     .......... RATIONAL QL TRANSFORMATION ..........
 121          G = D(M)
 122          IF (G .EQ. 0.0E0) G = B
 123          H = G
 124          S = 0.0E0
 125          MML = M - L
 126 C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
 127          DO 200 II = 1, MML
 128             I = M - II
 129             P = G * H
 130             R = P + E2(I)
 131             E2(I+1) = S * R
 132             S = E2(I) / R
 133             D(I+1) = H + S * (H + D(I))
 134             G = D(I) - E2(I) / G
 135             IF (G .EQ. 0.0E0) G = B
 136             H = G * P / R
 137   200    CONTINUE
 138 C
 139          E2(L) = S * G
 140          D(L) = H
 141 C     .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ..........
 142          IF (H .EQ. 0.0E0) GO TO 210
 143          IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
 144          E2(L) = H * E2(L)
 145          IF (E2(L) .NE. 0.0E0) GO TO 130
 146   210    P = D(L) + F
 147 C     .......... ORDER EIGENVALUES ..........
 148          IF (L .EQ. 1) GO TO 250
 149 C     .......... FOR I=L STEP -1 UNTIL 2 DO -- ..........
 150          DO 230 II = 2, L
 151             I = L + 2 - II
 152             IF (P .GE. D(I-1)) GO TO 270
 153             D(I) = D(I-1)
 154   230    CONTINUE
 155 C
 156   250    I = 1
 157   270    D(I) = P
 158   290 CONTINUE
 159 C
 160       GO TO 1001
 161 C     .......... SET ERROR -- NO CONVERGENCE TO AN
 162 C                EIGENVALUE AFTER 30 ITERATIONS ..........
 163  1000 IERR = L
 164  1001 RETURN
 165       END