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

   1 *DECK TQL2
   2       SUBROUTINE TQL2 (NM, N, D, E, Z, IERR)
   3 C***BEGIN PROLOGUE  TQL2
   4 C***PURPOSE  Compute the eigenvalues and eigenvectors of symmetric
   5 C            tridiagonal matrix.
   6 C***LIBRARY   SLATEC (EISPACK)
   7 C***CATEGORY  D4A5, D4C2A
   8 C***TYPE      SINGLE PRECISION (TQL2-S)
   9 C***KEYWORDS  EIGENVALUES, EIGENVECTORS, EISPACK
  10 C***AUTHOR  Smith, B. T., et al.
  11 C***DESCRIPTION
  12 C
  13 C     This subroutine is a translation of the ALGOL procedure TQL2,
  14 C     NUM. MATH. 11, 293-306(1968) by Bowdler, Martin, Reinsch, and
  15 C     Wilkinson.
  16 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
  17 C
  18 C     This subroutine finds the eigenvalues and eigenvectors
  19 C     of a SYMMETRIC TRIDIAGONAL matrix by the QL method.
  20 C     The eigenvectors of a FULL SYMMETRIC matrix can also
  21 C     be found if  TRED2  has been used to reduce this
  22 C     full matrix to tridiagonal form.
  23 C
  24 C     On Input
  25 C
  26 C        NM must be set to the row dimension of the two-dimensional
  27 C          array parameter, Z, as declared in the calling program
  28 C          dimension statement.  NM is an INTEGER variable.
  29 C
  30 C        N is the order of the matrix.  N is an INTEGER variable.
  31 C          N must be less than or equal to NM.
  32 C
  33 C        D contains the diagonal elements of the symmetric tridiagonal
  34 C          matrix.  D is a one-dimensional REAL array, dimensioned D(N).
  35 C
  36 C        E contains the subdiagonal elements of the symmetric
  37 C          tridiagonal matrix in its last N-1 positions.  E(1) is
  38 C          arbitrary.  E is a one-dimensional REAL array, dimensioned
  39 C          E(N).
  40 C
  41 C        Z contains the transformation matrix produced in the
  42 C          reduction by  TRED2, if performed.  If the eigenvectors
  43 C          of the tridiagonal matrix are desired, Z must contain
  44 C          the identity matrix.  Z is a two-dimensional REAL array,
  45 C          dimensioned Z(NM,N).
  46 C
  47 C      On Output
  48 C
  49 C        D contains the eigenvalues in ascending order.  If an
  50 C          error exit is made, the eigenvalues are correct but
  51 C          unordered for indices 1, 2, ..., IERR-1.
  52 C
  53 C        E has been destroyed.
  54 C
  55 C        Z contains orthonormal eigenvectors of the symmetric
  56 C          tridiagonal (or full) matrix.  If an error exit is made,
  57 C          Z contains the eigenvectors associated with the stored
  58 C          eigenvalues.
  59 C
  60 C        IERR is an INTEGER flag set to
  61 C          Zero       for normal return,
  62 C          J          if the J-th eigenvalue has not been
  63 C                     determined after 30 iterations.
  64 C
  65 C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
  66 C
  67 C     Questions and comments should be directed to B. S. Garbow,
  68 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  69 C     ------------------------------------------------------------------
  70 C
  71 C***REFERENCES  B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow,
  72 C                 Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen-
  73 C                 system Routines - EISPACK Guide, Springer-Verlag,
  74 C                 1976.
  75 C***ROUTINES CALLED  PYTHAG
  76 C***REVISION HISTORY  (YYMMDD)
  77 C   760101  DATE WRITTEN
  78 C   890831  Modified array declarations.  (WRB)
  79 C   890831  REVISION DATE from Version 3.2
  80 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  81 C   920501  Reformatted the REFERENCES section.  (WRB)
  82 C***END PROLOGUE  TQL2
  83 C
  84       INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR
  85       REAL D(*),E(*),Z(NM,*)
  86       REAL B,C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2
  87       REAL PYTHAG
  88 C
  89 C***FIRST EXECUTABLE STATEMENT  TQL2
  90       IERR = 0
  91       IF (N .EQ. 1) GO TO 1001
  92 C
  93       DO 100 I = 2, N
  94   100 E(I-1) = E(I)
  95 C
  96       F = 0.0E0
  97       B = 0.0E0
  98       E(N) = 0.0E0
  99 C
 100       DO 240 L = 1, N
 101          J = 0
 102          H = ABS(D(L)) + ABS(E(L))
 103          IF (B .LT. H) B = H
 104 C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
 105          DO 110 M = L, N
 106             IF (B + ABS(E(M)) .EQ. B) GO TO 120
 107 C     .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
 108 C                THROUGH THE BOTTOM OF THE LOOP ..........
 109   110    CONTINUE
 110 C
 111   120    IF (M .EQ. L) GO TO 220
 112   130    IF (J .EQ. 30) GO TO 1000
 113          J = J + 1
 114 C     .......... FORM SHIFT ..........
 115          L1 = L + 1
 116          L2 = L1 + 1
 117          G = D(L)
 118          P = (D(L1) - G) / (2.0E0 * E(L))
 119          R = PYTHAG(P,1.0E0)
 120          D(L) = E(L) / (P + SIGN(R,P))
 121          D(L1) = E(L) * (P + SIGN(R,P))
 122          DL1 = D(L1)
 123          H = G - D(L)
 124          IF (L2 .GT. N) GO TO 145
 125 C
 126          DO 140 I = L2, N
 127   140    D(I) = D(I) - H
 128 C
 129   145    F = F + H
 130 C     .......... QL TRANSFORMATION ..........
 131          P = D(M)
 132          C = 1.0E0
 133          C2 = C
 134          EL1 = E(L1)
 135          S = 0.0E0
 136          MML = M - L
 137 C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
 138          DO 200 II = 1, MML
 139             C3 = C2
 140             C2 = C
 141             S2 = S
 142             I = M - II
 143             G = C * E(I)
 144             H = C * P
 145             IF (ABS(P) .LT. ABS(E(I))) GO TO 150
 146             C = E(I) / P
 147             R = SQRT(C*C+1.0E0)
 148             E(I+1) = S * P * R
 149             S = C / R
 150             C = 1.0E0 / R
 151             GO TO 160
 152   150       C = P / E(I)
 153             R = SQRT(C*C+1.0E0)
 154             E(I+1) = S * E(I) * R
 155             S = 1.0E0 / R
 156             C = C * S
 157   160       P = C * D(I) - S * G
 158             D(I+1) = H + S * (C * G + S * D(I))
 159 C     .......... FORM VECTOR ..........
 160             DO 180 K = 1, N
 161                H = Z(K,I+1)
 162                Z(K,I+1) = S * Z(K,I) + C * H
 163                Z(K,I) = C * Z(K,I) - S * H
 164   180       CONTINUE
 165 C
 166   200    CONTINUE
 167 C
 168          P = -S * S2 * C3 * EL1 * E(L) / DL1
 169          E(L) = S * P
 170          D(L) = C * P
 171          IF (B + ABS(E(L)) .GT. B) GO TO 130
 172   220    D(L) = D(L) + F
 173   240 CONTINUE
 174 C     .......... ORDER EIGENVALUES AND EIGENVECTORS ..........
 175       DO 300 II = 2, N
 176          I = II - 1
 177          K = I
 178          P = D(I)
 179 C
 180          DO 260 J = II, N
 181             IF (D(J) .GE. P) GO TO 260
 182             K = J
 183             P = D(J)
 184   260    CONTINUE
 185 C
 186          IF (K .EQ. I) GO TO 300
 187          D(K) = D(I)
 188          D(I) = P
 189 C
 190          DO 280 J = 1, N
 191             P = Z(J,I)
 192             Z(J,I) = Z(J,K)
 193             Z(J,K) = P
 194   280    CONTINUE
 195 C
 196   300 CONTINUE
 197 C
 198       GO TO 1001
 199 C     .......... SET ERROR -- NO CONVERGENCE TO AN
 200 C                EIGENVALUE AFTER 30 ITERATIONS ..........
 201  1000 IERR = L
 202  1001 RETURN
 203       END