SUBROUTINE SLLDEF(ELAM,EPS,A,B,K,TOL,IFAIL)
      DOUBLE PRECISION ELAM,EPS,A,B,TOL
      INTEGER K,IFAIL
C
C This subroutine computes numerically the eigenvalues of a 
C left-definite right-indefinite regular Sturm-Liouville problem. 
C
C It is written in FORTRAN 77 and requires that the user already
C have installed LAPACK version 3. If you do not have LAPACK, go
C to
C       www.netlib.org/lapack
C
C and follow the instructions for installing the software on your
C system. SLLDEF must be linked with the LAPACK library at compile
C time.
C 
C 
C
C The problems solved by SLLDEF consist of a differential equation
C
C   -(p(x)y')' + q(x)y = lambda w(x) y,   x in (a,b),
C
C to be satisfied together with coupled boundary conditions of
C the form
C
C    c11 y(b) + c12 y(a) = d11 py'(b) - d12 py'(a)
C    c21 y(b) + c22 y(a) = d21 py'(b) - d22 py'(a)
C    
C Problems of this form have eigenvalues which are unbounded both
C above and below and which can be indexed according to the scheme
C 
C  \lambda_{-k} \leq \lambda{-k+1} \leq ... \leq \lambda_{-1}
C      \leq \lambda_{-0} < 0 < \lambda_{+0} \leq \lambda_{1} \leq ...
C         \leq \lambda_{k} \leq ...
C
C  
C Parameters:
C -----------
C 
C ELAM: (input/output, double precision real): On input with K
C     nonzero, ELAM can be assigned any value. On input with K=0,
C     ELAM should be assigned a negative or a positive value depending
C     on whether \lambda_{-0} or \lambda_{+0} is to be calculated.
C
C     On exit with IFAIL=0: the computed eigenvalue approximation.
C
C EPS:  (input/output, double precision real): On input, EPS should
C     be assigned a positive value. In general, the routine will work
C     faster when EPS is a good estimate of the distance between the
C     input value of EPS and the eigenvalue sought. On output, EPS is
C     overwritten.
C
C A,B:  (input, double precision reals): On input, the endpoints of
C     the interval on which the problem is posed. The interval must be
C     finite with A < B.
C
C K:    (input, integer): The index of the eigenvalue sought. This may
C     be any integer, but note that if K=0 then the input value of ELAM
C     must be nonzero.
C 
C TOL:  (input, double precision real): The tolerance for the numerical
C     integration. TOL must be strictly positive. In general, reducing 
C     the value of TOL will improve the accuracy of the computed 
C     eigenvalue.
C 
C IFAIL: (input/output, integer): On entry, IFAIL must have one of the
C     values -1, 0 or 1. On exit, IFAIL will be 0 unless an error is
C     detected. If IFAIL was 0 or 1 on entry then an error message will
C     be printed in the event of an error. If IFAIL was 0 on entry then
C     the execution of the program will be terminated.
C
C User Supplied Routines
C ----------------------
C
C The user must supply the coefficients p(x), q(x) and w(x) appearing
C in the differential equation by writing DOUBLE PRECISION FUNCTIONS 
C having the names P, Q and W as follows:
C
C     DOUBLE PRECISION FUNCTION P(X)
C     DOUBLE PRECISION X
C 
C     P = (assign value as a function of X)
C 
C     RETURN
C     END
C
C     DOUBLE PRECISION FUNCTION Q(X)
C     DOUBLE PRECISION X
C 
C     Q = (assign value as a function of X)
C 
C     RETURN
C     END
C
C     DOUBLE PRECISION FUNCTION W(X)
C     DOUBLE PRECISION X
C 
C     W = (assign value as a function of X)
C 
C     RETURN
C     END
C
C The functions P, Q and W must be such as to make the problem left
C definite and must be such that W is negative on a set of positive
C measure and positive on a set of positive measure. In practice the
C performance of the code will be much slower if the coefficients 
C P, Q and W are not smooth functions.
C
C The user must also write a subroutine called BCS to supply the 
C coefficients c11, c12, c21, c22, d11, d12, d21 and d22 appearing
C in the boundary conditions. Its specification must be as follows.
C
C     SUBROUTINE BCS(C,D,N)
C     INTEGER N
C     DOUBLE PRECISION C(N,N),D(N,N)
C
C     C(1,1) = c11
C     C(1,2) = c12
C 
C     (etc).
C
C     RETURN
C     END
C
C The integer N is an input parameter. Its value must not be changed
C by BCS. 
C
C
C The matrices C and D returned by this routine must satisfy the 
C condition
C
C   C(1,1)*D(2,1)+C(1,2)*D(2,2)-C(2,1)*D(1,1)-C(2,2)*D(1,2) = 0.
C
C This is necessary for self-adjointness.
C 
C ---------------------------------------------------------------------
C
      EXTERNAL USUB1,USUB2
      INTEGER N
      PARAMETER (N=2)
      INTEGER KNOBS,IWORK,ILWK
      PARAMETER (IWORK=N*(16+73*N)+65,ILWK=N,KNOBS=60)
      INTEGER LWK(ILWK)
      DOUBLE PRECISION WORK(IWORK)
      COMPLEX*16 CWORK(N*(N+65))

      CALL SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,WORK,
     &           IWORK,LWK,ILWK,CWORK,IFAIL)

      RETURN
      END

      SUBROUTINE USUB1(X,ELAM,S11,S12,S21,S22,N)
      DOUBLE PRECISION X,ELAM,S11(N,N),S12(N,N),S21(N,N),S22(N,N)
      DOUBLE PRECISION P,Q,W
      EXTERNAL P,Q,W
C Remark: The differential equation is -(Py')'+Qy = ELAM.W.y.
C The DOUBLE PRECISION functions providing the coefficient values
C MUST be called P, Q and W.

      DO 10 I = 1,N
         DO 20 J = 1,N
            S11(J,I) = 0.D0
            S12(J,I) = 0.D0
            S21(J,I) = 0.D0
            S22(J,I) = 0.D0
 20      CONTINUE
 10   CONTINUE

      S11(1,1) = ELAM*W(X)-Q(X)
      S22(1,1) = 1.D0/P(X)
      RETURN
      END

      SUBROUTINE USUB2(IEND,ISING,X,U,V,NDIM,N,ELAM)
      INTEGER IEND,NDIM,N
      LOGICAL ISING
      DOUBLE PRECISION U(NDIM,N),V(NDIM,N)
      DOUBLE PRECISION C(N,N),D(N,N)
      DOUBLE PRECISION TEST
      ISING = .FALSE.
      IF (IEND.EQ.0) THEN
C Provide boundary condition at left hand end:
         U(1,1) = 1.D0
         U(2,1) = 1.D0
         U(1,2) = 0.D0
         U(2,2) = 0.D0
         V(1,1) = 0.D0
         V(2,1) = 0.D0
         V(1,2) = 1.D0
         V(2,2) = -1.D0
      ELSE
         CALL BCS(C,D,N)
C Remark: the routine BCS returns the 2x2 matrices C and D which
C specify the boundary conditions in the form
C
C    c11 y(b) + c12 y(a) = d11 py'(b) - d12 py'(a)
C    c21 y(b) + c22 y(a) = d21 py'(b) - d22 py'(a)
C
C The routine MUST be called BCS and, for the time being, only real
C boundary conditions are supported.
C
C
         TEST = C(1,1)*D(2,1)+C(1,2)*D(2,2)
     &          -C(2,1)*D(1,1)-C(2,2)*D(1,2)
         IF (ABS(TEST).GT.1.0D-6) THEN
             WRITE (6,FMT=*)
     &       'Your boundary conditions are not selfadjoint'
             STOP
         END IF
         U(1,1) = D(1,1)
         U(2,1) = D(1,2)
         U(1,2) = D(2,1)
         U(2,2) = D(2,2)
         V(1,1) = C(1,1)
         V(2,1) = C(1,2)
         V(1,2) = C(2,1)
         V(2,2) = C(2,2)
      END IF

      RETURN
      END

C ---------------------------------------------------------------------

      SUBROUTINE SL11F(ELAM,EPS,A,B,K,N,TOL,USUB1,USUB2,KNOBS,WORK,
     &                 IWORK,LWK,ILWK,CWORK,IFAIL)
C
C Bug fixed on 8th September 1997
C
C Specification: SL11F solves Hamiltonian eigenvalue problems which are
C reducible to the form of a first order matrix differential equation
C Theta' = i Theta Omega(x,lambda,Theta),
C where Theta is a unitary matrix prescribed by boundary conditions at
C x=a and x=b.
C
C Variables:
C elam: double precision. On entry elam should be set to an approximation
C       to the eigenvalue to be found. The value of elam is not important
C       in determining the accuracy of the final result but it may affect
C       run times. On exit, elam specifies the approximation to the
C       eigenvalue computed by SL11F.
C eps:  double precision. On entry eps should be set to a positive real
C       number of an order of magnitude such that the true eigenvalue
C       lies within [elam-C*eps,elam+C*eps] where by elam we mean the
C       initial approximation provided by the user, and C is a constant
C       of modest size. Again the value of eps is not critical. In cases
C       where absolutely no prior information is available it would be
C       acceptable in almost all cases to set elam = 0 and eps = 10 on
C       entry.
C a,b:  double precision. On entry a and b specify the finite regular endpoints
C       of the interval on which the problem is posed. Unchanged on exit.
C k:    integer. The index of the eigenvalue sought. For many problems k=0
C       specifies the lowest eigenvalue (e.g. for problems derived from
C       formally self-adjoint even-order differential operators). However
C       for other problems there exist eigenvalues for any integer value
C       of the index k, positive or negative. Unchanged on exit.
C n:    integer. The dimension of the system, i.e. the size of the matrix
C       Theta. Unchanged on exit.
C tol:  double precision. On entry tol should have a positive value. It is
C       intended that the smaller tol is, the smaller will be the error in
C       the final approximation elam returned by the code. The value of
C       tol is unchanged on exit.
C knobs: integer. Used to control the maximum number of miss-distance
C       evaluations carried out by the code. The number of miss-distance
C       evaluations will be max(60,knobs). knobs is unchanged on exit.
C work: double precision array of dimension iwork. Used as workspace.
C iwork: integer. On entry iwork must be at least n*(16+n*(33+2*nsampl))+25+
C       2*nsampl. Unchanged on exit.
C lwk:  integer array of dimension ilwk. Used as workspace.
C ilwk: integer. On entry ilwk must be at least n. Unchanged on exit.
C cwork: complex*16 array of dimension at least N*(N+65). Used as workspace.
C ifail: integer. On entry, ifail must have one of the values -1,0,1.
C       On successful exit ifail will have the value 0. Any other value
C       indicates a failure.
C
C User Supplied Subroutines USUB1 and USUB2.
C
C USUB1:
C
C     SUBROUTINE usub1(x,elam,s11,s12,s21,s22,n)
C     INTEGER n
C     DOUBLE PRECISION x,elam,s11(n,n),s12(n,n),s21(n,n),s22(n,n)
C
C This routine must specify in s11, s12, s21 and s22 the values of the
C four coefficient matrix functions of x and lambda (elam) occurring
C in the Hamiltonian system. It must not change the values of x, elam,
C or n.
C
C USUB2:
C
C     SUBROUTINE usub2(iend,ising,x,u,v,ndim,n,elam)
C     INTEGER iend,ndim,n
C     LOGICAL ising
C     DOUBLE PRECISION x,elam,u(ndim,n),v(ndim,n)
C
C If iend=0 on entry then USUB2 must specify the boundary condition at
C x=a, while if iend=1 then USUB2 must specify the boundary condition
C at x=b. The parameter ising must be assigned the value .false. before
C returning. The parameters x, elam, ndim and n must not be changed by
C USUB2. The specification of boundary conditions is achieved as follows:
C a boundary condition
C                      transpose(A1).U + transpose(A2).V = 0
C is specified by setting
C                         U = A2, V = -A1
C and returning.
C
C END of brief documentation.
C
C     .. Parameters ..
      INTEGER NSAMPL
      PARAMETER (NSAMPL=20)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION C,EL,EU,TOLC05
      INTEGER IFLAG,IFO,IREST,ISAMPL,ISIGIL,ISIGIR,ISIGRL,ISIGRR,KNTMAX,
     &        NREC
      CHARACTER*6 SRNAME
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION A,B,ELAM,EPS,TOL
      INTEGER IFAIL,ILWK,IWORK,K,KNOBS,N
C     ..
C     .. Array Arguments ..
C IWORK must be at least n*(16+(33+2*nsampl)*n)+25+2*nsampl. With nsampl=10
C this gives n*(16+53*n)+45, with nsampl=20 we get n*(16+73*n)+65
      DOUBLE PRECISION WORK(IWORK)
      INTEGER LWK(ILWK)
      COMPLEX*16 CWORK(N*(N+65))
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL USUB1,USUB2
C     ..
C     .. External Subroutines ..
      EXTERNAL CHOOSC,SOLVSY
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC MAX,SQRT
C     ..
C     .. External Functions ..
      DOUBLE PRECISION XM2AJF
      INTEGER PM1ABF
      EXTERNAL XM2AJF,PM1ABF
C     ..
C     .. Local Arrays ..
      CHARACTER*80 REC(2)
C     ..
      SRNAME = 'SL11F'
      IFO = IFAIL

C First do the input checks.
      IF (A.GE.B .OR. N.LT.2 .OR. TOL.LE.0.0D0 .OR.
     &    IWORK.LT.N* (16+ (33+2*NSAMPL)*N)+25+2*NSAMPL .OR.
     &    ILWK.LT.N .OR. EPS.LE.0.0D0) THEN
C We have an input error
          IFLAG = 1
          IF (A.GE.B) THEN
              WRITE (REC,FMT=9000) A,B
              NREC = 2
              GO TO 10
          ELSE IF (N.LT.2) THEN
              WRITE (REC,FMT=9010) N
              NREC = 2
              GO TO 10
          ELSE IF (TOL.LE.0.0D0) THEN
              WRITE (REC,FMT=9020) TOL
              NREC = 2
              GO TO 10
          ELSE IF (IWORK.LT.N* (16+ (33+2*NSAMPL)*N)+25+2*NSAMPL) THEN
              WRITE (REC,FMT=9030) N* (16+ (33+2*NSAMPL)*N) + 25 +
     &          2*NSAMPL,IWORK
              NREC = 2
              GO TO 10
          ELSE IF (ILWK.LT.N) THEN
              WRITE (REC,FMT=9040) N,ILWK
              NREC = 2
              GO TO 10
          ELSE IF (EPS.LE.0.0D0) THEN
              WRITE (REC,FMT=9050) EPS
              NREC = 2
              GO TO 10
          ELSE

          END IF
      END IF

      ISIGRL = 1
      ISIGIL = N*N + 1
      ISIGRR = 2*N*N + 1
      ISIGIR = 3*N*N + 1
      IREST = 4*N*N + 1
      ISAMPL = 40*N*N + 4*N + 1
      KNTMAX = MAX(60,KNOBS)
      TOLC05 = 1.0D-2
      IFAIL = 0
      C = 0.50D0* (A+B)
      CALL SOLVSY(A,B,C,K,ELAM,EPS,EL,EU,USUB1,USUB2,TOL,TOLC05,
     &            WORK(IREST),WORK(ISIGRL),WORK(ISIGIL),N,WORK(ISIGRR),
     &            WORK(ISIGIR),N,LWK,N,KNTMAX,CWORK,IFAIL)
      IF (IFAIL.EQ.0) THEN
C Choose a better matching point and repeat
          CALL CHOOSC(A,B,C,EL,EU,TOL,USUB1,USUB2,WORK,LWK,N,K,NSAMPL,
     &                WORK(N* (16+33*N)+25+2*NSAMPL+1),CWORK,
     &                IFAIL)

          IF (IFAIL.EQ.0) THEN

              TOLC05 = MAX(1.0D-1*TOL,SQRT(XM2AJF(0.0D0)))
              ELAM = 0.50D0* (EU+EL)
              EPS = 0.50D0* (EU-EL)

              CALL SOLVSY(A,B,C,K,ELAM,EPS,EL,EU,USUB1,USUB2,TOL,TOLC05,
     &                    WORK(IREST),WORK(ISIGRL),WORK(ISIGIL),N,
     &                    WORK(ISIGRR),WORK(ISIGIR),N,LWK,N,KNTMAX,
     &                    CWORK,IFAIL)
              IF (IFAIL.EQ.0) THEN

                  RETURN
              END IF
          END IF
      END IF

      IFLAG = IFAIL
      IF (IFAIL.EQ.7) THEN
          WRITE (REC,FMT=9060)
          NREC = 2
          IFLAG = 2
      ELSE IF (IFAIL.EQ.12 .OR. IFAIL.EQ.14) THEN
          WRITE (REC,FMT=9070)
          NREC = 2
          IFLAG = 3
      ELSE IF (IFAIL.EQ.13 .OR. IFAIL.EQ.15) THEN
          WRITE (REC,FMT=9080)
          NREC = 2
          IFLAG = 4
      ELSE IF (IFAIL.EQ.10 .OR. IFAIL.EQ.11 .OR. IFAIL.EQ.16 .OR.
     &         IFAIL.EQ.17) THEN
          WRITE (REC,FMT=9090)
          NREC = 2
          IFLAG = 5
      ELSE IF (IFAIL.EQ.9 .OR. IFAIL.EQ.18) THEN
          WRITE (REC,FMT=9100)
          NREC = 2
          IFLAG = 6
      ELSE IF (IFAIL.EQ.20) THEN
          WRITE (REC,FMT=9110)
          NREC = 2
          IFLAG = 7
      ELSE IF (IFAIL.EQ.21) THEN
          WRITE (REC,FMT=9120)
          NREC = 2
          IFLAG = 8
      ELSE IF (IFAIL.EQ.19) THEN
          WRITE (REC,FMT=9130)
          NREC = 2
          IFLAG = 9
      END IF

   10 IFAIL = PM1ABF(IFO,IFLAG,SRNAME,NREC,REC)

 9000 FORMAT (' ** A must be less than B, you have',/,' ** A = ',g18.8,
     &       ' B = ',g18.8)
 9010 FORMAT (' ** N must be at least 2, you have',/,' ** N = ',i8)
 9020 FORMAT (' ** TOL must be strictly positive, you have',/,
     &       ' ** TOL = ',g18.8)
 9030 FORMAT (' ** IWORK must be at least ',i8,' you have',/,
     &       ' ** IWORK = ',i8)
 9040 FORMAT (' ** ILWK must be at least',i8,' you have',/,
     &       ' ** ILWK = ',i8)
 9050 FORMAT (' ** EPS must be strictly positive, you have',/,
     &       ' ** EPS = ',g18.8)
 9060 FORMAT (' ** Too many attempts have been made to find',/,
     &       ' ** required eigenvalue; try increasing EPS')
 9070 FORMAT (' ** Failure in FM2BJF when trying to convert boundary',/,
     &       ' ** conditions from user variables to H-matrix variables')
 9080 FORMAT (' ** The boundary conditions supplied do not',/,
     &       ' ** satisfy the rank condition')
 9090 FORMAT (' ** An H-matrix has occurred whose exponential',/,
     &       ' ** cannot be accurately computed')
 9100 FORMAT (' ** A failure has occurred on a call to FM2AJF;',/,
     &       ' ** cannot find eigenvalues of central miss-matrix')
 9110 FORMAT (' ** An H-matrix has arisen which cannot be',/,
     &       ' ** diagonalised; try reducing TOL')
 9120 FORMAT (' ** Too many unsuccessful attempts have been made',/,
     &       ' ** to find an initial stepsize')
 9130 FORMAT (' ** Too many unsuccessful attempts have been made',/,
     &       ' ** to find a suitable next stepsize')
      END


C ----------------------------------------------------------------------


      SUBROUTINE SOLVSY(A,B,C,K,ELAM,EPS,EL,EU,COFMAT,SYSBC,TOL,TOLC05,
     &                  WORK,SIGRL,SIGIL,ISIGL,SIGRR,SIGIR,ISIGR,LWK,N,
     &                  KNTMAX,CWORK,IFAIL)
C This subroutine is a modified version of the usual subroutine by the same
C name provided with the package SL11F. It has been modified to cope with
C left-definite 2nd order Sturm-Liouville problems with possibly-coupled
C boundary conditions and assumes that the miss distance function
C D(lambda) is monotone decreasing when lambda < 0, monotone increasing
C when lambda > 0, and zero when lambda lies in the range (lambda_{-0},
C lambda_{+0}), where lambda_{-0}<0<lambda_{+0}.
C
C Created by Marco Marletta on the 8th of April 2001.
C

C First find a bracket for the eigenvalue:
C     .. Scalar Arguments ..
      DOUBLE PRECISION A,B,C,EL,ELAM,EPS,EU,TOL,TOLC05
      INTEGER IFAIL,ISIGL,ISIGR,K,KNTMAX,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION SIGIL(ISIGL,N),SIGIR(ISIGR,N),SIGRL(ISIGL,N),
     &                 SIGRR(ISIGR,N),WORK(N* (13+25*N)+25)
      INTEGER LWK(N)
      COMPLEX*16 CWORK(N*(N+65))
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL COFMAT,SYSBC
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION DL,DU,E1,E2,SHIFT,TOLOC,E3,F3
      INTEGER IFLAG,IND,KNT
      LOGICAL DUAS,TREVAL,POSEIG,NEGEIG
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION C1(17)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION DLH13F
      EXTERNAL DLH13F
C     ..
C     .. External Subroutines ..
C      EXTERNAL C05AZF
      EXTERNAL C05TST
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,DBLE,MAX,MIN
C     ..
      IF (K.GT.0.OR.(K.EQ.0.AND.ELAM.GT.0.D0)) THEN
          SHIFT = DBLE(K) + 5.0D-1
          POSEIG = .TRUE.
      ELSE
          SHIFT = DBLE(ABS(K)) + 5.0D-1
          POSEIG = .FALSE.
      END IF
      NEGEIG = .NOT.POSEIG
      EL = ELAM - ABS(EPS)
      EU = ELAM + ABS(EPS)
      IF (POSEIG) THEN
          EL = MAX(0.D0,EL)
          EU = EL + 2.D0*ABS(EPS)
      ELSE
          EU = MIN(0.D0,EU)
          EL = EU - 2.D0*ABS(EPS)
      END IF

      DUAS = .false.
      KNT = 0
      IFAIL = 0
   10 CONTINUE

      WRITE (6,FMT=*) 'el=',el

      DL = DLH13F(A,B,C,EL,N,TOL,SIGRL,ISIGL,SIGRR,ISIGR,WORK(1),
     &     WORK(N*N+1),WORK(2*N*N+1),WORK(3*N*N+1),WORK(4*N*N+1),LWK,
     &     CWORK,COFMAT,SYSBC,SHIFT,IFAIL)

      IF (IFAIL.NE.0) THEN
          RETURN
C      WRITE (6,FMT=*) 'elam,dmiss:',el,dl
      ELSE IF ((DL.GT.0.0D0.AND.POSEIG).OR.(DL.LT.0.D0.AND.NEGEIG))
     &                                                         THEN
          WRITE (6,FMT=*) 'elam,dmiss:',el,dl
          EU = EL
          DU = DL
          DUAS = .true.
          EPS = 2.0D0*ABS(EPS)
          EL = EL - EPS
          GO TO 10
      END IF

      IF (.NOT.DUAS) THEN
          DU = DLH13F(A,B,C,EU,N,TOL,SIGRL,ISIGL,SIGRR,ISIGR,WORK(1),
     &         WORK(N*N+1),WORK(2*N*N+1),WORK(3*N*N+1),WORK(4*N*N+1),
     &         LWK,CWORK,COFMAT,SYSBC,SHIFT,IFAIL)

          IF (IFAIL.NE.0) THEN
              RETURN
          ELSE
   20         CONTINUE
      WRITE (6,FMT=*) 'elam,dmiss:',eu,du
              IF (DL*DU.GT.0.0D0) THEN
C We have not yet bracketed an eigenvalue
                  KNT = KNT + 1
                  IF (KNT.LE.KNTMAX) THEN

                      EPS = 2.0D0*ABS(EPS)
                      IF ((DU.LT.0.0D0.AND.POSEIG).OR.
     &                    (DU.GT.0.D0.AND.NEGEIG)) THEN
                          EL = EU
                          DL = DU
                          EU = EU + ABS(EPS)

                          DU = DLH13F(A,B,C,EU,N,TOL,SIGRL,ISIGL,SIGRR,
     &                         ISIGR,WORK(1),WORK(N*N+1),WORK(2*N*N+1),
     &                         WORK(3*N*N+1),WORK(4*N*N+1),LWK,CWORK,
     &                         COFMAT,SYSBC,SHIFT,IFAIL)

                          IF (IFAIL.NE.0) RETURN
                      ELSE
              WRITE (6,FMT=*) 'elam,dmiss:',eu,du

                          EU = EL
                          DU = DL
                          EL = EL - ABS(EPS)

                          DL = DLH13F(A,B,C,EL,N,TOL,SIGRL,ISIGL,SIGRR,
     &                         ISIGR,WORK(1),WORK(N*N+1),WORK(2*N*N+1),
     &                         WORK(3*N*N+1),WORK(4*N*N+1),LWK,CWORK,
     &                         COFMAT,SYSBC,SHIFT,IFAIL)

                          IF (IFAIL.NE.0) THEN
                              RETURN
                          ELSE
              WRITE (6,FMT=*) 'elam,dmiss:',el,dl
                          END IF
                      END IF

                      GO TO 20
                  END IF
              ELSE
                  GO TO 30
              END IF
              IFAIL = 7
              RETURN
          END IF
      END IF

C If we are here then it means that we have bracketed an eigenvalue and
C we can now think about locating it accurately with c05azf.
   30 KNT = 0
      IND = -1
      TOLOC = TOLC05
      C1(1) = DU
      IFLAG = 1
      E1 = EL
      E2 = EU
      TREVAL = .FALSE.
   40 CONTINUE

C      CALL C05AZF(E1,E2,DL,TOLOC,1,C1,IND,IFLAG)
      CALL C05TST(E1,E2,E3,F3,TREVAL,DL,TOLOC,1,C1,IND,IFLAG)

      IF (IND.EQ.0) THEN
          GO TO 50
      ELSE

          KNT = KNT + 1
          IF (KNT.LE.KNTMAX) THEN

              DL = DLH13F(A,B,C,E1,N,TOL,SIGRL,ISIGL,SIGRR,ISIGR,
     &             WORK(1),WORK(N*N+1),WORK(2*N*N+1),WORK(3*N*N+1),
     &             WORK(4*N*N+1),LWK,CWORK,COFMAT,SYSBC,SHIFT,
     &             IFAIL)
                      WRITE (6,FMT=*) 'elam,dmiss:',E1,DL

              IF (IFAIL.NE.0) THEN
                  RETURN
              ELSE

                  IF ((DL.GT.0.0D0.AND.POSEIG).OR.
     &                (DL.LT.0.D0.AND.NEGEIG)) EU = MIN(EU,E1)
                  IF ((DL.LT.0.0D0.AND.POSEIG).OR.
     &                (DL.GT.0.D0.AND.NEGEIG)) EL = MAX(EL,E1)
                  IF (ABS(EU-EL).LT.TOLOC*MAX(1.0D0,ABS(EL))) THEN
                      GO TO 50
                  ELSE
                      WRITE (6,FMT=*) 'elam,dmiss:',E1,DL

                      GO TO 40
                  END IF
              END IF
          END IF
      END IF
      IFAIL = 7
      RETURN

   50 ELAM = E1

      END


C ------------------------------------------------------------------------

        SUBROUTINE c05tst(x,y,z,fz,treval,fx,tolx,ir,c,ind,ifail)
C First written by Chris Graves (ESF Student, Leicester University) 
C on 13/10/92. Modified to cope with `nasty' functions by switching
C to bisection, by Marco Marletta on 7th June 1993. Specification
C similar to NAG routine C05AZF.

        DOUBLE PRECISION x,y,z,fx,fz,tolx,c(17)
        INTEGER ir,ind,ifail
        LOGICAL treval
        DOUBLE PRECISION x0,x1,valx0,valx1,new,new2,x1old
        SAVE x0,x1,valx0,valx1,new
        GOTO (100,110,120,130,140,150) ind+2
C       (IND=-1)
100     x0=x
        x1=y
        z = x0
        valx0=fx
        valx1=c(1)
        fz = valx0
        GOTO 150
c       (IND=0)
110     RETURN
c       (IND=1)
120     x0=x
        x1=y
        ind=2
        RETURN
c       (IND=2)
130     valx0=fx
        z = x
        fz = fx
        x=x1
        ind=3
        RETURN
c       (IND=3)
140     valx1=fx
c       (IND=4)
150     new=x1-(valx1*(x1-x0)/(valx1-valx0))
        IF (treval) THEN
            new2=x1-(valx1*(x1-z)/(valx1-fz))
        ELSE
            new2 = new
        END IF
        IF (sign(1.D0,valx0)*valx1.LT.0.D0) THEN
            z = x0
            fz = valx0
        END IF
        treval = .true.
        x0 = x1
        valx0=valx1
        IF (abs(new-new2).LT.0.25D0*abs(x1-z)) THEN
           x1=new
        ELSE
            x1 = 0.5D0*(x1+z)
        END IF
c       ..CHECK TOLERANCE..
        IF (ABS(x1-z).le.tolx*max(1.D0,x0,x1)) THEN
          ind=0
          RETURN
        END IF
        x=x1
        ind=3
        RETURN
        END




C ------------------------------------------------------------------------


      SUBROUTINE CHOOSC(A,B,C,EL,EU,TOL,COFMAT,SYSBC,WORK,LWK,N,K,
     &                  NSAMPL,WSAMPL,CWORK,IFAIL)
C Subroutine to determine the best matching point for a shooting algorithm
C to solve a Hamiltonian System eigenvalue problem. Uses a lot of storage,
C unfortunately.
C     .. Parameters ..
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION A,B,C,EL,EU,TOL
      INTEGER IFAIL,K,N,NSAMPL
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION WORK(N* (16+33*N)+25+2*NSAMPL),
     &                 WSAMPL(1:N,1:N,1:2*NSAMPL)
      INTEGER LWK(N)
      COMPLEX*16 CWORK(N*(N+65))
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL COFMAT,SYSBC
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ADJUST,ALFMAX,ALFMIN,ALPHA,ARGC,ARGL,ARGR,BETA,
     &                 DIFF,DLA13F,ELAM,FAC,HSAMPL,OMEGA,SHIFT,TWOPI,X1,
     &                 X2
      INTEGER I,IALFI,IALFR,IASMP,IBETA,IDLSMP,IHL,II,INDEX,ISAMPI,
     &        ISAMPR,IUL,IULO,IUR,IV,IVL,IVLO,IVR,J,JMATCH,JSAMPL,KNT,N2
      LOGICAL ISING
C     ..
C     .. External Functions ..
      DOUBLE PRECISION XM1AAF
      EXTERNAL XM1AAF
C     ..
C     .. External Subroutines ..
      EXTERNAL FM2AJF,FM2BJF,SMCOPY,SMLOAD,DGEMM,HINT,MATEXP,MVN,UDUT
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN2,DBLE,SIGN
C     ..
      IFAIL = 0
      SHIFT = DBLE(K) + 5.0D-1
      TWOPI = 2.0D0*XM1AAF(0.0D0)
      N2 = N*N
C Set up pointers for the workspace:

      IUL = N* (13+25*N) + 25 + 1
      IVL = IUL + N2
      IUR = IVL + N2
      IVR = IUR + N2
      IULO = IVR + N2
      IVLO = IULO + N2
      IALFR = IVLO + N2
      IALFI = IALFR + N
      IBETA = IALFI + N
      IV = IBETA + N
      IHL = IV + N2
      IASMP = IHL + N2
      IDLSMP = IASMP + NSAMPL

C Need workspace of length n*(16+33*n)+25+2*nsampl

      ELAM = EL
      DIFF = 100.0D0
   10 CONTINUE
C Set up initial condition at x=a:
      CALL SYSBC(0,ISING,A,WORK(IUL),WORK(IVL),N,N,ELAM)
C Now we convert the boundary condition information into information
C about the H-matrix.

      CALL SMCOPY('G',N,N,WORK(IUL),N,WORK(IULO),N)
      CALL SMCOPY('G',N,N,WORK(IVL),N,WORK(IVLO),N)
   20 CONTINUE
      CALL FM2BJF(N,WORK(IULO),N,WORK(IVLO),N,0.0D0,WORK(IALFR),
     &            WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,
     &            WORK(1),IFAIL)
      IF (IFAIL.NE.0) THEN
          KNT = KNT + 1
          IFAIL = 1
          IF (KNT.LT.4) THEN
              GO TO 20
          ELSE
              GO TO 230
          END IF
      END IF
C Now retrieve the eigenvectors.
      DO 30 J = 1,N
          ALPHA = WORK(IALFR+J-1)
          BETA = WORK(IBETA+J-1)
          IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN
              GO TO 220

          ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN
C Form the eigenvector as Uz and normalise
              CALL MVN(WORK(IUL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
          ELSE

C Form the eigenvector Vz and normalise
              CALL MVN(WORK(IVL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
          END IF
   30 CONTINUE

C Finished retrieval of eigenvectors; now get phase angles. Remember the
C funny things that ATAN2 can do in IEEE arithmetic.
      DO 40 J = 1,N
          ALPHA = WORK(IALFR+J-1)
          BETA = WORK(IBETA+J-1)
          IF (ALPHA.EQ.0.0D0) THEN
              WORK(IALFR+J-1) = 0.0D0
          ELSE

              FAC = SIGN(1.0D0,ALPHA)
              ALPHA = FAC*ALPHA
              BETA = FAC*BETA
              WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA)
          END IF
   40 CONTINUE

C
C Now we are ready to form the initial matrix H. Remember that:
C the eigenvectors are in the matrix work(iv+j),j=0,..,n*n-1;
C the eigenvalues are in the elements work(ialfr+j),j=0,..,n.
C The following routine call will overwrite the eigenvectors but
C not the eigenvalues. Also, notice that work(iul+j), j=0,..,n*n-1
C is used here purely for storage.
C
      CALL UDUT(WORK(IV),N,WORK(IALFR),WORK(IHL),N,WORK(IULO),N)

C Do case of lower eigenvalue approximation first:

C Integrate forward, storing the matrix Sigma and the scalar argl
C at each of the sampling points, which are equally spaced between
C x=a and x=b:

      HSAMPL = (B-A)/DBLE(NSAMPL)
      X1 = A
      X2 = X1 + HSAMPL
      ISAMPR = 1
      ISAMPI = N*N + 1

      DO 60 JSAMPL = 1,NSAMPL - 1
          IFAIL = 0
          CALL HINT(X1,X2,ELAM,WORK(IHL),N,N,ADJUST,TOL,WORK(1),COFMAT,
     &              IFAIL)
          IF (IFAIL.NE.0) THEN
              RETURN
          ELSE

C Store the outcome of the integration:
              CALL SMCOPY('G',N,N,WORK(IHL),N,WSAMPL(1,1,JSAMPL),N)
              ARGL = 0.0D0
              INDEX = IHL
              DO 50 II = 1,N
                  ARGL = ARGL + WORK(INDEX)
                  INDEX = INDEX + N + 1
   50         CONTINUE
              ARGL = ARGL + ADJUST
              WORK(IASMP+JSAMPL-1) = ARGL
              ISAMPR = ISAMPR + N*N
              X1 = X2
              X2 = X2 + HSAMPL
          END IF
   60 CONTINUE

C Set up initial condition at x=b:
      CALL SYSBC(1,ISING,B,WORK(IUL),WORK(IVL),N,N,ELAM)

      CALL SMCOPY('G',N,N,WORK(IUL),N,WORK(IULO),N)
      CALL SMCOPY('G',N,N,WORK(IVL),N,WORK(IVLO),N)
   70 CONTINUE
      CALL FM2BJF(N,WORK(IULO),N,WORK(IVLO),N,0.0D0,WORK(IALFR),
     &            WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,
     &            WORK(1),IFAIL)
      IF (IFAIL.NE.0) THEN
          KNT = KNT + 1
          IFAIL = 1
          IF (KNT.LT.4) THEN
              GO TO 70
          ELSE
              GO TO 210
          END IF
      END IF
C Now retrieve the eigenvectors.
      DO 80 J = 1,N
          ALPHA = WORK(IALFR+J-1)
          BETA = WORK(IBETA+J-1)
          IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN
              GO TO 200

          ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN
C Form the eigenvector as Uz and normalise
              CALL MVN(WORK(IUL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
          ELSE

C Form the eigenvector Vz and normalise
              CALL MVN(WORK(IVL),N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
          END IF
   80 CONTINUE

C Finished retrieval of eigenvectors; now get phase angles. Remember the
C funny things that ATAN2 can do in IEEE arithmetic.
      DO 90 J = 1,N
          ALPHA = WORK(IALFR+J-1)
          BETA = WORK(IBETA+J-1)
          IF (ALPHA.EQ.0.0D0) THEN
              WORK(IALFR+J-1) = 0.0D0
          ELSE

              FAC = SIGN(1.0D0,ALPHA)
              ALPHA = FAC*ALPHA
              BETA = FAC*BETA
              WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA)
          END IF
   90 CONTINUE

C
      CALL UDUT(WORK(IV),N,WORK(IALFR),WORK(IHL),N,WORK(IULO),N)

C Integrate backwards:
      X1 = B
      HSAMPL = (A-B)/DBLE(NSAMPL)
      X2 = X1 + HSAMPL
      JMATCH = NSAMPL

      DO 160 JSAMPL = NSAMPL - 1,1,-1
          IFAIL = 0
          CALL HINT(X1,X2,ELAM,WORK(IHL),N,N,ADJUST,TOL,WORK(1),COFMAT,
     &              IFAIL)
          IF (IFAIL.NE.0) THEN
              RETURN
          ELSE
              INDEX = IHL
              ARGR = 0.0D0
              DO 100 II = 1,N
                  ARGR = ARGR + WORK(INDEX)
                  INDEX = INDEX + N + 1
  100         CONTINUE
              ARGR = ARGR + ADJUST
C
C Form the miss matrix at the current point:
C

              CALL SMCOPY('G',N,N,WSAMPL(1,1,JSAMPL),N,WORK(IUL),N)

              ISAMPR = JSAMPL*N*N + 1
              CALL SMLOAD('G',N,N,0.0D0,0.0D0,WORK(IVL),N)
              CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,WORK(IULO),
     &                    WORK(IVLO),N,WORK(1),CWORK,IFAIL)
              IF (IFAIL.NE.0) THEN
                  GO TO 190
              ELSE

C That does exp(i.hl)
                  INDEX = IHL
                  DO 120 J = 1,N
                      DO 110 I = 1,N
                          WORK(INDEX) = -WORK(INDEX)
                          INDEX = INDEX + 1
  110                 CONTINUE
  120             CONTINUE

                  CALL SMCOPY('G',N,N,WORK(IHL),N,WORK(IUL),N)
                  CALL SMLOAD('G',N,N,0.0D0,0.0D0,WORK(IVL),N)
                  CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,WORK(IUR),
     &                        WORK(IVR),N,WORK(1),CWORK,IFAIL)
                  IF (IFAIL.NE.0) THEN
                      GO TO 180
                  ELSE
C Now form
C exp(-i.hr).exp(i.hl) = (ur + i.vr)*(ul + i.vl)

                      CALL DGEMM('N','N',N,N,N,1.0D0,WORK(IUR),N,
     &                            WORK(IULO),N,0.0D0,WORK(IUL),N)
                      CALL DGEMM('N','N',N,N,N,-1.0D0,WORK(IVR),N,
     &                            WORK(IVLO),N,1.0D0,WORK(IUL),N)

                      CALL DGEMM('N','N',N,N,N,1.0D0,WORK(IVR),N,
     &                            WORK(IULO),N,0.0D0,WORK(IVL),N)
                      CALL DGEMM('N','N',N,N,N,1.0D0,WORK(IUR),N,
     &                            WORK(IVLO),N,1.0D0,WORK(IVL),N)

                      IFAIL = 1
                      CALL FM2AJF(WORK(IUL),N,WORK(IVL),N,N,WORK(1),
     &                            WORK(1+N),LWK,CWORK,IFAIL)
                      IF (IFAIL.NE.0) THEN
                          GO TO 170
                      ELSE
C Now we get the phase angles and store them in work(ialfr+j) for
C j=0,..,n-1.
                          ARGC = 0.0D0
                          ALFMAX = -100.0D0
                          ALFMIN = 100.0D0
                          DO 130 J = 1,N
                              ALPHA = WORK(N+J)
                              BETA = WORK(J)
                              OMEGA = ATAN2(ALPHA,BETA)
                              IF (OMEGA.LT.0.0D0) OMEGA = OMEGA + TWOPI
                              WORK(J) = OMEGA
                              ARGC = ARGC + OMEGA
                              IF (OMEGA.GT.ALFMAX) ALFMAX = OMEGA
                              IF (OMEGA.LT.ALFMIN) ALFMIN = OMEGA
  130                     CONTINUE

                          IF (ELAM.EQ.EU) DLA13F = ALFMIN/TWOPI
                          IF (ELAM.EQ.EL) DLA13F = ALFMAX/TWOPI - 1.0D0

                          IF (ELAM.EQ.EL) THEN
                              WORK(IDLSMP+JSAMPL-1) = DLA13F
                          ELSE
C              WRITE (11,FMT=*) 'Stored miss distance:',dla13f
C              WRITE (11,FMT=*) 'elam=',elam
C Compare with the miss-distance at the lower lambda and see if the difference
C is smaller than at other mesh-points.
C              WRITE (6,FMT=*) 'x,du-dl:',dble(jsampl)* (b-a)/
C     &          dble(nsampl) + a,dla13f - work(idlsmp+jsampl-1)
                              IF (ABS(DLA13F-WORK(IDLSMP+JSAMPL-1)).LE.
     &                            DIFF) THEN
                                  JMATCH = JSAMPL
                                  DIFF = ABS(DLA13F-
     &                                   WORK(IDLSMP+JSAMPL-1))
                              END IF

                          END IF

                          X1 = X2
                          X2 = X2 + HSAMPL
                          INDEX = IHL
                          DO 150 J = 1,N
                              DO 140 I = 1,N
                                  WORK(INDEX) = -WORK(INDEX)
                                  INDEX = INDEX + 1
  140                         CONTINUE
  150                     CONTINUE
                      END IF
                  END IF
              END IF
          END IF
  160 CONTINUE


      IF (ELAM.EQ.EL) THEN
          ELAM = EU
          GO TO 10
      END IF

      C = A + DBLE(JMATCH)* (B-A)/DBLE(NSAMPL)

C      WRITE (6,FMT=*) 'Match-point chosen=',c

      RETURN
  170 IFAIL = 9
      RETURN
  180 IFAIL = 10
      RETURN
  190 IFAIL = 11
      RETURN
  200 IFAIL = 15
      RETURN
  210 IFAIL = 14
      RETURN
  220 IFAIL = 13
      RETURN
  230 IFAIL = 12

C$st$ Unreachable comments ...








      END

C ------------------------------------------------------------------------


      SUBROUTINE MATSUM(TYPE,A,IA,FACA,B,IB,FACB,C,IC,FACC,N)

C     .. Scalar Arguments ..
      DOUBLE PRECISION FACA,FACB,FACC
      INTEGER IA,IB,IC,N
      CHARACTER TYPE
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(IA,N),B(IB,N),C(IC,N)
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
      IF (FACC.NE.0.0D0) THEN
          IF (TYPE.EQ.'S') THEN

              DO 20 J = 1,N
                  DO 10 I = J,N
                      C(I,J) = FACC*C(I,J) + FACA*A(I,J) + FACB*B(I,J)
                      C(J,I) = C(I,J)
   10             CONTINUE
   20         CONTINUE

          ELSE IF (TYPE.EQ.'A') THEN

              DO 40 J = 1,N
                  DO 30 I = J + 1,N
                      C(I,J) = FACC*C(I,J) + FACA*A(I,J) + FACB*B(I,J)
                      C(J,I) = -C(I,J)
   30             CONTINUE
                  C(J,J) = 0.D0
   40         CONTINUE
          ELSE


              DO 60 J = 1,N
                  DO 50 I = 1,N
                      C(I,J) = FACC*C(I,J) + FACA*A(I,J) + FACB*B(I,J)
   50             CONTINUE
   60         CONTINUE

          END IF


C Special coding to take account of the fact that the array c might be
C un-initialised in this case.
      ELSE IF (TYPE.EQ.'S') THEN
          DO 80 J = 1,N
              DO 70 I = J,N
                  C(I,J) = FACA*A(I,J) + FACB*B(I,J)
                  C(J,I) = C(I,J)
   70         CONTINUE
   80     CONTINUE
      ELSE IF (TYPE.EQ.'A') THEN
          DO 100 J = 1,N
              DO 90 I = J + 1,N
                  C(I,J) = FACA*A(I,J) + FACB*B(I,J)
                  C(J,I) = -C(I,J)
   90         CONTINUE
              C(J,J) = 0.D0
  100     CONTINUE
      ELSE
          DO 120 J = 1,N
              DO 110 I = 1,N
                  C(I,J) = FACA*A(I,J) + FACB*B(I,J)
  110         CONTINUE
  120     CONTINUE
      END IF


      END


C ------------------------------------------------------------------------

      SUBROUTINE HSTEP(X1,X2,TOL,FSAL,ELAM,H,IH,N,COFMAT,ERR,WORK,PRN,
     &                 IFAIL)
C     .. Parameters ..
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION ELAM,ERR,TOL,X1,X2
      INTEGER IFAIL,IH,N
      LOGICAL FSAL,PRN
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION H(IH,N),WORK(N* (2+18*N))
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL COFMAT
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ELMAX,FAC1,FAC10,FAC11,FAC12,FAC13,FAC14,FAC15,
     &                 FAC16,FAC17,FAC18,FAC19,FAC2,FAC20,FAC21,FAC22,
     &                 FAC23,FAC3,FAC4,FAC5,FAC6,FAC7,FAC8,FAC9,HSTP,HN
      INTEGER COFF11,COFF12,COFF21,COFF22,HDOT1,HDOT2,HDOT3,HDOT4,HDOT5,
     &        HDOT6,HINDEX,IE,IR,N2,S11,S12,S21,S22,STORE,V
C     .. Local Arrays ..
      DOUBLE PRECISION WRAF(1)
C     ..
C     .. External Subroutines ..
      EXTERNAL SMCOPY,HDER,MATSUM,MONIT
C     ..
C     .. External Functions ..

      DOUBLE PRECISION DLANGE
      EXTERNAL DLANGE
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,REAL
C     ..

      N2 = N*N
      HINDEX = 1
      HDOT1 = HINDEX + N2
      HDOT2 = HDOT1 + N2
      HDOT3 = HDOT2 + N2
      HDOT4 = HDOT3 + N2
      HDOT5 = HDOT4 + N2
      HDOT6 = HDOT5 + N2
      S11 = HDOT6 + N2
      S12 = S11 + N2
      S21 = S12 + N2
      S22 = S21 + N2
      COFF11 = S22 + N2
      COFF12 = COFF11 + N2
      COFF21 = COFF12 + N2
      COFF22 = COFF21 + N2
      STORE = COFF22 + N2
      V = STORE + N2
      IR = V + N2
      IE = IR + N
      IEND = N* (2+17*N) + 1
C Need workspace of length (18*n+2)*n
      HN = DLANGE('M',N,N,H,IH,WRAF)
C Store current H:
      CALL SMCOPY('G',N,N,H,IH,WORK(HINDEX),N)
C Get hdot1:
      IF (FSAL) THEN
          CALL SMCOPY('G',N,N,WORK(IEND),N,WORK(HDOT1),N)
      ELSE
          CALL HDER(X1,ELAM,WORK(HINDEX),N,WORK(HDOT1),N,N,WORK(S11),
     &              WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),
     &              WORK(COFF12),WORK(COFF21),WORK(COFF22),WORK(STORE),
     &              WORK(V),N,WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
          IF (IFAIL.NE.0) RETURN
      END IF
C Force a stepsize reduction if hdot is too large.
      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT1),N,WRAF)).GT.1.0D1) THEN
          ERR = TOL*32.D0
          RETURN
      ENDIF
C Get hdot2:
      HSTP = (X2-X1)/5.D0
      CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,HSTP,WORK(HINDEX),N,
     &            0.0D0,N)
      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT2),N,N,
     &          WORK(S11),WORK(S12),WORK(S21),WORK(S22),
     &          WORK(COFF11),WORK(COFF12),WORK(COFF21),WORK(COFF22),
     &          WORK(STORE),WORK(V),N,WORK(IR),WORK(IE),COFMAT,
     &          ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN

      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT2),N,WRAF)).GT.1.0D1) THEN
          ERR = TOL*32.D0
          RETURN
      END IF

C Get hdot3:
      HSTP = 3.D0* (X2-X1)/10.D0
      CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,0.25D0*HSTP,WORK(STORE),
     &             N,0.0D0,N)
      CALL MATSUM('S',WORK(STORE),N,1.0D0,WORK(HDOT2),N,HSTP*0.75D0,
     &             WORK(HINDEX),N,0.0D0,N)
      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT3),N,N,WORK(S11),
     &          WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12),
     &          WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N,
     &          WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN

      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT3),N,WRAF)).GT.1.0D1) THEN
          ERR = TOL*32.D0
          RETURN
      END IF

C Get hdot4:

      HSTP = (X2-X1)
      FAC1 = 44.D0/45.D0
      FAC2 = -56.D0/15.D0
      FAC3 = 32.D0/9.D0
      CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,FAC1*HSTP,WORK(HINDEX),
     &            N,0.0D0,N)
      CALL MATSUM('S',WORK(HDOT2),N,FAC2*HSTP,WORK(HDOT3),N,FAC3*HSTP,
     &             WORK(HINDEX),N,1.D0,N)
      HSTP = 0.8D0*HSTP
      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT4),N,N,WORK(S11),
     &          WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12),
     &          WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N,
     &          WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN

      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT4),N,WRAF)).GT.1.0D1) THEN
          ERR = TOL*32.D0
          RETURN
      END IF

C Get hdot5:

      HSTP = X2 - X1
      FAC4 = 19372.D0/6561.D0
      FAC5 = -25360.D0/2187.D0
      FAC6 = 64448.D0/6561.D0
      FAC7 = -212.D0/729.D0
      CALL SMCOPY('G',N,N,H,IH,WORK(HINDEX),N)
      CALL MATSUM('S',WORK(HDOT1),N,FAC4*HSTP,WORK(HDOT2),N,FAC5*HSTP,
     &            WORK(HINDEX),N,1.D0,N)
      CALL MATSUM('S',WORK(HDOT3),N,FAC6*HSTP,WORK(HDOT4),N,FAC7*HSTP,
     &            WORK(HINDEX),N,1.D0,N)
      HSTP = HSTP*8.D0/9.D0
      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT5),N,N,WORK(S11),
     &          WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12),
     &          WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N,
     &          WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN
      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT5),N,WRAF)).GT.1.0D1) THEN
          ERR = TOL*32.D0
          RETURN
      END IF

C Get hdot6:

      HSTP = X2 - X1
      FAC8 = 9017.D0/3168.D0
      FAC9 = -355.D0/33.D0
      FAC10 = 46732.D0/5247.D0
      FAC11 = 49.D0/176.D0
      FAC12 = -5103.D0/18656.D0
      CALL SMCOPY('G',N,N,H,IH,WORK(HINDEX),N)
      CALL MATSUM('S',WORK(HDOT1),N,FAC8*HSTP,WORK(HDOT2),N,FAC9*HSTP,
     &             WORK(HINDEX),N,1.D0,N)
      CALL MATSUM('S',WORK(HDOT3),N,FAC10*HSTP,WORK(HDOT4),N,FAC11*HSTP,
     &             WORK(HINDEX),N,1.D0,N)
      CALL MATSUM('S',WORK(HINDEX),N,1.D0,WORK(HDOT5),N,FAC12*HSTP,
     &             WORK(S11),N,0.D0,N)
      CALL SMCOPY('G',N,N,WORK(S11),N,WORK(HINDEX),N)

      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT6),N,N,WORK(S11),
     &          WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),
     &          WORK(COFF12),WORK(COFF21),WORK(COFF22),WORK(STORE),
     &          WORK(V),N,WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN

C Bug fixed on 8th September 1997

      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT6),N,WRAF)).GT.1.0D1) THEN
           ERR = TOL*32.D0
           RETURN
      END IF

C Get hdot7 (store it in WORK(HDOT2) which is no longer required
C at this point):


      HSTP = X2 - X1
      FAC13 = 35.D0/384.D0
      FAC14 = 500.D0/1113.D0
      FAC15 = 125.D0/192.D0
      FAC16 = -2187.D0/6784.D0
      FAC17 = 11.D0/84.D0
      CALL MATSUM('S',H,IH,1.0D0,WORK(HDOT1),N,FAC13*HSTP,WORK(HINDEX),
     &             N,0.0D0,N)
      CALL MATSUM('S',WORK(HDOT3),N,FAC14*HSTP,WORK(HDOT4),N,FAC15*HSTP,
     &             WORK(HINDEX),N,1.0D0,N)
      CALL MATSUM('S',WORK(HDOT5),N,FAC16*HSTP,WORK(HDOT6),N,FAC17*HSTP,
     &             WORK(HINDEX),N,1.0D0,N)

      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT2),N,N,WORK(S11),
     &          WORK(S12),WORK(S21),WORK(S22),WORK(COFF11),WORK(COFF12),
     &          WORK(COFF21),WORK(COFF22),WORK(STORE),WORK(V),N,
     &          WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN

C Bug fixed on 8th September 1997

      IF (ABS((X2-X1)*DLANGE('M',N,N,WORK(HDOT2),N,WRAF)).GT.1.0D1) THEN
          ERR = TOL*32.D0
          RETURN
      END IF

C Now store the old value of H:
      CALL SMCOPY('G',N,N,H,IH,WORK(S11),N)
C Now update h:
      CALL SMCOPY('G',N,N,WORK(HINDEX),N,H,IH)
C Now store the old value of H:
      CALL SMCOPY('G',N,N,WORK(S11),N,WORK(HINDEX),N)

C Now form the error estimate
      HSTP = (X2-X1)
      FAC18 = 5179.D0/57600.D0
      FAC19 = 7571.D0/16695.D0
      FAC20 = 393.D0/640.D0
      FAC21 = -92097.D0/339200.D0
      FAC22 = 187.D0/2100.D0
      FAC23 = 1.D0/40.D0
      CALL MATSUM('S',WORK(HDOT1),N,FAC18*HSTP,WORK(HDOT3),N,FAC19*HSTP,
     &             WORK(HINDEX),N,1.D0,N)
      CALL MATSUM('S',WORK(HDOT4),N,FAC20*HSTP,WORK(HDOT5),N,FAC21*HSTP,
     &             WORK(HINDEX),N,1.D0,N)
      CALL MATSUM('S',WORK(HDOT6),N,FAC22*HSTP,WORK(HDOT2),N,FAC23*HSTP,
     &             WORK(HINDEX),N,1.D0,N)

      CALL HDER(X1+HSTP,ELAM,WORK(HINDEX),N,WORK(HDOT3),N,N,
     &          WORK(S11),WORK(S12),WORK(S21),WORK(S22),
     &          WORK(COFF11),WORK(COFF12),WORK(COFF21),
     &          WORK(COFF22),WORK(STORE),WORK(V),N,
     &          WORK(IR),WORK(IE),COFMAT,ELMAX,IFAIL)
      IF (IFAIL.NE.0) RETURN

      CALL MATSUM('S',WORK(HINDEX),N,1.D0,H,IH,-1.D0,WORK(S12),N,0.D0,N)
      CALL MATSUM('S',WORK(HDOT2),N,1.D0,WORK(HDOT3),IH,-1.D0,
     &             WORK(HDOT1),N,0.D0,N)

      ERR = SQRT( (DLANGE('M',N,N,WORK(S12),N,WRAF)**2
     &     + DLANGE('M',N,N,WORK(HDOT1),N,WRAF)**2)/(1.D0+HN**2) )

      CALL SMCOPY('G',N,N,WORK(HDOT2),N,WORK(IEND),N)

C Now return to the calling routine
      RETURN

      END



C ------------------------------------------------------------------------

      SUBROUTINE MATEXP(AR,AI,IA,N,ISWTCH,EXPR,EXPI,IEXP,WORK,CWORK,
     &                  IFAIL)
C This subroutine takes in a matrix ar+i*ai, which is either Hermitian or
C skew-Hermitian, diagonalises it, and hence computes its exponential.

C     .. Scalar Arguments ..
      INTEGER IA,IEXP,IFAIL,ISWTCH,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION AI(IA,N),AR(IA,N),EXPI(IEXP,N),EXPR(IEXP,N),
     &                 WORK(4*N* (N+1)),CWORK(N*(N+65))
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION C,MAXEL,S
      INTEGER I,INDEX,IR,IVI,IVIS,IVR,IVRS,IWK1,IWK2,IWK3,J
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION WRAF(1)

C     ..
C     .. External Subroutines ..
      EXTERNAL FM2AXF,DGEMM
C     ..
C     .. External Functions ..
      DOUBLE PRECISION DLANGE
      EXTERNAL DLANGE
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC COS,EXP,MAX,SIN
C     ..

      MAXEL = MAX(1.0D0,DLANGE('M',N,N,AR,IA,WRAF))
      MAXEL = MAX(MAXEL,DLANGE('M',N,N,AI,IA,WRAF))

      IF (ISWTCH.EQ.1) THEN
C The matrix is Hermitian and has real eigenvalues
          IR = 1
          IVR = N + 1
          IVI = IVR + N*N
          IWK1 = IVI + N*N
          IWK2 = IWK1 + N
          IWK3 = IWK2 + N
          IVRS = IWK3 + N
          IVIS = IVRS + N*N

          DO 20 J = 1,N
              DO 10 I = 1,N
                  EXPR(I,J) = AR(I,J)/MAXEL
                  EXPI(I,J) = AI(I,J)/MAXEL
   10         CONTINUE
   20     CONTINUE

          IFAIL = -1
          CALL FM2AXF(EXPR,IEXP,EXPI,IEXP,N,WORK(IR),WORK(IVR),N,
     &                WORK(IVI),N,WORK(IWK1),WORK(IWK2),WORK(IWK3),
     &                CWORK,IFAIL)
          IF (IFAIL.NE.0) THEN
              IFAIL = 1
              RETURN
          ELSE
C Now get the exponential of the eigenvalues:
              INDEX = IR
              DO 30 J = 1,N
                  WORK(INDEX) = EXP(MAXEL*WORK(INDEX))
                  INDEX = INDEX + 1
   30         CONTINUE
C Now the matrix exponential is given by
C
C     expr + i*expi = (rr+i*ri)*exp(D)*(transpose(rr)-i*transpose(ri))
C

              INDEX = 0
              DO 50 J = 1,N
                  DO 40 I = 1,N
                      WORK(IVRS+INDEX) = WORK(IVR+INDEX)*WORK(IR+J-1)
                      WORK(IVIS+INDEX) = WORK(IVI+INDEX)*WORK(IR+J-1)
                      INDEX = INDEX + 1
   40             CONTINUE
   50         CONTINUE

              CALL DGEMM('N','T',N,N,N,1.0D0,WORK(IVR),N,WORK(IVRS),N,
     &                    0.0D0,EXPR,IEXP)
              CALL DGEMM('N','T',N,N,N,1.0D0,WORK(IVI),N,WORK(IVIS),N,
     &                    1.0D0,EXPR,IEXP)

              CALL DGEMM('N','T',N,N,N,1.0D0,WORK(IVI),N,WORK(IVRS),N,
     &                    0.0D0,EXPI,IEXP)
              CALL DGEMM('N','T',N,N,N,-1.0D0,WORK(IVR),N,WORK(IVIS),N,
     &                    1.0D0,EXPI,IEXP)
          END IF
      ELSE


C     expr + i*expi = (rr+i*ri)*exp(D)*(transpose(rr)-i*transpose(ri))


C In this case the matrix A is anti-Hermitian. This means that -iA is
C Hermitian. Note that -iA = ai - i*ar. We find the eigenvalues of -iA,
C which are all real, and translate them back into the eigenvalues of A
C which follow from A = i(-iA).

C Negate all the elements of ar:

          DO 70 J = 1,N
              DO 60 I = 1,N
                  EXPR(I,J) = -AR(I,J)/MAXEL
                  EXPI(I,J) = AI(I,J)/MAXEL
   60         CONTINUE

              AR(J,J) = 0.0D0
              EXPR(J,J) = 0.0D0
   70     CONTINUE

C Get eigenvalues:

          IR = 1
          IVR = N + 1
          IVI = IVR + N*N
          IWK1 = IVI + N*N
          IWK2 = IWK1 + N
          IWK3 = IWK2 + N
          IVRS = IWK3 + N
          IVIS = IVRS + N*N
          IFAIL = -1
          CALL FM2AXF(EXPI,IEXP,EXPR,IEXP,N,WORK(IR),WORK(IVR),N,
     &                WORK(IVI),N,WORK(IWK1),WORK(IWK2),WORK(IWK3),
     &                CWORK,IFAIL)
          IF (IFAIL.NE.0) THEN
              IFAIL = 1
              RETURN
          ELSE

CC Now we can return the matrix A to its original form:
CC         DO 70 j = 1,n
C              DO 60 i = 1,n
C                  ar(i,j) = -ar(i,j)
C60            CONTINUE
C70        CONTINUE

C This time we need both sines and cosines of the eigenvalues in order to form
C the matrix exponential.

              INDEX = 0
              DO 90 J = 1,N
                  DO 80 I = 1,N
                      C = COS(MAXEL*WORK(IR+J-1))
                      S = SIN(MAXEL*WORK(IR+J-1))
                      WORK(IVRS+INDEX) = WORK(IVR+INDEX)*C +
     &                                   WORK(IVI+INDEX)*S
                      WORK(IVIS+INDEX) = WORK(IVR+INDEX)*S -
     &                                   WORK(IVI+INDEX)*C
                      INDEX = INDEX + 1
   80             CONTINUE
   90         CONTINUE

              CALL DGEMM('N','T',N,N,N,1.0D0,WORK(IVR),N,WORK(IVRS),N,
     &                    0.0D0,EXPR,IEXP)
              CALL DGEMM('N','T',N,N,N,-1.0D0,WORK(IVI),N,WORK(IVIS),N,
     &                    1.0D0,EXPR,IEXP)

              CALL DGEMM('N','T',N,N,N,1.0D0,WORK(IVR),N,WORK(IVIS),N,
     &                    0.0D0,EXPI,IEXP)
              CALL DGEMM('N','T',N,N,N,1.0D0,WORK(IVI),N,WORK(IVRS),N,
     &                    1.0D0,EXPI,IEXP)

          END IF
      END IF

      IFAIL = 0

C$st$ Unreachable comments ...


      END


C -------------------------------------------------------------------------

      DOUBLE PRECISION FUNCTION DLH13F(A,B,C,ELAM,N,TOL,HL,IHL,HR,IHR,
     &                 UL,VL,UR,VR,WORK,LWK,CWORK,COFMAT,SYSBC,SHIFT,
     &                 IFAIL)
C This subroutine computes the M(lambda) miss distance function
C     .. Scalar Arguments ..
      DOUBLE PRECISION A,B,C,ELAM,SHIFT,TOL
      INTEGER IFAIL,IHL,IHR,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION HL(IHL,N),HR(IHR,N),UL(N,N),UR(N,N),VL(N,N),
     &                 VR(N,N),WORK(N* (13+21*N)+25)
      INTEGER LWK(N)
      COMPLEX*16 CWORK(N*(N+65))
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL COFMAT,SYSBC
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ADJL,ADJR,ALFMAX,ALFMIN,ALPHA,ARGL,ARGR,BETA,FAC,
     &                 OMEGA,SUM,TWOPI
      INTEGER I,IALFI,IALFR,IBETA,IUL,IV,IVL,IWK,J,KNT,N2
      LOGICAL ISING
C     ..
C     .. External Functions ..
      DOUBLE PRECISION XM1AAF
      EXTERNAL XM1AAF
C     ..
C     .. External Subroutines ..
      EXTERNAL FM2AJF,FM2BJF,SMCOPY,SMLOAD,DGEMM,HINT,MATEXP,MVN,UDUT
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN2,DBLE,SIGN
C     ..
      TWOPI = 2.0D0*XM1AAF(0.0D0)
      N2 = N*N
C      WRITE (6,FMT=*) '**** In dlh13f, elam=',elam
C Set pointers for workspace array:
      IUL = 1
      IVL = IUL + N2
      IALFR = IVL + N2
      IALFI = IALFR + N
      IBETA = IALFI + N
      IV = IBETA + N
      IWK = IV + N2

C Set up boundary condition at x=a:
      CALL SYSBC(0,ISING,A,UL,VL,N,N,ELAM)
C Now we must transform this into a boundary condition on H by getting
C the eigenvalues and eigenvectors of the Theta-matrix associated with
C these boundary conditions. We do this by solving the simultaneous
C eigenvalue problem
C (cos(omega/2)U-sin(omega/2)V)z = 0,
C then recovering the eigenvectors w from a formula such as
C w = (V-iU)z.
C
      CALL SMCOPY('G',N,N,UL,N,WORK(IUL),N)
      CALL SMCOPY('G',N,N,VL,N,WORK(IVL),N)
      IFAIL = 1
      KNT = 0
   10 CONTINUE

      CALL FM2BJF(N,WORK(IUL),N,WORK(IVL),N,0.0D0,WORK(IALFR),
     &            WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,
     &            WORK(IWK),IFAIL)
      IF (IFAIL.NE.0) THEN
          KNT = KNT + 1
          IFAIL = 1
          IF (KNT.LT.4) THEN
              GO TO 10
          ELSE
              GO TO 140
          END IF
      END IF
C Now retrieve the eigenvectors.
      DO 20 J = 1,N
          ALPHA = WORK(IALFR+J-1)
          BETA = WORK(IBETA+J-1)
          IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN
              GO TO 130

          ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN
C Form the eigenvector as Uz and normalise
              CALL MVN(UL,N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
          ELSE

C Form the eigenvector Vz and normalise
              CALL MVN(VL,N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
          END IF
   20 CONTINUE

C Finished retrieval of eigenvectors; now get phase angles. Remember the
C funny things that ATAN2 can do in IEEE arithmetic.
      DO 30 J = 1,N
          ALPHA = WORK(IALFR+J-1)
          BETA = WORK(IBETA+J-1)
          IF (ALPHA.EQ.0.0D0) THEN
              WORK(IALFR+J-1) = 0.0D0
          ELSE

              FAC = SIGN(1.0D0,ALPHA)
              ALPHA = FAC*ALPHA
              BETA = FAC*BETA
              WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA)
          END IF
   30 CONTINUE

C
C Now we are ready to form the initial matrix H. Remember that:
C the eigenvectors are in the matrix work(iv+j),j=0,..,n*n-1;
C the eigenvalues are in the elements work(ialfr+j),j=0,..,n.
C The following routine call will overwrite the eigenvectors but
C not the eigenvalues. Also, notice that work(iul+j), j=0,..,n*n-1
C is used here purely for storage.
C
      CALL UDUT(WORK(IV),N,WORK(IALFR),HL,IHL,WORK(IUL),N)

C Now do the integration from x=a to x=c:

      CALL HINT(A,C,ELAM,HL,IHL,N,ADJL,TOL,WORK(1),COFMAT,IFAIL)
      IF (IFAIL.EQ.0) THEN

C That does the left-hand leg of the integration. Now do the right
C hand leg.

C Set up boundary condition at x=b:
          CALL SYSBC(1,ISING,B,UR,VR,N,N,ELAM)
C Now we must transform this into a boundary condition on H by getting
C the eigenvalues and eigenvectors of the Theta-matrix associated with
C these boundary conditions. We do this by solving the simultaneous
C eigenvalue problem
C (cos(omega/2)U-sin(omega/2)V)z = 0,
C then recovering the eigenvectors w from a formula such as
C w = (V-iU)z.
C
          IFAIL = 1
          KNT = 0
          CALL SMCOPY('G',N,N,UR,N,WORK(IUL),N)
          CALL SMCOPY('G',N,N,VR,N,WORK(IVL),N)
   40     CONTINUE

          CALL FM2BJF(N,WORK(IUL),N,WORK(IVL),N,0.0D0,WORK(IALFR),
     &                WORK(IALFI),WORK(IBETA),.true.,WORK(IV),N,LWK,
     &                WORK(IWK),IFAIL)
          IF (IFAIL.NE.0) THEN
              KNT = KNT + 1
              IFAIL = 1
              IF (KNT.LT.4) THEN
                  GO TO 40
              ELSE
                  GO TO 120
              END IF
          END IF
C Now retrieve the eigenvectors.
          DO 50 J = 1,N
              ALPHA = WORK(IALFR+J-1)
              BETA = WORK(IBETA+J-1)
              IF (ABS(ALPHA)+ABS(BETA).EQ.0.0D0) THEN
                  GO TO 110

              ELSE IF (ABS(ALPHA).GT.ABS(BETA)) THEN
C Form the eigenvector as Uz and normalise
                  CALL MVN(UR,N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
              ELSE

C Form the eigenvector Vz and normalise
                  CALL MVN(VR,N,WORK(IV+ (J-1)*N),WORK(IALFI),N)
              END IF
   50     CONTINUE

C Finished retrieval of eigenvectors; now get phase angles. Remember the
C funny things that ATAN2 can do in IEEE arithmetic.
          DO 60 J = 1,N
              ALPHA = WORK(IALFR+J-1)
              BETA = WORK(IBETA+J-1)
              IF (ALPHA.EQ.0.0D0) THEN
                  WORK(IALFR+J-1) = TWOPI
              ELSE

                  FAC = SIGN(1.0D0,ALPHA)
                  ALPHA = FAC*ALPHA
                  BETA = FAC*BETA
                  WORK(IALFR+J-1) = 2.0D0*ATAN2(ALPHA,BETA)
              END IF
   60     CONTINUE
C
C Now we are ready to form the initial matrix H. Remember that:
C the eigenvectors are in the matrix work(iv+j),j=0,..,n*n-1;
C the eigenvalues are in the elements work(ialfr+j),j=0,..,n.
C The following routine call will overwrite the eigenvectors but
C not the eigenvalues. Also, notice that work(iul+j), j=0,..,n*n-1
C is used here purely for storage.
C
          CALL UDUT(WORK(IV),N,WORK(IALFR),HR,IHR,WORK(IUL),N)

C Now do the integration from x=b to x=c:
          CALL HINT(B,C,ELAM,HR,IHR,N,ADJR,TOL,WORK(1),COFMAT,IFAIL)
          IF (IFAIL.EQ.0) THEN
C
C Now we have the left-hand and right-hand H-matrices at the centre.
C We must compute the phase-angles of the matrix exp(-i.hr).exp(i.hl).
C First we must compute the two matrices concerned, since the eigenvalues
C are not those of the matrix hl-hr.
C First get the traces since we will probably overwrite the matrices
C hl and hr and we need to know their traces.
              ARGL = ADJL
              ARGR = ADJR
              DO 70 I = 1,N
                  ARGL = ARGL + HL(I,I)
                  ARGR = ARGR + HR(I,I)
   70         CONTINUE

              CALL SMCOPY('G',N,N,HL,IHL,WORK(IUL),N)
              CALL SMLOAD('G',N,N,0.0D0,0.0D0,WORK(IVL),N)
              CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,UL,VL,N,WORK(IWK),
     &                    CWORK,IFAIL)
              IF (IFAIL.NE.0) THEN
                  IFAIL = 16
              ELSE

C That does exp(i.hl)

                  DO 90 J = 1,N
                      DO 80 I = 1,N
                          HR(I,J) = -HR(I,J)
   80                 CONTINUE
   90             CONTINUE

                  CALL SMCOPY('G',N,N,HR,IHR,WORK(IUL),N)
                  CALL SMLOAD('G',N,N,0.0D0,0.0D0,WORK(IVL),N)
                  CALL MATEXP(WORK(IVL),WORK(IUL),N,N,2,UR,VR,N,
     &                        WORK(IWK),CWORK,IFAIL)
                  IF (IFAIL.NE.0) THEN
                      IFAIL = 17
                  ELSE
C Now form
C exp(-i.hr).exp(i.hl) = (ur + i.vr)*(ul + i.vl)

                      CALL DGEMM('N','N',N,N,N,1.0D0,UR,N,UL,N,0.0D0,
     &                            WORK(IUL),N)
                      CALL DGEMM('N','N',N,N,N,-1.0D0,VR,N,VL,N,1.0D0,
     &                            WORK(IUL),N)

                      CALL DGEMM('N','N',N,N,N,1.0D0,VR,N,UL,N,0.0D0,
     &                            WORK(IVL),N)
                      CALL DGEMM('N','N',N,N,N,1.0D0,UR,N,VL,N,1.0D0,
     &                            WORK(IVL),N)

C We have got the real and imaginary parts; now we can find the
C eigenvalues. For this we use the NAG routine f02ajf

                      IFAIL = 1
                      CALL FM2AJF(WORK(IUL),N,WORK(IVL),N,N,WORK(IWK),
     &                            WORK(IWK+N),LWK,CWORK,IFAIL)
                      IF (IFAIL.NE.0) THEN
                          IFAIL = 18
                      ELSE
C Now we get the phase angles and store them in work(ialfr+j) for
C j=0,..,n-1.
                          SUM = 0.0D0
                          ALFMAX = -100.0D0
                          ALFMIN = 100.0D0
                          DO 100 J = 1,N
                              ALPHA = WORK(IWK+N+J-1)
                              BETA = WORK(IWK+J-1)
                              OMEGA = ATAN2(ALPHA,BETA)
                              IF (OMEGA.LT.0.0D0) OMEGA = OMEGA + TWOPI
                              WORK(IALFR+J-1) = OMEGA
                              SUM = SUM + OMEGA
                              IF (OMEGA.GT.ALFMAX) ALFMAX = OMEGA
                              IF (OMEGA.LT.ALFMIN) ALFMIN = OMEGA
  100                     CONTINUE

C      WRITE (6,FMT=*) 'argl,argr,argc:',argl,argr,sum

                          DLH13F = (ARGL-ARGR-SUM)/TWOPI + DBLE(N) -
     &                             SHIFT

C      WRITE (6,FMT=*) ' From DLH13F: ELAM,DMISS:',ELAM,DLH13F
                          IF (ABS(DLH13F).LE.0.6D0) THEN

                              IF (DLH13F.GT.0.0D0) THEN
                                  DLH13F = ALFMIN/TWOPI
                              ELSE

                                  DLH13F = ALFMAX/TWOPI - 1.0D0
                              END IF

                          END IF
                      END IF
                  END IF
              END IF
          END IF
          RETURN
  110     IFAIL = 15
          RETURN
  120     IFAIL = 14
      END IF
      RETURN
  130 IFAIL = 13
      RETURN
  140 IFAIL = 12

C$st$ Unreachable comments ...







      END

C ------------------------------------------------------------------------

      SUBROUTINE UDUT(U,IU,D,OUT,IOUT,STORE,N)
C     .. Scalar Arguments ..
      INTEGER IOUT,IU,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION D(N),OUT(IOUT,N),STORE(N,N),U(IU,N)
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
C     .. External Subroutines ..
      EXTERNAL DGEMM
C     ..
      DO 20 J = 1,N
          DO 10 I = 1,N
              STORE(I,J) = U(I,J)*D(J)
   10     CONTINUE
   20 CONTINUE
      CALL DGEMM('N','T',N,N,N,1.0D0,STORE,N,U,IU,0.0D0,OUT,IOUT)

      END

C -------------------------------------------------------------------------

      SUBROUTINE MVN(A,IA,V,STORE,N)
C     .. Scalar Arguments ..
      INTEGER IA,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(IA,N),STORE(N),V(N)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ELMAX,SUM
      INTEGER I,J
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX,SQRT
C     ..
      DO 10 I = 1,N
          STORE(I) = 0.0D0
   10 CONTINUE
      DO 30 J = 1,N
          DO 20 I = 1,N
              STORE(I) = STORE(I) + A(I,J)*V(J)
   20     CONTINUE
   30 CONTINUE
      ELMAX = 0.0D0
      DO 40 I = 1,N
          ELMAX = MAX(ELMAX,ABS(STORE(I)))
   40 CONTINUE
      SUM = 0.0D0
      DO 50 I = 1,N
          STORE(I) = STORE(I)/ELMAX
          SUM = SUM + STORE(I)**2
   50 CONTINUE
      SUM = SQRT(SUM)
      DO 60 I = 1,N
          V(I) = STORE(I)/SUM
   60 CONTINUE


      END

C -------------------------------------------------------------------------

      SUBROUTINE HINT(XO,XEND,ELAM,H,IH,N,ADJUST,TOL,WORK,COFMAT,IFAIL)
C     .. Scalar Arguments ..
      DOUBLE PRECISION ADJUST,ELAM,TOL,XEND,XO
      INTEGER IFAIL,IH,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION H(IH,N),WORK(N* (21*N+13)+25)
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL COFMAT
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ALFA,BETA,BND,ERR,GAP,GAPL,HL,HMAX,PI,RATIO,
     &                 TOLOC,TWOPI,U,X1,X2
      INTEGER I,IE,IERR,IFLAG,IGAP,IGOOD,IHLOC,INDEX,IR,ITRY,IV,IWORK,
     &        NSTEPS
      LOGICAL PHASE1,PRN,FSAL
C     ..
C     .. External Functions ..
      DOUBLE PRECISION XM1AAF,XM2AJF
      EXTERNAL XM1AAF,XM2AJF
C     ..
C     .. Parameters ..
      DOUBLE PRECISION SAFE,SMALL,FIFTH
      PARAMETER (SAFE=0.80D0,SMALL=0.70D0,FIFTH=0.20D0)
C     ..
C     .. External Subroutines ..
      EXTERNAL FM2ABF,SMCOPY,HSTEP,MM1CAF,UDUT
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,DBLE,MAX,MIN,NINT,SIGN
C     ..
      PRN = .false.
      FSAL = .FALSE.
      PI = XM1AAF(0.0D0)
      TWOPI = 2.0D0*XM1AAF(0.0D0)
      HL = TOL*SIGN(1.0D0,XEND-XO)
      HMAX = ABS(XEND-XO)/4.0D0
      ADJUST = 0.0D0
      TOLOC = MIN(0.1D0,TOL)
      NSTEPS = 0
      ITRY = 0
      IERR = 0
      IGOOD = 4
      IHLOC = 1

      IE = IHLOC + N*N
      IV = IE + N
      IR = IV + N*N

      IWORK = IHLOC + N*N
      PHASE1 = .true.
C Phase1 controls the special step-control algorithm which is used only
C on the first step to get the meshing `on scale'.
      X1 = XO
   10 CONTINUE
      CALL SMCOPY('G',N,N,H,IH,WORK(IHLOC),N)
C Before doing any stepping we subtract off a matrix M which commutes with
C H and is such that exp(iM) = I. This will avoid singularities later.

      IFLAG = 1
      CALL FM2ABF(WORK(IHLOC),N,N,WORK(IE),WORK(IV),N,WORK(IR),IFLAG)
      IF (IFLAG.NE.0) THEN
          GO TO 70
      ELSE

          CALL SMCOPY('G',N,1,WORK(IE),N,WORK(IR),N)
          INDEX = IR
          DO 20 I = 1,N
              ALFA = TWOPI*DBLE(NINT(WORK(INDEX)/TWOPI))
              WORK(INDEX) = WORK(INDEX) - ALFA
              IF (WORK(INDEX).LT.0.0D0) WORK(INDEX) = WORK(INDEX) +
     &            TWOPI
              INDEX = INDEX + 1
   20     CONTINUE

C Now order the eigenvalues
          IFLAG = 0
          CALL MM1CAF(WORK(IR),1,N,'A',IFLAG)

C Now find the greatest gap

          GAP = -WORK(IR+N-1) + WORK(IR) + TWOPI
          IGAP = N

          DO 30 I = 1,N - 1
              GAPL = WORK(IR+I) - WORK(IR+I-1)
              IF (GAPL.GT.GAP) THEN
                  GAP = GAPL
                  IGAP = I
              END IF
   30     CONTINUE

C Have now found the greatest gap.
          IF (IGAP.LT.N) THEN
              BETA = 0.50D0* (WORK(IR+IGAP)+WORK(IR+IGAP-1))
          ELSE
              BETA = 0.50D0* (WORK(IR)+WORK(IR+N-1)-TWOPI)
          END IF
          IF (BETA.GE.PI) BETA = BETA - TWOPI

C Now normalise eigenvalues so that those which correspond to close points on the
C unit circle are indeed close
          DO 40 I = 1,N
              ALFA = TWOPI*DBLE(NINT(WORK(IE+I-1)/TWOPI))
              U = WORK(IE+I-1) - ALFA
              IF (U.GE.TWOPI+BETA) THEN
                  U = U - TWOPI
              ELSE IF (U.LT.BETA) THEN
                  U = U + TWOPI
              END IF
              ADJUST = ADJUST + WORK(IE+I-1) - U
              WORK(IE+I-1) = U
   40     CONTINUE
C          WRITE (6,FMT=*) 'Eigenvalue:',I,WORK(IE+I-1)

C Now construct the new matrix H:

          CALL UDUT(WORK(IV),N,WORK(IE),WORK(IHLOC),N,WORK(IR),N)

          CALL SMCOPY('G',N,N,WORK(IHLOC),N,H,IH)

C Use the routine hstep to do the stepping
          X2 = X1 + HL

C      WRITE (6,FMT=*) 'H X1 X2 ELAM',X2-X1,X1,X2,ELAM

          IF (ABS(XEND-X1).LT.1.D1*XM2AJF(0.0D0)) THEN
              GO TO 60
          ELSE
C
C WARNING: The last N**2 elements of the array WORK must be
C preserved between one call to HSTEP and the next.
C
              CALL HSTEP(X1,X2,TOLOC,FSAL,ELAM,WORK(IHLOC),N,N,COFMAT,
     &                   ERR,WORK(IWORK),PRN,IFAIL)

              IF (IFAIL.NE.0) THEN
                  RETURN
              ELSE
                  FSAL = .FALSE.
                  IF (ERR.GT.TOLOC) THEN
C          WRITE (7,FMT=*) 'Failed step'
                      IERR = IERR + 1
                      IGOOD = 0
                      IF (IERR.GT.16 .AND. .NOT.PHASE1) THEN
                          GO TO 50
                      ELSE IF (.NOT.PHASE1 .OR. IERR.LE.64) THEN
                          RATIO = MIN((TOLOC/ERR)**FIFTH,SAFE)
                          HL = HL*RATIO
                          GO TO 10
                      END IF
                  ELSE

                      IF (ERR.GT.TOLOC*SAFE) THEN
                          RATIO = (TOLOC*SAFE/ERR)**FIFTH
                          HL = HL*RATIO
                      ELSE IF (ERR.LT.SMALL*TOLOC .AND. IERR.EQ.0) THEN
                          IGOOD = IGOOD + 1
                          RATIO = (ERR/TOLOC)**FIFTH
                          BND = MAX(FIFTH,1.0D0/ (DBLE(IGOOD+1)))
                          HL = HL/MAX(BND,RATIO)
                          IF (ABS(HL).GT.HMAX) HL = SIGN(1.0D0,HL)*HMAX
                          IF (PHASE1 .AND. ITRY.LT.8 .AND.
     &                        ABS(HL).LT.HMAX) THEN
                              ITRY = ITRY + 1
C              WRITE (7,FMT=*) 'Failed step'
                              GO TO 10
                          END IF
                      END IF

                      FSAL = .TRUE.
C      WRITE (7,FMT=*) 'Successful step'
                      IERR = 0
                      PHASE1 = .false.
                      IF (ABS(X2-X1).GE.1.0D-3) PRN = .false.
                      X1 = X2

                      IF ((XEND-X2-HL)*SIGN(1.0D0,HL).LE.
     &                    0.0D0) HL = XEND - X1
C Update h
                      CALL SMCOPY('G',N,N,WORK(IHLOC),N,H,IH)
                      NSTEPS = NSTEPS + 1
                      GO TO 10
                  END IF
              END IF
          END IF
      END IF
      IFAIL = 21
      RETURN
   50 IFAIL = 19
      RETURN
C          WRITE (6,FMT=*) 'No. of steps:',nsteps
   60 IFAIL = 0
      RETURN
   70 IFAIL = 20

      END


C ----------------------------------------------------------------------

      SUBROUTINE HDER(X,ELAM,H,IH,HDOT,IHDOT,N,S11,S12,S21,S22,COFF11,
     &                COFF12,COFF21,COFF22,STORE,V,IV,R,E,COFMAT,ELMAX,
     &                IFAIL)
C     .. Scalar Arguments ..
      DOUBLE PRECISION ELAM,ELMAX,X
      INTEGER IFAIL,IH,IHDOT,IV,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION COFF11(N,N),COFF12(N,N),COFF21(N,N),COFF22(N,N),
     &                 E(N),H(IH,N),HDOT(IHDOT,N),R(N),S11(N,N),
     &                 S12(N,N),S21(N,N),S22(N,N),STORE(N,N),V(IV,N)
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL COFMAT
C     ..
C     .. External Functions ..
      DOUBLE PRECISION XM1AAF
      EXTERNAL XM1AAF
C     ..
C     .. External Subroutines ..
      EXTERNAL FM2ABF,SMCOPY,DGEMM,FGH
C     ..

C Step 1: find the eigenvalues and eigenvectors of H:
C     .. Local Scalars ..
      DOUBLE PRECISION ELOC,TWOPI,F,G,GT,HL
      INTEGER I,IFLAG,J,P,Q
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX
C     ..
      IFLAG = 1
      CALL FM2ABF(H,IH,N,E,V,IV,R,IFLAG)
      IF (IFLAG.NE.0) THEN
          IFAIL = 20
      ELSE
          TWOPI = 2.0D0*XM1AAF(0.0D0)

C Eigenvalues and eigenvectors now found. Now suppose that Theta =
C exp(iH), and Theta' = i.Theta.Omega. From the known expression for Omega
C we can compute dH/dt.

          CALL COFMAT(X,ELAM,S11,S12,S21,S22,N)

          ELOC = 0.0D0
          DO 20 J = 1,N
              DO 10 I = 1,N
                  ELOC = MAX(ELOC,ABS(S11(I,J)),ABS(S12(I,J)),
     &                   ABS(S21(I,J)),ABS(S22(I,J)))
   10         CONTINUE
   20     CONTINUE
          ELMAX = MAX(ELOC,ELMAX)

          CALL DGEMM('T','N',N,N,N,1.0D0,V,IV,S11,N,0.0D0,STORE,N)
          CALL DGEMM('N','N',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,COFF11,N)

          CALL DGEMM('T','N',N,N,N,1.0D0,V,IV,S12,N,0.0D0,STORE,N)
          CALL DGEMM('N','N',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,COFF12,N)

          CALL DGEMM('T','N',N,N,N,1.0D0,V,IV,S22,N,0.0D0,STORE,N)
          CALL DGEMM('N','N',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,COFF22,N)

          DO 40 P = 1,N
              DO 30 Q = 1,N
                  COFF21(Q,P) = COFF12(P,Q)
   30         CONTINUE
   40     CONTINUE

          DO 60 Q = 1,N
              DO 50 P = 1,Q
                  CALL FGH(E(P),E(Q),F,G,GT,HL)
                  STORE(P,Q) = 0.5D0*( F*COFF11(P,Q)+G*COFF12(P,Q)+
     &                                 GT*COFF21(P,Q)+HL*COFF22(P,Q) )
                  STORE(Q,P) = STORE(P,Q)
   50         CONTINUE
   60     CONTINUE

          CALL SMCOPY('G',N,N,STORE,N,COFF11,N)
          CALL DGEMM('N','N',N,N,N,1.0D0,V,IV,COFF11,N,0.0D0,STORE,N)
          CALL DGEMM('N','T',N,N,N,1.0D0,STORE,N,V,IV,0.0D0,HDOT,IHDOT)

C We have now computed the derivative required.

      END IF

      END

C ----------------------------------------------------------------------

      SUBROUTINE ANSYM(A,IA,B,IB,N,C,IC)
C     .. Scalar Arguments ..
      INTEGER IA,IB,IC,N
C     ..
C     .. Array Arguments ..

      DOUBLE PRECISION A(IA,N),B(IB,N),C(IC,N)
C     ..
C     .. Local Scalars ..
      INTEGER I,J,K
C     ..
      DO 20 J = 1,N
          DO 10 I = 1,N
              C(I,J) = 0.D0
   10     CONTINUE
   20 CONTINUE

      DO 30 J = 1,N
         DO 40 I = 1,J
            DO 50 K = 1,N
               C(I,J) = C(I,J) + A(I,K)*B(K,J)
50          CONTINUE
40       CONTINUE
30    CONTINUE

C      DO 50 K = 1,N
C          DO 40 J = 1,N
C              DO 30 I = 1,J
C                  C(I,J) = C(I,J) + A(I,K)*B(K,J)
C   30         CONTINUE
C   40     CONTINUE
C   50 CONTINUE

      DO 70 J = 1,N
          DO 60 I = 1,J - 1
              C(J,I) = C(I,J)
   60     CONTINUE
   70 CONTINUE

      END


C ----------------------------------------------------------------------


      SUBROUTINE FGH(X,Y,F,G,GT,H)
      DOUBLE PRECISION X,Y,F,G,GT,H
      DOUBLE PRECISION T,PHI,RATIO,DIFF,SUM,COSDIF,COSSUM,SINDIF,SINSUM
      EXTERNAL PHI
      
      DIFF = (X-Y)*0.5D0
      SUM = DIFF + Y
      T = DIFF**2

      SINDIF = SIN(DIFF)

      IF (T.LT.0.1D0) THEN
C          RATIO = 1.D0/PHI(-T)
           RATIO = PHI(-T)
      ELSE
C          RATIO = DIFF/(SINDIF+SIGN(1.D0,SINDIF)*1.0D-14)
          RATIO = (SINDIF+SIGN(1.D0,SINDIF)*1.0D-14)/DIFF
      END IF

C      RATIO = 2.D0*RATIO
      RATIO = 0.5D0*RATIO

      COSDIF = COS(DIFF)
      COSSUM = COS(SUM)
      SINSUM = SIN(SUM)      

C      F = RATIO*(COSDIF-COSSUM)
C      G = RATIO*(SINDIF+SINSUM)
C      GT = RATIO*(SINSUM-SINDIF)
C      H = RATIO*(COSDIF+COSSUM)

      F = (COSDIF-COSSUM)/RATIO
      G = (SINDIF+SINSUM)/RATIO
      GT = (SINSUM-SINDIF)/RATIO
      H = (COSDIF+COSSUM)/RATIO

      RETURN
      END


C -------------------------------------------------------------------


      DOUBLE PRECISION FUNCTION phi(v)
C     .. Parameters ..
      DOUBLE PRECISION a0,a1,a2,a3,a4,a5
      PARAMETER (a0=1.D0,a1=a0/6.D0,a2=a1/20.D0,a3=a2/42.D0,a4=a3/72.D0,
     &          a5=a4/110.D0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION v
C     ..
      phi = a0 + v* (a1+v* (a2+v* (a3+v* (a4+v*a5))))
      RETURN

      END


C -------------------------------------------------------------------


      DOUBLE PRECISION FUNCTION TRACE(A,IA,N)
C     .. Scalar Arguments ..
      INTEGER IA,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(IA,N)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION SUM
      INTEGER J
C     ..
      SUM = 0.D0
      DO 10 J = 1,N
          SUM = SUM + A(J,J)
   10 CONTINUE
      TRACE = SUM
      END

C -------------------------------------------------------------------

      SUBROUTINE MONIT(A,IA,N,IPRN)

C     .. Scalar Arguments ..
      INTEGER IA,IPRN,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(IA,N)
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC REAL
C     ..
      DO 10 I = 1,N
          WRITE (IPRN,FMT=*) (REAL(A(I,J)),J=1,N)
   10 CONTINUE
      END

C -----------------------------------------------------


      SUBROUTINE FM2AJF(AR,IAR,AI,IAI,N,RR,RI,INTGER,
     &                  CWORK,IFAIL)
C
C Interface to replacement for withdrawn NAG routine F02AJF.
C
      IMPLICIT NONE
      INTEGER IAR,IAI,N,INTGER(N),I,J,INDEX,IFAIL
      DOUBLE PRECISION AR(IAR,N),AI(IAI,N),RR(N),RI(N)
      COMPLEX*16 CWORK(N*(N+65)),Z(1,1),VL(1,1),VR(1,1)


      INDEX = 0
      DO 25 J = 1,N
         DO 35 I = 1,N
            INDEX = INDEX + 1
            CWORK(INDEX) = DCMPLX(AR(I,J),AI(I,J))
35       CONTINUE
25    CONTINUE
    
      CALL ZGEEV('N','N',N,CWORK(1),N,CWORK(N*N+1),VL,1,VR,1,
     &           CWORK(N*N+N+1),64*N,AR,IFAIL)
      IF (IFAIL.NE.0) THEN
        WRITE (6,FMT=*) 'Failure exit from ZGEEV, INFO=',IFAIL
        RETURN
      END IF
      IFAIL = 0

      DO 10 I = 1,N
         RR(I) = DBLE(CWORK(N*N+I))
         RI(I) = DIMAG(CWORK(N*N+I))
C         WRITE (6,FMT=*) 'Eigenvalue (LAP):',RR(I),RI(I)
10    CONTINUE


      RETURN
      END

      DOUBLE PRECISION FUNCTION XM2AJF(X)
      DOUBLE PRECISION X
      DOUBLE PRECISION DUMMY
C Replacement for the copyright-protected NAG routine by the
C name X02AJF. Provides machine epsilon.
 
      DUMMY = 0.125D0
20    IF (ABS(1.D0+DUMMY).GT.1.D0) THEN
        DUMMY = DUMMY/2.D0
        GOTO 20
      END IF
      XM2AJF = DUMMY

      RETURN
      END

      DOUBLE PRECISION FUNCTION XM1AAF(X)
      DOUBLE PRECISION X
      DOUBLE PRECISION DUMMY
C Replacement for the copyright-protected NAG routine by the
C name X01AAF. Provides the value of PI.
 
      XM1AAF = 4.D0*ATAN(1.D0)

      RETURN
      END


      SUBROUTINE FM2AXF(AR,IAR,AI,IAI,N,R,VR,IVR,VI,IVI,WK1,WK2,WK3,
     &                  CWORK,IFAIL)
      IMPLICIT NONE
C
C Replacement for NAG routine F02AXF
C
      INTEGER IAR,IAI,N,IVR,IVI,IFAIL,INDEX,I,J
      DOUBLE PRECISION AR(IAR,N),AI(IAI,N),R(N),VR(IVR,N),VI(IVI,N),
     & WK1(N),WK2(N),WK3(N)
      COMPLEX*16 CWORK(N*(N+65))
      EXTERNAL F02AXF,ZHEEV

      INDEX = 0
      DO 25 J = 1,N
         DO 35 I = 1,N
            INDEX = INDEX + 1
            IF (J.LE.I) CWORK(INDEX) = DCMPLX(AR(I,J),AI(I,J))
35       CONTINUE
25    CONTINUE

      CALL ZHEEV('V','L',N,CWORK(1),N,R,CWORK(N*N+N+1),
     &           64*N,VR,IFAIL)

      INDEX = 0
      DO 20 J = 1,N
C         WRITE (6,FMT=*) 'Eigenvalue (LAP):',R(J)
         DO 30 I = 1,N
            INDEX = INDEX + 1
            VR(I,J) = DBLE(CWORK(INDEX))
            VI(I,J) = DIMAG(CWORK(INDEX))
30       CONTINUE
20    CONTINUE

      RETURN
      END


      SUBROUTINE FM2ABF(A,IA,N,R,V,IV,E,IFAIL)
      INTEGER IA,N,IV,IFAIL
      DOUBLE PRECISION A(IA,N),R(N),V(IV,N),E(N)
      DOUBLE PRECISION WORK(5)

      DO 10 I = 1,N
         DO 20 J = I,N
            V(J,I) = A(J,I)
20       CONTINUE
10    CONTINUE

      IF (N.GT.2) THEN
          CALL DSYEV('V','L',N,V,IV,R,A,N*N,IFAIL)
      ELSE
          CALL DSYEV('V','L',N,V,IV,R,WORK,5,IFAIL)
      END IF

      RETURN
      END


      SUBROUTINE MM1CAF(RV, M1, M2, ORDER, IFAIL)
C Replacement for NAG routine M01CAF
      INTEGER M1,M2,IFAIL
      DOUBLE PRECISION RV(*)
      CHARACTER*1 ORDER,ORDLOC

      ORDLOC = ORDER
      IF (ORDLOC.EQ.'A'.OR.ORDLOC.EQ.'a') ORDLOC = 'I'

      CALL DLASRT(ORDLOC,M2-M1+1,RV,IFAIL)

      RETURN
      END


      SUBROUTINE SMCOPY(MATRIX,M,N,A,LDA,B,LDB)
C
C Replacement for NAG routine F06QFF
C
      CHARACTER*1 MATRIX
      INTEGER M,N,LDA,LDB
      DOUBLE PRECISION A(LDA,N), B(LDB,N)
      INTEGER            I, J
C
C  SMCOPY copies the M by N matrix A to the M by N matrix B.
C
C  If MATRIX = 'G' or 'g' then A and B are general matrices;
C
C  if MATRIX = 'U' or 'u' then A and B are upper triangular,  
C                         and only elements in the upper triangle
C                         of each (including the diagonal) are
C                         accessed;
C
C  if MATRIX = 'L' or 'l' then A and B are lower triangular,  
C                         and only elements in the lower triangle
C                         of each (including the diagonal) are
C                         accessed.

      IF((MATRIX.EQ.'G').OR.(MATRIX.EQ.'g'))THEN
         DO 20 J = 1, N
            DO 10 I = 1, M
               B(I,J) = A(I,J)
   10       CONTINUE
   20    CONTINUE
      ELSE IF((MATRIX.EQ.'U').OR.(MATRIX.EQ.'u'))THEN
         DO 40 J = 1, N
            DO 30 I = 1, MIN(M,J)
               B(I,J) = A(I,J)
   30       CONTINUE
   40    CONTINUE
      ELSE IF((MATRIX.EQ.'L').OR.(MATRIX.EQ.'l'))THEN
         DO 60 J = 1, MIN(M,N)
            DO 50 I = J, M
               B(I,J) = A(I,J)
   50       CONTINUE
   60    CONTINUE
      END IF
C
      RETURN
      END

      SUBROUTINE SMLOAD(MATRIX,M,N,CONST,DIAG,A,LDA)
C
C Replacement for NAG routine F06QHF
C
      CHARACTER*1 MATRIX
      DOUBLE PRECISION CONST,DIAG
      INTEGER LDA,M,N
      DOUBLE PRECISION A(LDA,N)
      INTEGER            I, J
C     
C
C  SMLOAD forms the m by n matrix A whose diagonal elements are DIAG
C         and whose off-diagonal elements are CONST.
C
C
C  If MATRIX = 'G' or 'g' then A is regarded as a general matrix;
C
C  if MATRIX = 'U' or 'u' then A is regarded as upper triangular, and only
C                         elements A(i,j) for which  i.le.j are referenced;
C
C  if MATRIX = 'U' or 'u' then A is regarded as lower triangular, and only
C                         elements A(i,j) for which  i.ge.j are referenced;
C
C

      IF( ( MATRIX.EQ.'G' ).OR.( MATRIX.EQ.'g' ) )THEN
         DO 20 J = 1, N
            DO 10 I = 1, M
               A( I, J ) = CONST
   10       CONTINUE
   20    CONTINUE
         IF( CONST.NE.DIAG )THEN
            DO 30 I = 1, MIN( M, N )
               A( I, I ) = DIAG
   30       CONTINUE
         END IF
      ELSE IF( ( MATRIX.EQ.'U' ).OR.( MATRIX.EQ.'u' ) )THEN
         DO 50 J = 1, N
            DO 40 I = 1, MIN( M, J )
               A( I, J ) = CONST
   40       CONTINUE
   50    CONTINUE
         IF( CONST.NE.DIAG )THEN
            DO 60 I = 1, MIN( M, N )
               A( I, I ) = DIAG
   60       CONTINUE
         END IF
      ELSE IF( ( MATRIX.EQ.'L' ).OR.( MATRIX.EQ.'l' ) )THEN
         DO 80 J = 1, MIN( M, N )
            DO 70 I = J, M
               A( I, J ) = CONST
   70       CONTINUE
   80    CONTINUE
         IF( CONST.NE.DIAG )THEN
            DO 90 I = 1, MIN( M, N )
               A( I, I ) = DIAG
   90       CONTINUE
         END IF
      END IF
C
      RETURN
      END


      INTEGER FUNCTION PM1ABF(IFO,IFLAG,SRNAME,NREC,REC)
C
C Replacement for NAG routine P01ABF
C
      INTEGER IFO,IFLAG,NREC
      CHARACTER*80 REC(2)
      CHARACTER*6 SRNAME
      INTEGER I

      IF (IFLAG.EQ.0) THEN
          PM1ABF = 0
          RETURN
      END IF

      IF (IFO.EQ.-1) THEN
C Silent exit
          PM1ABF = IFLAG
      ELSE IF (IFO.EQ.1) THEN
          WRITE (6,FMT=*) 'Failure in routine',SRNAME
          WRITE (6,FMT=*) 'IFAIL = ',IFLAG
          WRITE (6,FMT=*) (REC(I),I=1,NREC)
          PM1ABF = IFLAG
          RETURN
      ELSE
          WRITE (6,FMT=*) 'Failure in routine',SRNAME
          WRITE (6,FMT=*) 'IFAIL = ',IFLAG
          WRITE (6,FMT=*) (REC(I),I=1,NREC)
          PM1ABF = IFLAG
          STOP
      END IF

      END

      SUBROUTINE FM2BJF(N,A,IA,B,IB,EPS1,ALFR,ALFI, 
     &                  BETA,MATV,V,IV,ITER,
     &                  WORK,IFAIL) 
      INTEGER N,IA,IB,IV,ITER(N),IFAIL 
      DOUBLE PRECISION A(IA,N),B(IB,N),EPS1,ALFR(N), 
     &                 ALFI(N),BETA(N),V(IV,N),
     &                 WORK(8*N)
      CHARACTER*1 JOBVL,JOBVR
      LOGICAL MATV 
      DOUBLE PRECISION VL(1,1)
      
      JOBVL = 'N'
      IF (MATV) THEN
          JOBVR = 'V'
      ELSE
          JOBVR = 'N'
      END IF

      CALL DGGEV(JOBVL,JOBVR,N,A,IA,B,IB,ALFR,ALFI,
     &           BETA,VL,1,V,IV,WORK,8*N,IFAIL)

      RETURN
      END