C      FGHEVEN
C      F. Gogtas, G.G. Balint-Kurti and C.C. Marston,
C      QCPE Program No. 647, (1993).
C
C     ***************************************************************
C     This program solves one dimensional Schrodinger equation for
C     bound state eigenvalues and eigenfunctions corresponding to a
C     potential V(x).
C
C     The Fourier Grid Hamiltonian method is used in the program.
C     This method is fully described in:
C     C.C. Marston and G.G.Balint-Kurti, J.Chem. Phys.,91,3571(1989).
C     and
C     G.G. Balint-Kurti, R.N. Dixon and C.C. Marston, Internat. Rev. 
C     Phys. Chem.,11, 317 (1992).
C
C     The program uses an even number of grid points. The Hamiltonian
C     matrix elements are calculated using analytic formula described
C     in the second of the above references.
C
C     The analytical Hamiltonian expression given in the reference
C     contains a small error.  The formula should read: 
C
C                H(i,j) = {(h**2)/(4*m*(L**2)} *
C                          [(N-1)(N-2)/6 + (N/2)] + V(Xi),   if i=j
C
C                H(i,j) = {[(-1)**(i-j)] / m } *
C                          { h/[2*L*sin(pi*(i-j)/N)]}**2 ,     if i#j
C
C     The eigenvalues of the Hamiltonian matrix which lie below the
C     asymptotic (large x) value of V(x) correspond to the bound state
C     energies. The corresponding eigenvectors are the representation
C     of the bound state wavefunctions on a regular one dimensional grid.
C
C     The user must supply a self contained subroutine VSUB(X,V)
C     which is called with a value of X and returns the value of the
C     potential in V.
C
C     As supplied the program uses a Morse potential with parameters chosen 
C     to represent the HF molecule. 
C
C     The Program should work on any computer
C     Programming language used                    : Fortran77
C     Memory required to execute with typical data : 1 Mbyte
C
C     To run the program on a SUN workstation (or any Unix based system)
C     First type:
C        f77 -o fgheven fgheven.f
C     to create an executable program (fgheven).
C     Then type:
C        fgheven > fgheven.out &   
C     to run the program.  
C     
C     The dimensions of the arrays generally depend on the number of grid
C     points used.  This number NX is set in a PARAMETER statement.
C     The following arrays represent  :
C
C        XA  : X-coordinates i.e. position on a grid.
C
C        AR  : The Hamiltonian Matrix
C         
C        WCH : The eigenvalues for respective energy levels.
C
C        ZR  : The eigenvectors (X,Y); where
C                                      X : Wavefunction.
C                                      Y : Energy level.       
C    The folowing data constans represent :
C        ZMA, ZMB : Mass of atoms.
C        R0       : Equilibrium bond separation
C        NPRIN    : 0 or 1 for eigenvalue or eigenvector respectively
C        NWRIT    : Number of eigenvales/eigenvectors to be printed
C        RMIN     : Starting point of grid
C        RMAX     : End point of grid
C        ZL       : Grid length
C        DX       : Grid spacings
C
C       All quantities in au (unless otherwise stated)
C     ****************************************************************
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (NX=1800)
C
C.....Declare arrays.
      DIMENSION XA(NX),AR(NX,NX),WCH(NX),ZR(NX,NX) 
      DIMENSION FV1(NX),FV2(NX)
      DIMENSION POT(NX)
C
C.....Variable data input.
C     DATA ZMA/1836.9822D0/
C     DATA ZMB/34629.61319D0/
C     DATA R0/1.7329D0/
C     DATA DHART/0.2250073497D0/
C     DATA BETA/1.1741D0/
      DATA PI/3.141592653589793D0/
      DATA NPRIN/0/
      DATA NWRIT/20/
      INTEGER :: L
      LOGICAL :: MORLJ  

C     NAMELIST DATA
      NAMELIST /TTDATA/ ZMA, ZMB, R0, BETA, DHART, NPRIN, NWRIT, L,
     &RMAX, RMIN, MORLJ, C6, C8, C10, C12
      READ(5,NML=TTDATA)

C
C....Test that NX is even
      ITEST=MOD(NX,2)
      IF(ITEST.NE.0) THEN
        WRITE(6,*)' **** NX MUST BE EVEN-FATAL ERROR  ****'
        STOP
      END IF
C
C.....Set up grid
      WRITE(6,*)'Grid paremeters:'
      WRITE(6,*)'   Number of grid poins    = ',NX
C      RMIN=R0*0.2D0
C      RMAX=R0*3.0D0
      ZL=(RMAX-RMIN)
      WRITE(6,*)'   Grid length             = ',ZL
      DX=ZL/DFLOAT(NX)
      WRITE(6,*)'   Grid spacings           = ',DX
C.....Print mass of atoms A and B
      WRITE(6,*)'Mass of atoms:'
      WRITE(6,*)'   Atomic mass of atom A      = ',ZMA
      WRITE(6,*)'   Atomic mass of atom B      = ',ZMB
      WRITE(6,*)'   R0                         = ',R0
      WRITE(6,*)'   BETA                       = ',BETA
      WRITE(6,*)'   DHART                      = ',DHART
C.....Compute reduced mass.
      ZMU=ZMA*ZMB/(ZMA+ZMB)
      WRITE(6,*)'   Reduced mass of AB         = ',ZMU
C
C.....Compute contants:
      PSQ=PI*PI
      CONST1=PSQ/(ZMU*(ZL**2))
      NFACT1=(NX-1)*(NX-2)
      NFACT2=(NX-2)/2
      CONST2=CONST1*(DFLOAT(NFACT1)/6.0D0 +1.0D0+DFLOAT(NFACT2))
      DARG=PI/DFLOAT(NX)
      VMIN=0.0D10
C
C.....Now compute Hamiltonian matrix:
      X=RMIN
      DO 003 I=1,NX
         XA(I)=X
       DO 004 J=1,I
          IJD=(I-J)
        IF(IJD.EQ.0)THEN
          AR(I,J)=CONST2
        ELSE
          RATIO=1.0D0/DSIN(DARG*DFLOAT(IJD))
          AR(I,J)=((-1)**IJD)*CONST1*(RATIO**2)
        END IF
 004  CONTINUE
C.....Find the potential value at x
      IF(MORLJ) THEN
       CALL VSUB_M(X,VV,R0,BETA,DHART,L,ZMU)
      ELSE
       CALL VSUB_LJ(X,VV,R0,C6,C8,C10,C12,L,ZMU)
      ENDIF
      POT(I) = VV
      IF(VMIN .gt. VV) VMIN = VV
C.....Add the potential value when the kronecker delta function
C.....equals to one, i.e. when I and J are equal
        AR(I,I)=AR(I,I)+VV
        X=X+DX
 003  CONTINUE
C.....Now fill out Hamiltonian matrix.
      DO 006 I=1,NX
       DO 007 J=1,I
        AR(J,I)=AR(I,J)
007    CONTINUE
006   CONTINUE
C.....Now call eigenvalue solver.
      call flush(6)
      CALL RS(NX,NX,AR,WCH,NPRIN,ZR,FV1,FV2,IERR)
C.....Print out eigenvalues and eigenvectors.
      WRITE(6,12)
12    FORMAT(/)
      WRITE(6,*)' The first',NWRIT,' energy levels for AB molecule '
      DO 009 I=1,NWRIT  
       WRITE(6,*)' ENERGY LEVEL NO ',I,' EIGENVALUE= ',WCH(I)
       IF (WCH(I) .LT. 0.0D0 .AND. WCH(I) .GT. VMIN) THEN
        WRITE(10,*) I-1.0, WCH(I)
       ENDIF
 009  CONTINUE
      WRITE(6,12)
      WRITE(6,*)' THE CORRESPONDING EIGENFUNCTIONS ARE:'
      DO 010 I=1,NWRIT  
       WRITE(6,*)' ENERGY LEVEL NO ',I,' EIGENVALUE= ',WCH(I)
       IF (NPRIN.EQ.1) THEN
       DO 011 J=1,NX
         WRITE(6,*)' R = ',XA(J),' WAVEFUNCTION=',ZR(J,I)
 011   CONTINUE
       ENDIF
 010  CONTINUE
      OPEN(unit=99,file="pot.dat")
      DO 013 I=1,NX
         WRITE(99,*)XA(I),POT(I)
 013  CONTINUE
      CLOSE(99)
      STOP
      END
C
C
      SUBROUTINE VSUB_M(X,V,R0,BETA,DHART,L,MU)
C     **************************************************************
C.....SUBROUTINE TO COMPUTE POTENTIAL ENERGY AT SEPERATION X (IN   
C.....ATOMIC UNITS). THE VALUE OF THE POTENTIAL IS RETURNED IN V,  
C.....AGAIN IN ATOMIC UNITS.THE POTENTIAL ENERGY IS ASSUMED TO     
C..... BE A MORSE CURVE:V(X) = DHART * (EXP(-BETA*(X-R0)) - 1)**2  
C..... THE FOLLOWING DATA CONSTANTS REPRESENT :                    
C.....    DHART : THE DISSOCIATION ENERGY FOR THE GIVEN MOLECULE   
C.....    BETA  : THE EXPONENTIAL PARAMETER                        
C.....    R0    : THE BOND LENGTH AT EQUILIBRIUM FOR GIVEN MOLECULE
C     **************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION :: MU
      INTEGER :: L
C      DATA DHART/0.2250073497D0/
C      DATA BETA/1.1741D0/
C      DATA R0/1.7329D0/
C.....Set morse curve parameters.
      ARG=X-R0
      ARG=-BETA*ARG
      EF=DEXP(ARG)
      EF=1.0D0-EF
      EF=EF*EF
      V=DHART*EF-DHART
      V = V + DBLE(L)*(DBLE(L)+1.0)/(2.0*MU*X**2)
      RETURN
      END
CC
      SUBROUTINE VSUB_LJ(X,V,R0,C6,C8,C10,C12,L,MU)
C     **************************************************************
C.....SUBROUTINE TO COMPUTE POTENTIAL ENERGY AT SEPERATION X (IN   
C.....ATOMIC UNITS). THE VALUE OF THE POTENTIAL IS RETURNED IN V,  
C.....AGAIN IN ATOMIC UNITS.THE POTENTIAL ENERGY IS ASSUMED TO     
C..... BE A MORSE CURVE:V(X) = DHART * (EXP(-BETA*(X-R0)) - 1)**2  
C..... THE FOLLOWING DATA CONSTANTS REPRESENT :                    
C.....    DHART : THE DISSOCIATION ENERGY FOR THE GIVEN MOLECULE   
C.....    BETA  : THE EXPONENTIAL PARAMETER                        
C.....    R0    : THE BOND LENGTH AT EQUILIBRIUM FOR GIVEN MOLECULE
C     **************************************************************
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION :: MU
      INTEGER :: L
C      DATA DHART/0.2250073497D0/
C      DATA BETA/1.1741D0/
C      DATA R0/1.7329D0/
C.....Set morse curve parameters.
      V =C12/X**12-C10/X**6-C8/X**6-C6/X**6
      V = V + DBLE(L)*(DBLE(L)+1.0)/(2.0*MU*X**2)
      RETURN
      END
C
C      
      SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
C     ****************************************************************
C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF                     
C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)         
C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
C     OF A REAL SYMMETRIC MATRIX.
C
C     ON INPUT : 
C
C        N,M  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 
C        DIMENSION STATEMENT,
C        N  IS THE ORDER OF THE MATRIX  A,
C        A  CONTAINS THE REAL SYMMETRIC MATRIX,
C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
C        ONLY EIGENVALUES ARE DESIRED,  OTHERWISE IT IS SET TO
C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
C
C     ON OUTPUT :
C
C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER,
C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO,
C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN
C        ERROR COMPLETION CODE DESCRIBED IN SECTION 2B OF THE
C        DOCUMENTATION.  THE NORMAL COMPLETION CODE IS ZERO,
C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C    *******************************************************************
      IF (N .LE. NM) GO TO 10
      IERR = 10 * N
      GO TO 50
C
   10 IF (MATZ .NE. 0) GO TO 20
C     ********** FIND EIGENVALUES ONLY **********
      CALL  TRED1(NM,N,A,W,FV1,FV2)
      CALL  TQLRAT(N,W,FV2,IERR)
      GO TO 50
C     ********** FIND BOTH EIGENVALUES AND EIGENVECTORS **********
   20 CALL  TRED2(NM,N,A,W,FV1,Z)
      CALL  TQL2(NM,N,W,FV1,Z,IERR)
   50 RETURN
C     ********** LAST CARD OF RS **********
      END
C
C
      SUBROUTINE TRED1(NM,N,A,D,E,E2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NM,N),D(N),E(N),E2(N)
C     ***************************************************************
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT :
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT :
C
C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
C          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
C          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
C
C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C 
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     ****************************************************************
C
      DO 100 I = 1, N
  100 D(I) = A(I,I)
C     ********** FOR I=N STEP -1 UNTIL 1 DO -- ********** 
      DO 300 II = 1, N 
         I = N + 1 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 1) GO TO 130
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** 
         DO 120 K = 1, L 
  120    SCALE = SCALE + DABS(A(I,K)) 
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = 0.0D0
         E2(I) = 0.0D0
         GO TO 290
C
  140    DO 150 K = 1, L
            A(I,K) = A(I,K) / SCALE 
            H = H + A(I,K) * A(I,K)
  150    CONTINUE
C
         E2(I) = SCALE * SCALE * H
         F = A(I,L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         A(I,L) = F - G
         IF (L .EQ. 1) GO TO 270
         F = 0.0D0
C
        DO 240 J = 1, L 
            G = 0.0D0
C     ********** FORM ELEMENT OF A*U ********** 
            DO 180 K = 1, J
  180       G = G + A(J,K) * A(I,K)
C
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
  200       G = G + A(K,J) * A(I,K)
C     ********** FORM ELEMENT OF P **********
  220       E(J) = G / H
            F = F + E(J) * A(I,J)
  240    CONTINUE
C
         H = F / (H + H)
C     ********** FORM REDUCED A **********
         DO 260 J = 1, L 
            F = A(I,J)
            G = E(J) - H * F
            E(J) = G
C
            DO 260 K = 1, J
               A(J,K) = A(J,K) - F * E(K) - G * A(I,K)
  260    CONTINUE 
C
  270    DO 280 K = 1, L
  280    A(I,K) = SCALE * A(I,K)
C
  290    H = D(I)
         D(I) = A(I,I)
         A(I,I) = H
  300 CONTINUE
C 
      RETURN
C     ********** LAST CARD OF TRED1 **********
      END
C
C
      SUBROUTINE TRED2(NM,N,A,D,E,Z)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(NM,N),D(N),E(N),Z(NM,N) 
C     **************************************************************** 
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A 
C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT :
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE 
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C     ON OUTPUT :
C 
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX, 
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
C 
C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX 
C          PRODUCED IN THE REDUCTION,
C
C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
C 
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     *************************************************************
      DO 100 I = 1, N 
         DO 100 J = 1, I
            Z(I,J) = A(I,J)
  100 CONTINUE 
C
      IF (N .EQ. 1) GO TO 320
C     ********** FOR I=N STEP -1 UNTIL 2 DO -- ********** 
      DO 300 II = 2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 2) GO TO 130
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
         DO 120 K = 1, L
  120    SCALE = SCALE + DABS(Z(I,K))
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = Z(I,L)
         GO TO 290
C
  140    DO 150 K = 1, L 
            Z(I,K) = Z(I,K) / SCALE 
            H = H + Z(I,K) * Z(I,K)
  150    CONTINUE
C 
         F = Z(I,L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         Z(I,L) = F - G
         F = 0.0D0
C 
         DO 240 J = 1, L 
            Z(J,I) = Z(I,J) / H
            G = 0.0D0
C     ********** FORM ELEMENT OF A*U ********** 
            DO 180 K = 1, J
  180       G = G + Z(J,K) * Z(I,K)
C
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
  200       G = G + Z(K,J) * Z(I,K)
C     ********** FORM ELEMENT OF P **********
  220       E(J) = G / H 
            F = F + E(J) * Z(I,J)
  240    CONTINUE
C  
         HH = F / (H + H) 
C     ********** FORM REDUCED A **********
         DO 260 J = 1, L 
            F = Z(I,J) 
            G = E(J) - HH * F
            E(J) = G
C
            DO 260 K = 1, J
               Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260    CONTINUE
C 
  290    D(I) = H 
  300 CONTINUE
C 
  320 D(1) = 0.0D0
      E(1) = 0.0D0
C     ********** ACCUMULATION OF TRANSFORMATION MATRICES **********
      DO 500 I = 1, N
         L = I - 1
         IF (D(I) .EQ. 0.0D0) GO TO 380
C
         DO 360 J = 1, L
            G = 0.0D0
C 
            DO 340 K = 1, L
  340       G = G + Z(I,K) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * Z(K,I)
  360    CONTINUE
C
  380    D(I) = Z(I,I)
         Z(I,I) = 1.0D0
         IF (L .LT. 1) GO TO 500
C
         DO 400 J = 1, L
            Z(I,J) = 0.0D0
            Z(J,I) = 0.0D0
  400    CONTINUE
C                
  500 CONTINUE
C
      RETURN 
C     ********** LAST CARD OF TRED2 **********
      END
C
C
      SUBROUTINE TQLRAT(N,D,E2,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N),E2(N) 
      REAL*8 MACHEP
C     *****************************************************************
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
C     ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
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,
C 
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX'
C
C        E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE
C          INPUT MATRIX IN ITS LAST N-1 POSITIONS.  E2(1) IS ARBITRARY.
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        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     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C     **************************************************************
C 
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C 
C
      MACHEP = 2.D0**(-26)
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.0D0
      B = 0.0D0
      C = 0.0D0
      E2(N) = 0.0D0
C
      DO 290 L = 1, N
         J = 0
         H = MACHEP * (DABS(D(L)) + DSQRT(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
         WRITE(6,*)' **** FATAL ERROR IN TQLRAT **** ' 
         WRITE(6,*)' **** FALLEN THROUGH BOTTOM OF LOOP 110 *** '
         STOP
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 = DSQRT(E2(L))
         G = D(L)
         P = (D(L1) - G) / (2.0D0 * S)
         R = DSQRT(P*P+1.0D0)
         D(L) = S / (P + DSIGN(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.0D0) G = B
         H = G
         S = 0.0D0
         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.0D0) 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.0D0) GO TO 210
         IF (DABS(E2(L)) .LE.DABS(C/H)) GO TO 210
         E2(L) = H * E2(L)
         IF (E2(L) .NE. 0.0D0) 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 
C     ********** LAST CARD OF TQLRAT **********
      END 
C 
C  
      SUBROUTINE TQL2(NM,N,D,E,Z,IERR) 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N),E(N),Z(NM,N)
      REAL*8 MACHEP
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 TWO-DIMENSIONAL 
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT,
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        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.
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 SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, 
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ***************************************************************
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.

      MACHEP = 2.D0**(-26)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C 
      DO 100 I = 2, N 
  100 E(I-1) = E(I) 
C
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
C
      DO 240 L = 1, N
         J = 0 
         H = MACHEP * (DABS(D(L)) + DABS(E(L)))
         IF (B .LT. H) B = H
C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
         DO 110 M = L, N
            IF (DABS(E(M)) .LE. B) GO TO 120
C     ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT 
C                THROUGH THE BOTTOM OF THE LOOP **********
  110    CONTINUE 
  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
         G = D(L)
         P = (D(L1) - G) / (2.0D0 * E(L)) 
         R = DSQRT(P*P+1.0D0)
         D(L) = E(L) / (P + DSIGN(R,P))
         H = G - D(L)
         DO 140 I = L1, N
  140    D(I) = D(I) - H 
C
         F = F + H 
C     ********** QL TRANSFORMATION ********** 
         P = D(M)
         C = 1.0D0
         S = 0.0D0
         MML = M - L 
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
         DO 200 II = 1, MML 
            I = M - II
            G = C * E(I)
            H = C * P 
            IF (DABS(P) .LT. DABS(E(I))) GO TO 150 
            C = E(I) / P
            R = DSQRT(C*C+1.0D0)
            E(I+1) = S * P * R
            S = C / R 
            C = 1.0D0 / R
            GO TO 160 
  150       C = P / E(I)
            R = DSQRT(C*C+1.0D0)
            E(I+1) = S * E(I) * R
            S = 1.0D0 / 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
  200    CONTINUE
         E(L) = S * P
         D(L) = C * P 
         IF (DABS(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
         DO 280 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P 
  280    CONTINUE
  300 CONTINUE
      GO TO 1001
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ********** 
 1000 IERR = L 
 1001 RETURN
C     ********** LAST CARD OF TQL2 **********
      END
