imtqlv.f

C
C     ------------------------------------------------------------------
C
      SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
C
      INTEGER I,J,K,L,M,N,II,MML,TAG,IERR
      REAL*8 D(N),E(N),E2(N),W(N),RV1(N)
      REAL*8 B,C,F,G,P,R,S
      REAL*8 MACHEP
C	  PARAMETER (MACHEP = $3410000000000000)
	  PARAMETER (MACHEP = 16.0D0**(-13))
      REAL*8 DSQRT,DABS,DSIGN
      INTEGER IND(N)
C
C     THIS SUBROUTINE IS A VARIANT OF  IMTQL1  WHICH IS A TRANSLATION OF
C     ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
C     WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
C     MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
C     THEIR CORRESPONDING SUBMATRIX INDICES.
C
C     ON INPUT
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        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C          E2(1) IS ARBITRARY.
C
C     ON OUTPUT
C
C        D AND E ARE UNALTERED;
C
C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
C          E2(1) IS ALSO SET TO ZERO;
C
C        W 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        IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
C          CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
C          BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
C          2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.;
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        RV1 IS A TEMPORARY STORAGE ARRAY.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C                MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C                MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC
C                ON S360
C
      IERR = 0
      K = 0
      TAG = 0
C
      DO 10 I = 1, N
      W(I) = D(I)
      IF (I .NE. 1) RV1(I-1) = E(I)
   10 CONTINUE
C
      E2(1) = 0.0D0
      RV1(N) = 0.0D0
C
      DO 140 L = 1, N
      J = 0
C                LOOK FOR SMALL SUB-DIAGONAL ELEMENT
   20 DO 30 M = L, N
      IF (M .EQ. N) GO TO 40
      IF (DABS(RV1(M)) .LE. MACHEP * (DABS(W(M)) + DABS(W(M+1)))) GO TO
     X     40
C                GUARD AGAINST UNDERFLOWED ELEMENT OF E2
      IF (E2(M+1) .EQ. 0.0D0) GO TO 50
   30 CONTINUE
C
   40 IF (M .LE. K) GO TO 60
      IF (M .NE. N) E2(M+1) = 0.0D0
   50 K = M
      TAG = TAG + 1
   60 P = W(L)
      IF (M .EQ. L) GO TO 100
      IF (J .EQ. 30) GO TO 150
      J = J + 1
C                FORM SHIFT
      G = (W(L+1) - P) / (2.0D0 * RV1(L))
      R = DSQRT(G*G+1.0D0)
      G = W(M) - P + RV1(L) / (G + DSIGN(R,G))
      S = 1.0D0
      C = 1.0D0
      P = 0.0D0
      MML = M - L
C                FOR I=M-1 STEP -1 UNTIL L DO --
      DO 90 II = 1, MML
      I = M - II
      F = S * RV1(I)
      B = C * RV1(I)
      IF (DABS(F) .LT. DABS(G)) GO TO 70
      C = G / F
      R = DSQRT(C*C+1.0D0)
      RV1(I+1) = F * R
      S = 1.0D0 / R
      C = C * S
      GO TO 80
   70 S = F / G
      R = DSQRT(S*S+1.0D0)
      RV1(I+1) = G * R
      C = 1.0D0 / R
      S = S * C
   80 G = W(I+1) - P
      R = (W(I) - G) * S + 2.0D0 * C * B
      P = S * R
      W(I+1) = G + P
      G = C * R - B
   90 CONTINUE
C
      W(L) = W(L) - P
      RV1(L) = G
      RV1(M) = 0.0D0
      GO TO 20
C                ORDER EIGENVALUES
  100 IF (L .EQ. 1) GO TO 120
C                FOR I=L STEP -1 UNTIL 2 DO --
      DO 110 II = 2, L
      I = L + 2 - II
      IF (P .GE. W(I-1)) GO TO 130
      W(I) = W(I-1)
      IND(I) = IND(I-1)
  110 CONTINUE
C
  120 I = 1
  130 W(I) = P
      IND(I) = TAG
  140 CONTINUE
C
      GO TO 160
C                SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS
  150 IERR = L
  160 RETURN
C                LAST CARD OF IMTQLV
      END