C **************************  IMPORTANT *****************************
C                                                                   *
C All software in this file is (C) Copyright Marco Marletta, 1992.  *
C All rights reserved.                                              *
C                                                                   *
C *******************************************************************

      SUBROUTINE sl12f(elam1,eps,a,b,k,n,tol,hamat,sysbc1,xmesh,mesh,
     &                 work,iwork,lwk,ilwk,nxtrap,ifail)
C This subroutine solves eigenvalue problems for linear Hamiltonian
C systems in which the sub-matrix S22 is positive definite. The
C problems concerned are reducible to the form of a first order matrix
C differential equation
C Theta' = i Theta Omega(x,lambda,Theta)
C where Theta is a unitary matrix which must satisfy prescribed
C boundary conditions at two endpoints x=a and x=b.
C The method use by SL12F is to approximate the coefficient matrices in
C the Hamiltonian system by piecewise constants and transform the problem
C to one in Sturm Liouville form using a transformation due to Reid (1973).
C
C Variables:
C elam1: double precision. On entry elam1 should be set to an approximation
C     to the eigenvalue to be found. The value of elam1 is not important in
C     determining the accuracy of the final result but it may affect the
C     run times. On exit, elam1 specifies the approximation to the eigenvalue
C     computed by SL13F.
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 lies
C     within [elam1-C*eps,elam1+C*eps], where by elam1 we mean the initial
C     approximation provided by the user, and C is a constant of modest size.
C     Again the value of eps is not critical. On exit, eps specifies a
C     (generally pessimistic) assessment of the error in the final
C     approximation elam1.
C a,b: double precision. On entry, a and b specify the finite regular
C     endpoints of the interval on which the problem is posed. Unchanged
C     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 of
C     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.
C     By a process of Wynn extrapolation, SL13F will endeavour to ensure
C     that the error in the final approximation elam1 is less than tol.
C     Unchanged on exit.
C xmesh: double precision array of dimension (0:mesh). Used as workspace.
C mesh: integer. On entry, the first dimension of xmesh as declared in the
C     calling (sub)program. The user should assign mesh a large enough value
C     to ensure that SL12F is able to store a sufficiently fine mesh to meet
C     the user's accuracy requirements. For most purposes a mesh of 2000
C     intervals is adequate (i.e. setting mesh = 2000 should be sufficient).
C     On exit, mesh specifies the number of mesh intervals actually used.
C work: double precision array of dimension iwork. Used as workspace.
C iwork: integer. The first dimension of work as declared in the calling
C     (sub)program, iwork >= n*(12+64*n). Unchanged on exit.
C lwk: integer array of dimension ilwk. Used as workspace.
C ilwk: integer. The first dimension of lwk as declared in the calling
C     (sub)program, ilwk >= n. Unchanged on exit.
C nxtrap: integer. On entry, nxtrap need not be set. On exit, nxtrap
C     specifies the number of Wynn extrapolations which were required
C     to achieve the user's specified accuracy.
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 HAMAT and SYSBC1:
C
C Hamat:
C
C     SUBROUTINE hamat(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 Sysbc1:
C
C     SUBROUTINE sysbc1(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 sysbc1 must specify the boundary condition at
C x=a, while if iend=1 then sysbc1 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 sysbc1. 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 ismax,ismall,nsampl
      PARAMETER (ismax=10,ismall=300,nsampl=20)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION a,b,elam1,eps,tol
      INTEGER ifail,ilwk,iwork,k,mesh,n,nxtrap
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION work(((24+2*nsampl)*n+12)*n),xmesh(0:mesh)
      INTEGER lwk(ilwk)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION c,el,epso,err,eu,h,toloc
      INTEGER i,ic,ifo,imatch,istor,j,kloc,kntmax,kp,newmsh,nmesh,nrec,
     &        nsmall
      CHARACTER srname*6
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION argsmp(nsampl),dsampl(nsampl),elam(ismax),
     &                 xsmall(0:ismall)
      INTEGER nj(ismax)
      CHARACTER rec(3)*80
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x02ajf
      INTEGER p01abf
      EXTERNAL x02ajf,p01abf
C     ..
C     .. External Subroutines ..
      EXTERNAL chosc2,extrap,hamesh,solvsy
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,max,sqrt,float
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat,sysbc1
C     ..
      srname = 'SL12F'
      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*(12+64*n).OR.
     &    ilwk.LT.n.OR.eps.LE.0.0D0) THEN
C We have an input error
          ifail = 1
          IF (a.GE.b) THEN
              WRITE (REC,FMT=9000) a,b
9900          FORMAT ('** A must be less than B, you have',/,
     &                '** A = ',g18.8,' B = ',g18.8)
              nrec = 2
              GOTO 70
          END IF
          IF (n.LT.2) THEN
              WRITE (REC,FMT=9920) n
9920          FORMAT ('** N must be at least 2, you have',/,
     &                '** N = ',i8)
              nrec = 2
              GOTO 70
          END IF
          IF (tol.LE.0.0D0) THEN
              WRITE (REC,FMT=9930) tol
9930          FORMAT ('** TOL must be strictly positive, you have',/,
     &                '** TOL = ',g18.8)
              nrec = 2
              GOTO 70
          END IF
          IF (iwork.LT.n*(12+64*n)) THEN
              WRITE (REC,FMT=9940)
     &        n*(12+64*n),iwork
9940          FORMAT ('** IWORK must be at least ',i8,' you have',/,
     &                '** IWORK = ',i8)
              nrec = 2
              GOTO 70
          END IF
          IF (ilwk.LT.n) THEN
              WRITE (REC,FMT=9950) n,ilwk
9950          FORMAT ('** ILWK must be at least',i8,' you have',/,
     &                '** ILWK = ',i8)
              nrec = 2
              GOTO 70
          END IF
          IF (eps.LE.0.0D0) THEN
              WRITE (REC,FMT=9955) eps
9955          FORMAT ('** EPS must be strictly positive, you have',/,
     &                '** EPS = ',g18.8)
              nrec = 2
              GOTO 70
          END IF

      END IF

      nxtrap = 0
      n2 = n*n

      iul = 1
      ivl = iul+n2
      iur = ivl+n2
      ivr = iur+n2
      iwk = ivr+n2
      isampl = iwk+20*n2+12*n

      epso = eps
      err = epso
      c = 0.50d0* (a+b)
      toloc = min(1.0D-2,sqrt(tol))
      ifail = 0
C Set up a mesh using a fairly ad-hoc method
      CALL hamesh(a,b,c,n,elam1,xmesh,mesh,nmesh,ic,toloc,work,hamat,
     &            ifail)

      IF (ifail.NE.0) THEN
C There are two cases:
C IFAIL = 10: mesh refined as much as possible and still cannot meet
C             accuracy requirement.
C IFAIL = 11: mesh routine ran out of space.
          IF (ifail.EQ.11) THEN
              WRITE (rec,FMT=9000)

9000          FORMAT ('** Meshing routine ran out of space,',/,
     &               '** increase MESH or increase TOL.')

              ifail = 5

          ELSE
              WRITE (rec,FMT=9010)

9010          FORMAT ('** Meshing routine making little progress,',/,
     &               '** stepsize too small. Increase TOL.')

              ifail = 4
          END IF

          GO TO 70

      END IF

C If there are too many mesh points then we simply thin the mesh out by
C the crudest method possible.

10    IF (nmesh.GT.ismall) THEN
          newmsh = nmesh/2
          ic = ic/2
          DO 20 j = 1,newmsh
              xmesh(j) = xmesh(2*j)
20        CONTINUE
          xmesh(newmsh) = xmesh(nmesh)
          nmesh = newmsh
          GO TO 10
      END IF

      imatch = ic

      nj(1) = 1
      nj(2) = 2
      nj(3) = 3
      nj(4) = 4
      nj(5) = 6
      nj(6) = 8
      nj(7) = 10
      nj(8) = 14
      nj(9) = 18
      nj(10) = 24

C      nj(1) = 1
C      nj(2) = 2
C      nj(3) = 4
C      nj(4) = 8
C      nj(5) = 16
C      nj(6) = 32
C      nj(7) = 64
C      nj(8) = 128
C      nj(9) = 256
C      nj(10) = 512

C Store the mesh in the small array:
      DO 30 j = 0,nmesh
          xsmall(j) = xmesh(j)
30    CONTINUE
      nsmall = nmesh

      istor = 1
C Solve the problem

40    ifail = 0
      toloc = sqrt(x02ajf(0.0D0))
      kntmax = 100
      eps = max(epso/100.D0,10.0D0*tol,ABS(err))

      CALL solvsy(xmesh,nmesh,ic,k,elam1,eps,el,eu,hamat,sysbc1,toloc,
     &            work,lwk,n,kntmax,ifail)
      IF (ifail.NE.0) THEN
C Some useful messages for the user:
          nrec = 2
          IF (ifail.EQ.6) THEN
              ifail = 2
              WRITE (rec,FMT=9020)

9020          FORMAT ('** Error finding Theta matrix from boundary',/,
     &               '** conditions supplied.')

          ELSE IF (ifail.EQ.7) THEN
              ifail = 3
              WRITE (rec,FMT=9030)

9030          FORMAT ('** Too many miss-distance evaluations.')

              nrec = 1

          ELSE IF (ifail.EQ.14) THEN
              ifail = 6
              WRITE (rec,FMT=9040)

9040          FORMAT ('** The matrix V+iU has become rank deficient',/,
     &               '** Try reducing TOL.')

          ELSE IF (ifail.EQ.16) THEN
              ifail = 8
              WRITE (rec,FMT=9050)

9050          FORMAT ('** The eigenvalues of a Theta matrix cannot',/,
     &               '** be found, possibly due to ill-conditioning of',
     &               /,'** V+iU. Try reducing TOL.')

              nrec = 3

          ELSE IF (ifail.EQ.17) THEN
              ifail = 9
              WRITE (rec,FMT=9060)

9060          FORMAT ('** An H-matrix has occured whose eigenvalues',/,
     &               '** cannot be computed accurately.')

          ELSE IF (ifail.EQ.18) THEN
              ifail = 10
              WRITE (rec,FMT=9070)

9070          FORMAT ('** A non-Hermitian H-matrix has occurred',/,
     &               '** possibly due to roundoff error.')

          ELSE IF (ifail.EQ.19) THEN
              ifail = 7
              WRITE (rec,FMT=9080)

9080          FORMAT ('** The matrix V+iU has become ill-conditioned')

              nrec = 1

          ELSE IF (ifail.EQ.20) THEN
              ifail = 11
              WRITE (rec,FMT=9090)

9090          FORMAT ('** The coefficient S22 cannot be diagonalised',/,
     &               '** with sufficient accuracy.')

          ELSE IF (ifail.EQ.21) THEN
              ifail = 12
              WRITE (rec,FMT=9100)

9100          FORMAT ('** The coefficient S22 is not positive definite.'
     &               )

              nrec = 1
          END IF

          GO TO 70

      END IF

      IF (istor.EQ.1) THEN
C Now get a better matchpoint by calling the new matchpoint routine.

      CALL chosc2(a,b,ic,work(iul),work(ivl),n,work(iur),work(ivr),
     &            n,el,eu,xmesh,nmesh,hamat,sysbc1,work(iwk),lwk,
     &            n,k,nsampl,work(isampl),dsampl,argsmp,
     &            ifail)
      imatch = ic
      END IF

C      WRITE (6,*) ' nmesh,ic:',nmesh,ic
C      WRITE (6,FMT=*) 'Eigenvalue found was',elam1
      elam(istor) = elam1
      istor = istor + 1
      err = epso
      IF (istor.GT.2) THEN
C We have some eigenvalues available and we can do some Richardson
C extrapolation:

          CALL extrap(elam,work,work(istor),nj,istor-1,elam1,err)
C          WRITE (6,FMT=*) 'Best approx so far:',elam1
C          WRITE (6,FMT=*) 'Current error:',err
      END IF

      IF (abs(err).LT.tol*max(1.0D0,abs(elam1)) .OR. istor.GT.ismax.OR.
     &    2*nmesh.GT.mesh) THEN
C         WRITE (6,FMT=*) 'Number of meshpoints used:',nmesh
C         WRITE (6,FMT=*) 'Number of extrapolations:',istor - 2
          nxtrap = istor - 2
          eps = err
          mesh = nmesh
          RETURN

      END IF
C Now add more meshpoints and see what happens
      nmesh = nsmall*nj(istor)
      ic = imatch*nj(istor)
      kp = 1
      DO 60 i = 1,nsmall
          h = (xsmall(i)-xsmall(i-1))/float(nj(istor))
          DO 50 kloc = 1,nj(istor)
              xmesh(kp) = xmesh(kp-1) + h
              kp = kp + 1
50        CONTINUE
60    CONTINUE

      GO TO 40

70    ifail = p01abf(ifo,ifail,srname,nrec,rec)
      RETURN

      END

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

      SUBROUTINE extrap(a,b,c,nj,n,best,err)

C     .. Scalar Arguments ..
      DOUBLE PRECISION best,err
      INTEGER n
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a(n),b(n),c(n,n)
      INTEGER nj(n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION fac1
      INTEGER j,k,p
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC float
C     ..
      DO 10 j = 1,n
          b(j) = a(j)
10    CONTINUE

C      DO 30 k = 1,n - 1
C          fac1 = 2.0D0**(2*k)
C          DO 20 j = n,k + 1,-1
C              b(j) = (b(j)-b(j-1)/fac1)/ (1.0D0-1.0D0/fac1)
C20        CONTINUE
C30    CONTINUE
C
      DO 30 p = 1,n
          DO 20 j = 1,n
              c(j,p) = (1.0D0/float(nj(j)))** (2*p)
20        CONTINUE
30    CONTINUE

      DO 70 k = 1,n - 1
          DO 40 j = n,k + 1,-1
              best = (b(j)- (c(j,1)/c(j-1,1))*b(j-1))/
     &               (1.0D0-c(j,1)/c(j-1,1))
              IF (j.EQ.n.AND.k.EQ.n-1) err = best-b(j)
              b(j) = best
40        CONTINUE
          DO 60 p = n - k,1,-1
              DO 50 j = n,k + 1,-1
                  c(j,p) = (c(j,p+1)- (c(j,1)/c(j-1,1))*c(j-1,p+1))/
     &                     (1.0D0-c(j,1)/c(j-1,1))
50            CONTINUE
60        CONTINUE
70    CONTINUE

      best = b(n)

      RETURN

      END

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

      SUBROUTINE hamesh(a,b,c,n,elam,xmesh,maxmsh,nmesh,imatch,tol,work,
     &                  hamat,ifail)
C     .. Scalar Arguments ..
      DOUBLE PRECISION a,b,c,elam,tol
      INTEGER ifail,imatch,maxmsh,n,nmesh
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION work(12*n*n),xmesh(0:maxmsh)
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat
C     ..
C     .. Local Scalars ..
      INTEGER ia1,ia2,ia3,ib1,ib2,ib3,ic1,ic2,ic3,id1,id2,id3
C     ..
C     .. External Subroutines ..
      EXTERNAL hamish
C     ..
      ia1 = 1
      ia2 = ia1 + n*n
      ia3 = ia2 + n*n
      ib1 = ia3 + n*n
      ib2 = ib1 + n*n
      ib3 = ib2 + n*n
      ic1 = ib3 + n*n
      ic2 = ic1 + n*n
      ic3 = ic2 + n*n
      id1 = ic3 + n*n
      id2 = id1 + n*n
      id3 = id2 + n*n
      CALL hamish(a,b,c,n,elam,xmesh,maxmsh,nmesh,imatch,tol,work(ia1),
     &            work(ia2),work(ia3),work(ib1),work(ib2),work(ib3),
     &            work(ic1),work(ic2),work(ic3),work(id1),work(id2),
     &            work(id3),hamat,ifail)
      RETURN

      END

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

      SUBROUTINE hamish(a,b,c,n,elam,xmesh,maxmsh,nmesh,imatch,tol,a1,
     &                  a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3,hamat,ifail)

C This subroutine forms an initial mesh for a Hamiltonian system using
C the Pruess equidistribution heuristic.
C
C Equidistribution heuristic: (h**3)*max(coefficient'') = const.
C
C     .. Parameters ..
      DOUBLE PRECISION half,two,one,safe,trd,five,ten,safe3,p008,sixth,
     &                 zero,tenth,twelf
      PARAMETER (half=5.d-1,two=2.d0,one=1.d0,safe=0.8D0,trd=one/3.d0,
     &          five=5.d0,ten=10.d0,safe3=7.29d-1,p008=8.d-3,
     &          sixth=one/6.d0,zero=0.d0,tenth=one/ten,twelf=half*sixth)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION a,b,c,elam,tol
      INTEGER ifail,imatch,maxmsh,n,nmesh
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a1(n,n),a2(n,n),a3(n,n),b1(n,n),b2(n,n),b3(n,n),
     &                 c1(n,n),c2(n,n),c3(n,n),d1(n,n),d2(n,n),d3(n,n),
     &                 xmesh(0:maxmsh)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION di,err,h,hmax,qmax,ratio,tolow,x1,x2,x3,xend,xo
      INTEGER i,ic,icnt,idi,imesh,irf,j,kntr,ii,jj
      LOGICAL phase1,repeat
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x02ajf
      EXTERNAL x02ajf
C     ..
C     .. External Subroutines ..
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,log,max,min,sign,int
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat
C     ..
      tolow = safe*tol
      irf = int(log(x02ajf(1.0D0))*0.50D0/log(safe))

C First do the left hand half of the interval:
      xo = c/ (10.0D0+abs(c))
      xend = a/ (10.0D0+abs(a))
      imesh = maxmsh
      x1 = xo
      xmesh(imesh) = x1
      imesh = imesh - 1
      h = abs(xend-xo)/20
      di = -1.0D0
      idi = -1
      hmax = h
      kntr = 0
      icnt = 0
      phase1 = .true.
      repeat = .false.

10    x2 = x1 + h*di
      IF (x2*di.GT.xend*di) x2 = xend
C Get the error measure for the interval [x1,x2]:
C      CALL hamat(10.D0*x1/ (1.0D0-abs(x1)),elam,a1,b1,c1,d1,n)
      CALL hamat(10.0D0*x1/ (1.0D0-abs(x1)),elam,c1,d1,a1,b1,n)
C      CALL hamat(10.D0*x2/ (1.0D0-abs(x2)),elam,a2,b2,c2,d2,n)
      CALL hamat(10.0D0*x2/ (1.0D0-abs(x2)),elam,c2,d2,a2,b2,n)
      x3 = 0.50D0* (x1+x2)
C      CALL hamat(10.D0*x3/(1.0D0-abs(x3)),elam,a3,b3,c3,d3,n)
      CALL hamat(10.0D0*x3/ (1.0D0-abs(x3)),elam,c3,d3,a3,b3,n)

      DO 7000 jj = 1,n
         DO 7010 ii = 1,n
            c1(ii,jj)=-c1(ii,jj)
            c2(ii,jj)=-c2(ii,jj)
            c3(ii,jj)=-c3(ii,jj)
7010     CONTINUE
7000  CONTINUE

C Form the three-point difference in the usual way:
      err = 0.0D0
      qmax = 1.0D0
      DO 30 j = 1,n
          DO 20 i = 1,n
              qmax = max(qmax,abs(a3(i,j)),abs(b3(i,j)),abs(c3(i,j)),
     &               abs(d3(i,j)))
              a3(i,j) = a1(i,j) + a2(i,j) - 2.0D0*a3(i,j)
              b3(i,j) = b1(i,j) + b2(i,j) - 2.0D0*b3(i,j)
              c3(i,j) = c1(i,j) + c2(i,j) - 2.0D0*c3(i,j)
              d3(i,j) = d1(i,j) + d2(i,j) - 2.0D0*d3(i,j)
              err = max(err,abs(a3(i,j)),abs(b3(i,j)),abs(c3(i,j)),
     &              abs(d3(i,j)))
20        CONTINUE
30    CONTINUE
      err = ((1.0D0-abs(x3))**2)*h*err/sqrt(max(one,qmax))

      DO 35 j = 1,n
         DO 25 i = 1,n
            c3(i,j) = 0.50D0*(c2(i,j)-c1(i,j))
            d3(i,j) = 0.50D0*(d2(i,j)-d1(i,j))
25       CONTINUE
35    CONTINUE

      CALL f06yaf('N','N',n,n,n,1.0D0,a3,n,c3,n,0.0D0,b1,n)
      CALL f06yaf('T','N',n,n,n,1.0D0,a3,n,c3,n,0.0D0,b2,n)
      CALL f06yaf('N','N',n,n,n,1.0D0,a3,n,d3,n,0.0D0,d1,n)
      CALL f06yaf('T','N',n,n,n,1.0D0,a3,n,d3,n,0.0D0,d2,n)
      CALL f06yaf('N','N',n,n,n,1.0D0,c3,n,d3,n,0.0D0,b3,n)

      qmax = 0.0D0
      DO 37 j = 1,n
         DO 27 i = 1,n
            qmax = max(qmax,abs(b1(i,j)),abs(b2(i,j)),abs(b3(i,j)),
     &                 abs(d1(i,j)),abs(d2(i,j)))
27       CONTINUE
37    CONTINUE

      err = max(err,h**2*qmax/sqrt(max(one,qmax)))

C Now decide what to do at the next step:

      IF (err.GT.tol) THEN
          repeat = .true.
          ratio = safe
          IF (tol.LT.err*safe3) ratio = (tol/err)**trd
          h = h*ratio
          icnt = icnt + 1
          IF (icnt.LT.irf) GO TO 10
          ifail = 10
          RETURN

      END IF

      IF (err.LT.tolow .AND. icnt.EQ.0) THEN
          ratio = five
          IF (err.GT.tolow*p008) ratio = (tolow/err)**trd
          h = h*ratio
          h = min(h,hmax)
C The next piece of code gets the meshing routine 'on-scale' at the first step.
          IF (phase1) THEN
              repeat = .true.
              kntr = kntr + 1
              IF (kntr.LT.5) GO TO 10
          END IF

      END IF
C General step has been completed successfully, prepare for next step.
      phase1 = .false.
      xmesh(imesh) = x2
      IF (x2*di.LT.xend*di) THEN
          imesh = imesh + idi
          IF (imesh.LT.0.0D0 .OR. imesh.GT.maxmsh) THEN
              ifail = 11
              RETURN

          END IF

          x1 = x2
          GO TO 10
C Look to see if we have done the whole of the interval yet:
      ELSE IF (di*sign(1.0D0,a-c).GT.0.0D0) THEN
C Shift what we have already got of the mesh into the first half of the array
C xmesh:
          ic = maxmsh - imesh
          imatch = ic
          IF (ic.GE.maxmsh) THEN
              ifail = 11
              RETURN

          END IF

          DO 40 j = 0,ic
              xmesh(j) = xmesh(imesh+j)
40        CONTINUE
          xo = c/ (10.0D0+abs(c))
          xend = b/ (10.0D0+abs(b))
          x1 = xo
          imesh = ic + 1
          di = 1.0D0
          idi = 1
          kntr = 0
          icnt = 0
          phase1 = .true.
          repeat = .false.
          GO TO 10

      END IF

      nmesh = imesh
C Convert mesh back to original coordinates
      DO 50 i = 0,nmesh
          xmesh(i) = 10.D0*xmesh(i)/ (1.0D0-abs(xmesh(i)))
50    CONTINUE
      ifail = 0
      RETURN

      END

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

      SUBROUTINE chosc2(a,b,ic,ul,vl,ndiml,ur,vr,ndimr,el,eu,xmesh,
     &           nmesh,cofmat,sysbc,work,lwk,n,k,nsampl,wsampl,dlsamp,
     &           argsmp,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 ..
      INTEGER iparam
      PARAMETER (iparam=1)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION a,b,el,eu
      INTEGER ic,ifail,k,n,ndiml,ndimr,nmesh,nsampl
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION argsmp(nsampl),dlsamp(nsampl),ul(ndiml,n),
     &                 ur(ndimr,n),vl(ndiml,n),vr(ndimr,n),
     &                 work(n*(18*n+12+ndiml+ndimr)),
     &                 wsampl(2*nsampl*n**2),xmesh(0:nmesh)
      INTEGER lwk(n)
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL cofmat,sysbc
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION alfaj,alfmax,alfmin,alfsum,argc,argl,argr,c,
     &                 diff,dla11f,elam,shift,twopi
      INTEGER i,imissi,imissr,iri,irr,isampi,isampr,ix,iy,j,jmatch,
     &        jsampl,jstep,n2
      LOGICAL ising,lnor
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION params(iparam)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x01aaf
      EXTERNAL x01aaf
C     ..
C     .. External Subroutines ..
      EXTERNAL f02ajf,f06qff,f06yaf,spthta,sumvec,thstp4
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,atan2,dble
C     ..
      ifail = 0
      shift = dble(k) + 5.0D-1
      twopi = 2.0D0*x01aaf(0.0D0)
      n2 = n*n
      ic = nmesh/2
      jstep = max(1,nmesh/nsampl)
      nsmplo = min(nsampl,nmesh)
C      WRITE (6,*) 'nmesh,nsmplo,jstep:',nmesh,nsmplo,jstep

      elam = el
      diff = 100.0D0
C Set up initial condition at x=a:
10    CALL sysbc(0,ising,a,ul,vl,ndiml,n,elam)
C Here we have used ul and vl as storage for y and pdy respectively.
C Get the phase-angles of the initial Theta-matrix:
      CALL f06qff('G',n,n,ul,ndiml,work(n+1),n)
      CALL f06qff('G',n,n,vl,ndiml,work(n2+n+1),n)

      lnor = .false.
      ifail = 0
      CALL spthta(work(n+1),work(n2+n+1),n,n,work(1),work(2*n2+n+1),
     &            work(3*n2+n+1),lwk,lnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 5
          RETURN

      END IF

      CALL sumvec(work,n,argl)
      argl = -argl
      alfsum = argl

C Do case of lower eigenvalue approximationm first:
      params(1) = elam
C Integrate from x=a to x=b 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:

      isampr = 1
      isampi = n*n + 1
      DO 30 jsampl = 1,nsmplo
          ifail = 0
          DO 20 i = (jsampl-1)*jstep,jsampl*jstep - 1
C              WRITE (7,*) 'Stepping from ',i,' to ',i+1
              CALL thstp4(xmesh(i),xmesh(i+1),ul,vl,ndiml,n,elam,argl,
     &                    alfsum,work,lwk,cofmat,ifail)
20        CONTINUE
          IF (ifail.NE.0) RETURN
C Store the outcome of the integration:
          CALL f06qff('G',n,n,ul,ndiml,wsampl(isampr),n)
          CALL f06qff('G',n,n,vl,ndiml,wsampl(isampi),n)
          argsmp(jsampl) = argl
          isampr = isampi + n*n
          isampi = isampr + n*n
30    CONTINUE

C Set up initial condition at x=b:
      CALL sysbc(1,ising,b,ur,vr,ndimr,n,elam)
C Here we have used ur and vr as storage for y and pdy respectively.
C Get the phase-angles of the initial Theta-matrix:

      CALL f06qff('G',n,n,ur,ndimr,work(n+1),n)
      CALL f06qff('G',n,n,vr,ndimr,work(n2+n+1),n)

      lnor = .true.
      ifail = 0
      CALL spthta(work(n+1),work(n2+n+1),n,n,work(1),work(2*n2+n+1),
     &            work(3*n2+n+1),lwk,lnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 5
          RETURN

      END IF

      CALL sumvec(work,n,argr)
      argr = -argr
      alfsum = 0.0D0
      DO 25 i = 1,n
          IF (work(i).NE.twopi) alfsum = alfsum + work(i)
25    CONTINUE
      alfsum = -alfsum

C Integrate from x=b to x=c:
      jmatch = nsmplo
      DO 40 i = nmesh,nsmplo*jstep + 1,-1
C          WRITE (7,*) 'Stepping from ',i,' to ',i-1
          CALL thstp4(xmesh(i),xmesh(i-1),ur,vr,ndimr,n,elam,argr,
     &                alfsum,work,lwk,cofmat,ifail)
40    CONTINUE
      DO 70 jsampl = nsmplo - 1,1,-1
          ifail = 0
          DO 50 i = (jsampl+1)*jstep,jsampl*jstep + 1,-1
C              WRITE (7,*) 'Stepping from ',i,' to ',i-1
              CALL thstp4(xmesh(i),xmesh(i-1),ur,vr,ndimr,n,elam,argr,
     &                    alfsum,work,lwk,cofmat,ifail)
50        CONTINUE
          IF (ifail.NE.0) RETURN
C
C Form the product of the sigma matrices at the centre. What we actually want
C is Conj(ThetaR)*ThetaL = SigmaR*Conj(SigmaL)
C      = (ur + i*vr)*(transp(ul)-i*transp(vl)).
C
          isampi = 1+n2+2*(jsampl-1)*n2
          isampr = isampi - n*n
          imissr = 1
          CALL f06yaf('N','T',n,n,n,1.0D0,ur,ndimr,wsampl(isampr),n,
     &                0.0D0,work(imissr),n)
          CALL f06yaf('N','T',n,n,n,1.0D0,vr,ndimr,wsampl(isampi),n,
     &                1.0D0,work(imissr),n)

          imissi = imissr + n2
          CALL f06yaf('N','T',n,n,n,1.0D0,vr,ndimr,wsampl(isampr),n,
     &                0.0D0,work(imissi),n)
          CALL f06yaf('N','T',n,n,n,-1.0D0,ur,ndimr,wsampl(isampi),n,
     &                1.0D0,work(imissi),n)

C Get eigenvalues using F02AJF:

          irr = imissi + n2
          iri = irr + n

          ifail = -1
          CALL f02ajf(work(imissr),n,work(imissi),n,n,work(irr),
     &                work(iri),lwk,ifail)
          IF (ifail.NE.0) THEN
              ifail = 4
              RETURN
          END IF
C Now get the phase angles which are required to form the miss distance
C function.

          alfaj = 0.0D0
          argc = 0.0D0
          alfmax = -100.0D0
          alfmin = 100.0D0
          DO 60 j = 1,n
              ix = irr + j - 1
              iy = iri + j - 1
              alfaj = atan2(work(iy),work(ix))
              IF (alfaj.LT.0.0D0) alfaj = alfaj + twopi
              IF (alfaj.GT.alfmax) alfmax = alfaj
              IF (alfaj.LT.alfmin) alfmin = alfaj
C              WRITE (6,FMT=*) 'Eigenvalue from F02AJF',work(ix),work(iy)
              argc = argc + alfaj
60        CONTINUE

          argl = argsmp(jsampl)

          dla11f = (-argl+argr-argc)/twopi + n - shift

          IF (dla11f.GT.0.0D0 .AND. dla11f.LT.0.6D0) dla11f = alfmin/
     &        twopi

          IF (dla11f.LT.0.0D0 .AND. dla11f.GT.
     &        (-0.6D0)) dla11f = alfmax/twopi - 1.0D0

          IF (elam.EQ.el) THEN
              dlsamp(jsampl) = dla11f

          ELSE
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:',xmesh(i-1),
C     &                         dla11f - dlsamp(jsampl)
              IF (abs(dla11f-dlsamp(jsampl)).LE.diff) THEN
                  jmatch = jsampl
                  diff = abs(dla11f-dlsamp(jsampl))
                  ic = i-1
              END IF

          END IF

70    CONTINUE

      IF (elam.EQ.el) THEN
          elam = eu
          GO TO 10

      END IF

C      WRITE (6,FMT=*) 'jmatch,nsmplo,a,b:',jmatch,nsmplo,a,b

      c = xmesh(ic)
C      ic = jmatch*jstep
C      WRITE (6,FMT=*) 'Match-point chosen (c,ic): ',c,ic

      RETURN

      END

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

      SUBROUTINE solvsy(xmesh,nmesh,ic,k,elam,eps,el,eu,hamat,sysbc,tol,
     &                  work,lwk,n,kntmax,ifail)
C This subroutine solves a regular matrix-vector Sturm-Liouville system,
C provided there are no bugs in it.

C First find a bracket for the eigenvalue:
C     .. Scalar Arguments ..
      DOUBLE PRECISION el,elam,eps,eu,tol
      INTEGER ic,ifail,k,kntmax,n,nmesh
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION work(24*n**2+12*n),xmesh(0:nmesh)
      INTEGER lwk(n)
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat,sysbc
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION dl,du,e1,e2,shift,toloc,twopi
      INTEGER iflag,ind,knt
      LOGICAL duas
C     ..
C     .. Local Arrays ..
      DOUBLE PRECISION c1(17)
C     ..
C     .. External Functions ..
      DOUBLE PRECISION dlhsys,x01aaf
      EXTERNAL dlhsys,x01aaf
C     ..
C     .. External Subroutines ..
      EXTERNAL c05azf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,max,min,dble
C     ..
      twopi = 2.0D0*x01aaf(0.0D0)
      shift = dble(k) + 5.0D-1
      el = elam - eps
      eu = elam + eps
      duas = .false.
      knt = 0
      ifail = 0
10    dl = dlhsys(xmesh,nmesh,ic,el,hamat,sysbc,work,lwk,n,ifail) -
     &     shift
      IF (ifail.NE.0) RETURN
C      WRITE (6,FMT=*) 'elam,dmiss:',el,dl
      IF (dl.GT.0.0D0) THEN
          eu = el
          du = dl
          duas = .true.
          eps = 2.0D0*eps
          el = el - eps
          knt = knt + 1
          IF (knt.GT.kntmax) THEN
              ifail = 7
              RETURN
          END IF
          GO TO 10

      END IF

      IF (duas) GO TO 30
      knt = 0
      du = dlhsys(xmesh,nmesh,ic,eu,hamat,sysbc,work,lwk,n,ifail) -
     &     shift
      IF (ifail.NE.0) RETURN
C      WRITE (6,FMT=*) 'elam,dmiss:',eu,du
20    IF (dl*du.GT.0.0D0) THEN
C We have not yet bracketed an eigenvalue
          knt = knt + 1
          IF (knt.GT.kntmax) THEN
              ifail = 7
              RETURN

          END IF

          eps = 2.0D0*eps
          IF (du.LT.0.0D0) THEN
              el = eu
              dl = du
              eu = eu + eps
              du = dlhsys(xmesh,nmesh,ic,eu,hamat,sysbc,work,lwk,n,
     &             ifail) - shift
              IF (ifail.NE.0) RETURN
C              WRITE (6,FMT=*) 'elam,dmiss:',eu,du

          ELSE
C I think that this next chunk of code can never be executed unless
C by error in the function values!
              eu = el
              du = dl
              el = el - eps
              dl = dlhsys(xmesh,nmesh,ic,el,hamat,sysbc,work,lwk,n,
     &             ifail) - shift
              IF (ifail.NE.0) RETURN
C              WRITE (6,FMT=*) 'elam,dmiss:',el,dl
          END IF

          GO TO 20

      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 = tol
      c1(1) = du
      iflag = 1
      e1 = el
      e2 = eu

40    CALL c05azf(e1,e2,dl,toloc,1,c1,ind,iflag)

      IF (ind.EQ.0) GO TO 50

      knt = knt + 1
      IF (knt.GT.kntmax) THEN
          ifail = 7
          RETURN

      END IF

      dl = dlhsys(xmesh,nmesh,ic,e1,hamat,sysbc,work,lwk,n,ifail) -
     &     shift
      IF (ifail.NE.0) RETURN
C      WRITE (6,FMT=*) 'elam,dmiss:',e1,dl

      IF (dl.GT.0.0D0) eu = min(eu,e1)

      IF (dl.LT.0.0D0) el = max(el,e1)

      IF (abs(eu-el).LT.toloc*max(1.0D0,abs(el))) GO TO 50
C Modify to a continuous miss distance function

      IF (abs(dl).LT.0.6D0) THEN
          IF (dl.LT.0.0D0) THEN
              dl = work(n) - twopi
          ELSE
              dl = work(1)
          END IF
C          WRITE (6,FMT=*) 'dmiss (modified):',dl
      END IF

      GO TO 40

50    elam = el
C      WRITE (6,FMT=*) 'el,eu,eu-el:',el,eu,eu - el
      RETURN

      END

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

      DOUBLE PRECISION FUNCTION dlhsys(xmesh,nmesh,ic,elam,hamat,sysbc,
     &                 work,lwk,n,ifail)
C This routine computes the Atkinson miss-distance function for a
C Hamiltonian system.
C     .. Parameters ..
      DOUBLE PRECISION zero,one,two
      PARAMETER (zero=0.0D0,one=1.0D0,two=2.0D0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION elam
      INTEGER ic,ifail,n,nmesh
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION work(n* (24*n+12)),xmesh(0:nmesh)
      INTEGER lwk(n)
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat,sysbc
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION alfsum,argc,argl,argr,twopi
      INTEGER i,ir,iu,iuo,iv,ivo,iwork,j
      LOGICAL altnor,ising,matv
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x01aaf
      EXTERNAL x01aaf
C     ..
C     .. External Subroutines ..
      EXTERNAL ethtex,f06qff,m01caf,spthta,sumvec,thstp4
C     ..
      ifail = 0
      twopi = two*x01aaf(one)
      iu = 1
      iv = n*n + iu
      iuo = n*n + iv
      ivo = n*n + iuo
      ir = n*n + ivo
      iwork = ir + n*n
C Step 1: set up the boundary condition at the left-hand end.
      CALL sysbc(0,ising,xmesh(0),work(iu),work(iv),n,n,elam)

C Step 2: Compute the eigenvalues of the Atkinson theta-matrix associated with
C the initial conditions and store in work(1,..,n)
      altnor = .false.
      CALL spthta(work(iu),work(iv),n,n,work(iwork),work(ir),
     &            work(iwork+n),lwk,altnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 6
          RETURN

      END IF

      CALL sumvec(work(iwork),n,argl)
      alfsum = argl
C Step 3: Integrate the system, following the eigenvalues as we go.
      DO 10 i = 1,ic
C Integrate from xmesh(i-1) to xmesh(i):
          CALL thstp4(xmesh(i-1),xmesh(i),work(iu),work(iv),n,n,elam,
     &                argl,alfsum,work(iwork),lwk,hamat,ifail)
          IF (ifail.NE.0) RETURN
10    CONTINUE

C Completed left-hand leg. Now turn to right hand leg.
C Step 4: set up the boundary condition at the right-hand end.
      CALL sysbc(1,ising,xmesh(nmesh),work(iuo),work(ivo),n,n,elam)

C Step 5: Compute the eigenvalues of the Atkinson theta-matrix associated with
C the initial conditions.
      altnor = .true.
      CALL spthta(work(iuo),work(ivo),n,n,work(iwork),work(ir),
     &            work(iwork+n),lwk,altnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 6
          RETURN

      END IF

      CALL sumvec(work(iwork),n,argr)
      alfsum = 0.0D0
      DO 20 i = iwork,iwork + n - 1
          IF (work(i).NE.twopi) alfsum = alfsum + work(i)
20    CONTINUE

C Step 6: Integrate the system, following the eigenvalues as we go.
      DO 30 i = nmesh,ic + 1,-1
C Integrate from xmesh(i) to xmesh(i-1):
          CALL thstp4(xmesh(i),xmesh(i-1),work(iuo),work(ivo),n,n,elam,
     &                argr,alfsum,work(iwork),lwk,hamat,ifail)
          IF (ifail.NE.0) RETURN
30    CONTINUE

C We now have the eigenvalues for both left and right-hand theta matrices.
C Use these to compute the miss-distance function.
      matv = .true.
      CALL ethtex(work(iu),work(iv),work(iuo),work(ivo),n,n,work(iwork),
     &            work(ir),work(iwork+n),n* (2*n+3),lwk,n,matv,ifail)
      IF (ifail.NE.0) RETURN
C Order the eigenvalues:
      CALL m01caf(work(iwork),1,n,'A',ifail)

C      WRITE (6,FMT=*) 'Eigenvalues of theta_{miss}:',
C     &  (work(j),j=iwork,iwork+n-1)
      CALL sumvec(work(iwork),n,argc)
C We may now compute the required miss-distance function.
      dlhsys = (argl-argr-argc)/twopi + n
C Also for convenience of the calling program let us copy the phase angles
C into the first n elements of the array work.
      CALL f06qff('G',n,1,work(iwork),n,work(1),n)
      RETURN

      END

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

      SUBROUTINE thstp4(xo,xend,u,v,ndim,n,elam,argdet,alfsum,work,lwk,
     &                  hamat,ifail)
C     .. Scalar Arguments ..
      DOUBLE PRECISION alfsum,argdet,elam,xend,xo
      INTEGER ifail,n,ndim
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION u(ndim,n),v(ndim,n),work(n* (18*n+12+ndim))
      INTEGER lwk(n)
C     ..
C     .. Local Scalars ..
      INTEGER ia,ia0,ia1,ia2,ib,ibet,ibs,ic,id,ievals,im,ir,iscl,ist1,
     &        ist2,iwork,n2
C     ..
C     .. External Subroutines ..
      EXTERNAL thstp3
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat
C     ..
      n2 = n*n
      ia = 1
      ib = ia + n2
      ibs = ib + n2
      ic = ibs + n2
      ia0 = ic + n2
      ia1 = ia0 + n2
      ia2 = ia1 + n2
      ist1 = ia2 + n2
      ist2 = ist1 + n2
      im = ist2 + n2
      ir = im + n2
      ibet = ir + ndim*n
      id = ibet + n
      iscl = id + n
      ievals = iscl + n
      iwork = ievals + n

      ifail = 0
      CALL thstp3(xo,xend,u,v,ndim,n,elam,argdet,alfsum,work(ia),
     &            work(ib),work(ibs),work(ic),work(ia0),work(ia1),
     &            work(ia2),work(ist1),work(ist2),work(im),work(ir),
     &            work(ibet),work(id),work(iscl),work(ievals),
     &            work(iwork),lwk,hamat,ifail)

      RETURN

      END

      SUBROUTINE thstp3(xo,xend,u,v,ndim,n,elam,argdet,alfsum,a,b,bsqrt,
     &                  c,a0,a1,a2,store1,store2,m,r,beta,d,sclrow,
     &                  evals,work,lwk,hamat,ifail)
C The purpose of this subroutine is to integrate a Hamiltonian
C system approximately over an interval [xo,xend] on which the
C coefficients are constant.

C     .. Parameters ..
      DOUBLE PRECISION zero,qtr,half,one,two
      PARAMETER (zero=0.0D0,qtr=2.5D-1,half=0.50D0,one=1.0D0,two=2.0D0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION alfsum,argdet,elam,xend,xo
      INTEGER ifail,n,ndim
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a(n,n),a0(n,n),a1(n,n),a2(n,n),b(n,n),beta(n),
     &                 bsqrt(n,n),c(n,n),d(n),evals(n),m(n,n),r(ndim,n),
     &                 sclrow(n),store1(n,n),store2(n,n),u(ndim,n),
     &                 v(ndim,n),work(8*n*(n+1))
      INTEGER lwk(n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION alfnew,alpha,arg0,argbet,betah,csdk,diff,difmax,
     &                 dk,dkssdk,dyold,fac,h,ht,psdk,scmax,scmin,sdk,
     &                 sgn,ssdk,t,twopi,yold
      INTEGER i,imax,j,k,kpjm1n
      LOGICAL altnor
C     ..
C     .. External Functions ..
      DOUBLE PRECISION armax,armin,phi,psi,scl,spexp,x01aaf
      INTEGER x02bkf
      EXTERNAL armax,armin,phi,psi,scl,spexp,x01aaf,x02bkf
C     ..
C     .. External Subroutines ..
      EXTERNAL acomp,f02abf,f04aaf,f06qff,f06yaf,m01caf,matexp,spthta,
     &         sqrmat,sumvec
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,atan2,cos,max,min,sign,sin,sqrt,float
C     ..
C     .. Subroutine Arguments ..
      EXTERNAL hamat
C     ..
      ifail = 0

C Evaluate the coefficients and compute any square roots of
C matrices which may be required.
C      CALL hamat(half* (xo+xend),elam,a,b,c,store1,n)

      CALL hamat(half*(xo+xend),elam,c,store1,a,b,n)

      DO 7020 j = 1,n
         DO 7030 i = 1,n
            c(i,j) = -c(i,j)
7030     CONTINUE
7020  CONTINUE

C Get the matrix sqrt(b):

      CALL sqrmat(b,n,bsqrt,n,n,work,ifail)
      IF (ifail.NE.0) THEN
          IF (ifail.EQ.1) THEN
              ifail = 20

          ELSE
              ifail = 21
          END IF

      END IF

C Now form the matrix A0 = inv(bsqrt).a.bsqrt using an appropriate NAG
C routine. Note that bsqrt has to be copied into a storage array as
C this is overwritten by f04aaf
      CALL f06qff('G',n,n,bsqrt,n,store1,n)
      ifail = 1
      CALL f04aaf(store1,n,a,n,n,n,store2,n,work,ifail)
      IF (ifail.NE.0) THEN
C          WRITE (7,*) 'Matrix bsqrt (1):'
C          CALL monit(bsqrt,n,n,7)
          ifail = 21
          RETURN

      END IF
C The matrix in store2 is inv(bsqrt).a, now post-multiply by bsqrt
      CALL f06yaf('N','N',n,n,n,one,store2,n,bsqrt,n,zero,store1,n)
C The matrix in store1 is now A0
      CALL f06qff('G',n,n,store1,n,a0,n)
C Now form the matrix A1 = 0.5*(A0+transpose(A0)) and the matrix
C A2 = 0.5*(A0-transpose(A0)).
C (This is a pretty crude piece of code but it will do until such time as
C  I can think of something better
      DO 20 j = 1,n
          DO 10 i = j,n
              a1(i,j) = half* (a0(i,j)+a0(j,i))
              a1(j,i) = a1(i,j)
              a2(i,j) = half* (a0(i,j)-a0(j,i))
              a2(j,i) = -a2(i,j)
10        CONTINUE
20    CONTINUE
C Form the matrix C0 = bsqrt.c.bsqrt and overwrite it on c:
      CALL f06yaf('N','N',n,n,n,one,c,n,bsqrt,n,zero,store1,n)
      CALL f06yaf('N','N',n,n,n,one,bsqrt,n,store1,n,zero,c,n)
C
C Now we must form the matrix M(x), which is given by
C M(x) = exp(A2(x-xmid)) where xmid = (xo+xend)/2. We will compute this
C matrix for the case where x = xend; the case x=xo will be the inverse of
C the case x=xend. The routine matexp computes the exponential for us.
C
      fac = half* (xend-xo)
      DO 40 j = 1,n
          DO 30 i = 1,n
              store1(i,j) = a2(i,j)*fac
              store2(i,j) = zero
30        CONTINUE
40    CONTINUE

      CALL matexp(store1,store2,n,n,2,m,work(1),n,work(n*n+1),ifail)
      IF (ifail.NE.0) THEN
          IF (ifail.EQ.1) THEN
              ifail = 17

          ELSE
              ifail = 18
          END IF

          RETURN

      END IF
C      WRITE (6,*) 'M matrix is:'
C      DO 17 i = 1,n
C         WRITE (6,*) (real(m(i,j)),j=1,n)
C17    CONTINUE

C The matrix m has now been computed. Now carry out the initial transformation
C from the matrices u and v to the new matrices u2 and v2 given by
C u2 = inv(M(xo))*inv(bsqrt)*u, v2 = inv(M(xo))*(bsqrt*v+a1*inv(bsqrt)*u)
C Note that inv(M(xo)) = M(xend) = m, which saves us some work.
      CALL f06qff('G',n,n,bsqrt,n,store1,n)
C Use the array A0 to contain inv(bsqrt)*u, since we don't need its current
C contents any more.
      ifail = 1
      CALL f04aaf(store1,n,u,ndim,n,n,a0,n,work,ifail)
      IF (ifail.NE.0) THEN
C          WRITE (7,*) 'Matrix bsqrt (2):'
C          CALL monit(bsqrt,n,n,7)
          ifail = 21
          RETURN
      END IF
C Now form u2 (overwrite it on u):
      CALL f06yaf('N','N',n,n,n,one,m,n,a0,n,zero,u,ndim)

C Now form v2 (overwrite it on v):
      CALL f06yaf('N','N',n,n,n,one,a1,n,a0,n,zero,store1,n)
      CALL f06yaf('N','N',n,n,n,one,bsqrt,n,v,ndim,one,store1,n)
      CALL f06yaf('N','N',n,n,n,one,m,n,store1,n,zero,v,ndim)

C Now find the phase angles of the matrix (v2+iu2)*inv(v2-iu2) using the
C routine spthta.
      ifail = 0
      CALL spthta(u,v,n,ndim,work(1),r,work(n+1),lwk,.false.,ifail)
      IF (ifail.NE.0) THEN
          ifail = 16
          RETURN

      END IF

      CALL sumvec(work(1),n,argbet)
C In perparation for the integration step we compute the matrix
C which Reid claims is given by
C C2 = C0 + A1*A2 -A2*A1
C but is actually given by
C C2 = C0 + A1*A1 (which we overwrite on C)
C
      CALL f06qff('G',n,n,a1,n,a2,n)
      CALL f06yaf('N','N',n,n,n,one,a1,n,a2,n,one,c,n)
C      CALL f06yaf('N','N',n,n,n,-one,a2,n,a1,n,one,c,n)
C Now we are ready to go. The next chunk of code we lift straight from the
C old coefficient approximation routine.

      twopi = two*x01aaf(one)
      h = xend - xo
      ht = xend/ (one+abs(xend)) - xo/ (one+abs(xo))
      sgn = one
C Use a NAG matrix eigenvalue routine to reduce c to diagonal
C form. Note: it might be a good idea to do some balancing first.

C      WRITE (7,*) 'xo,xend:',xo,xend
C      WRITE (7,*) 'Q matrix after Reid''s transformation:'
C      CALL monit(c,n,n,7)

      CALL f02abf(c,n,n,d,r,ndim,work,ifail)
      IF (ifail.NE.0) THEN
          ifail = 9
          RETURN
      END IF

C      WRITE (7,*) 'Eigenvalues of Q matrix:'
C      DO 17 jj = 1,n
C         WRITE (7,*) jj,d(jj)
C17    CONTINUE

C The columns of the array r are the eigenvectors of c. We may put these
C into c for safe storage: having diagonalised c we do not need it any more.
      CALL f06qff('G',n,n,r,ndim,c,n)
C Transform dependent variables:
      CALL f06qff('G',n,n,u,ndim,work(1),n)
      CALL f06yaf('T','N',n,n,n,one,r,ndim,work,n,zero,u,ndim)
      CALL f06qff('G',n,n,v,ndim,work(1),n)
      CALL f06yaf('T','N',n,n,n,one,r,ndim,work,n,zero,v,ndim)
C Finished transformation.

C Integrate from xo to xend:

      DO 60 k = 1,n
          dk = -d(k)
          sdk = sqrt(abs(dk))
          psdk = sdk
          t = dk* (h**2)
          IF (t.GE.0.1) THEN
              csdk = cos(sdk*h)
              ssdk = sin(sdk*h)/sdk
              dkssdk = -dk*ssdk
              sclrow(k) = zero
              evals(k) = two*scl(sdk*h,one/sdk)

          ELSE IF (abs(t).LT.0.1) THEN
              csdk = psi(-t)
              ssdk = h*phi(-t)
              dkssdk = -dk*ssdk
              sclrow(k) = zero
              evals(k) = two*atan2(ssdk,csdk)

          ELSE
C             IF (t.LE.-0.1) THEN
              sclrow(k) = abs(sdk*h)
              arg0 = spexp(-two*sclrow(k))
              csdk = half* (one+arg0)
              ssdk = half*sign(one,h)* (one-arg0)/sdk
              dkssdk = -dk*ssdk
              evals(k) = two*atan2(ssdk,csdk)
          END IF

          kpjm1n = k
          DO 50 j = 1,n
              yold = u(k,j)
              dyold = v(k,j)
              u(k,j) = csdk*yold + ssdk*dyold
              v(k,j) = dkssdk*yold + csdk*dyold
C Additional matrix needed for technical reasons:
              work(kpjm1n) = csdk*dyold - ssdk*yold
              kpjm1n = kpjm1n + n
50        CONTINUE
60    CONTINUE
C Finished integration.

      yold = zero
      DO 80 j = 1,n
          DO 70 i = 1,n
              yold = max(yold,abs(u(i,j)),abs(v(i,j)))
70        CONTINUE
80    CONTINUE

      scmax = armax(sclrow,n)
      scmin = armin(sclrow,n)

      fac = qtr*float(min(127,-x02bkf(one)))

C      IF (scmax-scmin.GT.fac*abs(ht).OR.yold.GT.toobig.OR.
C     &    .NOT.cheap) THEN

C Get argdet((C+iS)*inv(C-iS)):
      CALL sumvec(evals,n,arg0)

C Get eigenvalues of phi (stored in evals):
      altnor = .false.
      CALL f06qff('G',n,n,u,ndim,work(n**2+1),n)
      ifail = 1
      CALL spthta(work(n**2+1),work(1),n,n,evals,r,work(2* (n**2)+1),
     &            lwk,altnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 16
          RETURN

      END IF

C Form the sum of these phase angles:
      CALL sumvec(evals,n,betah)

C Now get the eigenvalues of theta2 correct including multiples of 2pi.
C This is a bit of a pain because we don't actually need these eigenvalues
C except to find a well-conditioned normalising matrix.
      CALL spthta(u,v,n,ndim,d,r,work,lwk,altnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 16
          RETURN

      END IF
C Now prepare for possible renormalisation of y and dy. First step is to
C select alpha such that cos(alpha/2)*y-sin(alpha/2)*dy is nowhere near
C singular.

C First order the eigenvalues using m01caf:
      CALL m01caf(d,1,n,'A',ifail)

      difmax = twopi + d(1) - d(n)
      imax = n
      DO 90 i = 1,n - 1
          diff = d(i+1) - d(i)
          IF (diff.GT.difmax) THEN
              imax = i
              difmax = diff
          END IF

90    CONTINUE

C Select phase-shift alpha for next step:
      IF (imax.LT.n) THEN
          alpha = half* (d(imax+1)+d(imax))

      ELSE
          alpha = half* (twopi+d(1)+d(n))
      END IF

C work has to be of dimension at least n*(2*n+1) here
      altnor = .true.
C Re-assign R to be the orthogonal diagonalising matrix:

      CALL f06qff('G',n,n,c,n,r,ndim)
      CALL acomp(u,v,ndim,r,ndim,sclrow,alpha,work,n,altnor,sgn,ifail)
      IF (ifail.NE.0) THEN
          ifail = 19
          RETURN

      END IF

C We have now post-multiplied u2 and v2 by a suitable normalising
C matrix. We can transform them back to u and v now. The formulae are
C u = bsqrt*M(xend)*u2 and v = inv(bsqrt)*(M(xend)*v2-a1*M(xend)*u2)
      CALL f06yaf('N','N',n,n,n,one,m,n,u,ndim,zero,store1,n)
      CALL f06yaf('N','N',n,n,n,one,bsqrt,n,store1,n,zero,u,ndim)

      CALL f06yaf('N','N',n,n,n,-one,a1,n,store1,ndim,zero,store2,n)
      CALL f06yaf('N','N',n,n,n,one,m,n,v,ndim,one,store2,n)
C Now we premultiply the contents of store2 by inv(bsqrt):
      CALL f06qff('G',n,n,bsqrt,n,store1,n)
      ifail = 1
      CALL f04aaf(store1,n,store2,n,n,n,v,ndim,work,ifail)
      IF (ifail.NE.0) THEN
C          WRITE (7,*) 'Matrix bsqrt (3):'
C          CALL monit(bsqrt,n,n,7)
          ifail = 21
          RETURN
      END IF

C Now get the eigenvalues of the theta matrix:
      altnor = .false.
      ifail = 0
      CALL spthta(u,v,n,ndim,d,r,work,lwk,altnor,ifail)
      IF (ifail.NE.0) THEN
          ifail = 16
          RETURN

      END IF

      CALL sumvec(d,n,alfnew)

      argdet = argdet - alfsum + alfnew + arg0 + argbet - betah
      alfsum = alfnew

      RETURN

      END

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

      SUBROUTINE matexp(ar,ai,ia,n,iswtch,expr,expi,iexp,work,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))
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION c,s
      INTEGER i,index,ir,ivi,ivis,ivr,ivrs,iwk1,iwk2,iwk3,j
C     ..
C     .. External Subroutines ..
      EXTERNAL f02axf,f06yaf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,cos,exp,sin
C     ..
      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
          ifail = -1
          CALL f02axf(ar,ia,ai,ia,n,work(ir),work(ivr),n,work(ivi),n,
     &                work(iwk1),work(iwk2),work(iwk3),ifail)
          IF (ifail.NE.0) THEN
              ifail = 1
              RETURN

          END IF
C Now get the exponential of the eigenvalues:
          index = ir
          DO 10 j = 1,n
              work(index) = exp(work(index))
              index = index + 1
10        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 30 j = 1,n
              DO 20 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
20            CONTINUE
30        CONTINUE

          CALL f06yaf('N','T',n,n,n,1.0D0,work(ivr),n,work(ivrs),n,
     &                0.0D0,expr,iexp)
          CALL f06yaf('N','T',n,n,n,1.0D0,work(ivi),n,work(ivis),n,
     &                1.0D0,expr,iexp)

          CALL f06yaf('N','T',n,n,n,1.0D0,work(ivi),n,work(ivrs),n,
     &                0.0D0,expi,iexp)
          CALL f06yaf('N','T',n,n,n,-1.0D0,work(ivr),n,work(ivis),n,
     &                1.0D0,expi,iexp)

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

      ELSE
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 50 j = 1,n
              DO 40 i = 1,n
                  ar(i,j) = -ar(i,j)
40            CONTINUE
              IF (abs(ar(j,j)).GT.1.0D-6) THEN
                  ifail = 2
                  RETURN

              END IF

              ar(j,j) = 0.0D0
50        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 f02axf(ai,ia,ar,ia,n,work(ir),work(ivr),n,work(ivi),n,
     &                work(iwk1),work(iwk2),work(iwk3),ifail)
          IF (ifail.NE.0) THEN
              ifail = 1
              RETURN

          END IF

C Now we can return the matrix A to its original form:
          DO 70 j = 1,n
              DO 60 i = 1,n
                  ar(i,j) = -ar(i,j)
60            CONTINUE
70        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(work(ir+j-1))
                  s = sin(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 f06yaf('N','T',n,n,n,1.0D0,work(ivr),n,work(ivrs),n,
     &                0.0D0,expr,iexp)
          CALL f06yaf('N','T',n,n,n,-1.0D0,work(ivi),n,work(ivis),n,
     &                1.0D0,expr,iexp)

          CALL f06yaf('N','T',n,n,n,1.0D0,work(ivr),n,work(ivis),n,
     &                0.0D0,expi,iexp)
          CALL f06yaf('N','T',n,n,n,1.0D0,work(ivi),n,work(ivrs),n,
     &                1.0D0,expi,iexp)

      END IF

      ifail = 0
      RETURN

      END

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

      SUBROUTINE sqrmat(a,ia,as,ias,n,work,ifail)
C The purpose of this routine is to compute the square root of positive
C definite Hermitian matrix A.
C     .. Scalar Arguments ..
      INTEGER ia,ias,ifail,n
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a(ia,n),as(ias,n),work(n* (2*n+1))
C     ..
C     .. Local Scalars ..
      INTEGER i
C     ..
C     .. External Subroutines ..
      EXTERNAL f02abf,rsmble
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC sqrt
C     ..
      ifail = 0
C First we obtain the Schur decomposition of the matrix:
      CALL f02abf(a,ia,n,work(1),work(2*n+1),n,work(n+1),ifail)
      IF (ifail.NE.0) THEN
          ifail = 1
          RETURN

      END IF
C Now we take the square roots of the eigenvalues
      DO 10 i = 1,n
          IF (work(i).LT.0.0D0) THEN
              ifail = 2
              RETURN

          END IF

          work(i) = sqrt(work(i))
10    CONTINUE
C Now reassemble the matrix sqrt(A):
      CALL rsmble(as,ias,n,work(1),work(2*n+1),n)
      RETURN

      END

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

      SUBROUTINE rsmble(as,is,n,d,r,ir)

C     .. Scalar Arguments ..
      INTEGER ir,is,n
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION as(is,n),d(n),r(ir,n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION fac1,fac2
      INTEGER i,j,k
C     ..
      DO 20 j = 1,n
          DO 10 i = 1,n
              as(i,j) = 0.0D0
10        CONTINUE
20    CONTINUE

      DO 50 k = 1,n
          fac1 = d(k)
          DO 40 j = 1,n
              fac2 = r(j,k)
              DO 30 i = 1,n
                  as(i,j) = as(i,j) + r(i,k)*fac1*fac2
30            CONTINUE
40        CONTINUE
50    CONTINUE

      RETURN

      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
      RETURN

      END

C ---------------------------------------------------------------------
C ---------------------------------------------------------------------
C
C Additional subroutines required by SL12F.
C (C) Copyright Marco Marletta, 1991. All rights reserved.
C

      SUBROUTINE spthta(y,dy,n,ndim,omega,r,work,lwk,lnor,ifail)
C Transfer y and dy values in order to avoid overwriting.
C     .. Parameters ..
      DOUBLE PRECISION zero,one
      PARAMETER (zero=0.d0,one=1.d0)
C     ..
C     .. Scalar Arguments ..
      INTEGER ifail,n,ndim
      LOGICAL lnor
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION dy(ndim,n),omega(n),r(ndim,n),work(n* (2*n+3)),
     &                 y(ndim,n)
      INTEGER lwk(n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION alpha,beta,fac,eps,twopi
      INTEGER i,knt,n2
      LOGICAL matv
C     ..
C     .. External Subroutines ..
      EXTERNAL f06qff,f02bjf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC atan2,sign
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x01aaf,x02ajf
      EXTERNAL x01aaf,x02ajf
C     ..
      LOGICAL prin
      COMMON/printb/prin
      n2 = n**2
      CALL f06qff('G',n,n,y,ndim,work(1),n)
      CALL f06qff('G',n,n,dy,ndim,work(n*n+1),n)
      ifail = 1
      matv = .false.
      eps = x02ajf(zero)
      twopi = 2.0D0*x01aaf(zero)
      knt = 0
100   CALL f02bjf(n,work(1),n,work(n*n+1),n,eps,work(2*n2+1),
     &            work(2*n2+n+1),work(2*(n2+n)+1),matv,r,ndim,lwk,
     &            ifail)
      IF (ifail.NE.0) THEN
          knt = knt + 1
          ifail = 1
          IF (knt.LT.4) GO TO 100
          ifail = 8
          RETURN
      END IF
      DO 30 i = 1,n
          alpha = work(2*n2+i)
          beta = work(2*(n2+n)+i)
          IF ((abs(alpha)+abs(beta)).EQ.zero) THEN
             ifail = 7
             RETURN
          ENDIF
          IF (alpha.EQ.zero) THEN
              omega(i) = 0.0D0
              IF (lnor) omega(i) = twopi
              GO TO 30
          END IF
          fac = sign(one,alpha)
          alpha = fac*alpha
          beta = fac*beta
          omega(i) = 2.d0*atan2(alpha,beta)
30    CONTINUE

      RETURN

      END

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

      SUBROUTINE ethtex(ul,vl,ur,vr,nsys,ndim,omega,vmat,wk,iwk,lwk,
     &                  ilwk,matv,ifail)

C This routine acts as the acceptable face of SUBROUTINE sptht0.
C Its purpose is described in the in-code comment for sptht0.

C iwk must be at least 2*(nsys**2) + 3*nsys.

C     .. Scalar Arguments ..
      INTEGER ifail,ilwk,iwk,ndim,nsys
      LOGICAL matv
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION omega(nsys),vmat(1:ndim,1:nsys),ul(ndim,ndim),
     &                 ur(ndim,ndim),vl(ndim,ndim),vr(ndim,ndim),wk(iwk)
      INTEGER lwk(ilwk)
C     ..
C     .. Local Scalars ..
      INTEGER i,i1,i2,i3,i4,i5,i6
C     ..
C     .. External Subroutines ..
      EXTERNAL sptht0
C     ..
      i1 = 1
      i = nsys**2
      i2 = i + 1
      i3 = i2 + i
      i4 = i3 + nsys
      i5 = i4 + nsys

      CALL sptht0(ul,vl,ur,vr,nsys,ndim,omega,wk(i1),wk(i2),vmat,
     &            wk(i3),wk(i4),wk(i5),lwk,matv,ifail)

      RETURN

      END

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

      SUBROUTINE sptht0(ul,vl,ur,vr,nsys,ndim,omega,a,b,v,alfar,alfai,
     &                  beta,iter,matv,ifail)
      DOUBLE PRECISION zero,one
      PARAMETER (zero=0.d0,one=1.d0)
C Subroutine to find the eigenvalues of the unitary matrix thetar^{*}.thetal
C where thetal and thetar are unitary matrices given by the Cayley transform
C          thetal(x) = (vl(x)+iul(x))(vl(x)-iul(x))^{-1}
C          thetar(x) = (vr(x)+iur(x))(vr(x)-iur(x))^{-1}
C
C The eigenvalues are all of the form cos(w)+isin(w), where w satisfies
C
C   det( cos(w/2)(vr^{*}ul-ur^{*}vl) - sin(w/2)(vr^{*}vl+ur^{*}ul) ) = 0
C We use the QZ algorithm to find cos(w/2) and sin(w/2), and thus recover w.
C
C     .. Scalar Arguments ..
      INTEGER ifail,ndim,nsys
      LOGICAL matv
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION a(nsys,nsys),alfai(nsys),alfar(nsys),
     &                 b(nsys,nsys),beta(nsys),omega(nsys),
     &                 ul(ndim,nsys),ur(ndim,nsys),v(ndim,nsys),
     &                 vl(ndim,nsys),vr(ndim,nsys)
      INTEGER iter(nsys)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION eps1,fac,sum,ulmax,urmax,vlmax,vrmax
      INTEGER i
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x02aaf,x01aaf
      EXTERNAL x02aaf,x01aaf
C     ..
C     .. External Subroutines ..
      EXTERNAL f06yaf,f06qff,f06qhf,f02bjf,m01caf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,atan2
C     ..
      eps1 = 0.d0
      ifail = 0
C First set A and B to zero:
      CALL f06qhf('G',nsys,nsys,zero,zero,a,nsys)
      CALL f06qhf('G',nsys,nsys,zero,zero,b,nsys)

      IF (matv) THEN
C Form the matrices A and B given by
C       A = vr^{*}ul - ur^{*}vl
C       B = vr^{*}vl + ur^{*}ul
C (where * denotes transpose)
        fac = 1.d0
        CALL f06yaf('T','N',nsys,nsys,nsys,fac,vr,ndim,ul,ndim,fac,
     &               a,nsys)
        CALL f06yaf('T','N',nsys,nsys,nsys,-fac,ur,ndim,vl,ndim,fac,
     &               a,nsys)
        CALL f06yaf('T','N',nsys,nsys,nsys,fac,vr,ndim,vl,ndim,fac,
     &               b,nsys)
        CALL f06yaf('T','N',nsys,nsys,nsys,fac,ur,ndim,ul,ndim,fac,
     &               b,nsys)
      ELSE

C Set a = ul and b = vl
       CALL f06qff('G',nsys,nsys,ul,ndim,a,nsys)
       CALL f06qff('G',nsys,nsys,vl,ndim,b,nsys)

      ENDIF

C Step 2: find the simultaneous eigenvalues of A and B.

      ifail = 1
      CALL f02bjf(nsys,a,nsys,b,nsys,eps1,alfar,alfai,beta,matv,v,ndim,
     &            iter,ifail)
      IF (ifail.NE.0) THEN
         IFAIL = 8
         RETURN
      END IF

      sum = 0.d0
      DO 13 i = 1,nsys
          fac = sign(one,alfar(i))
          omega(i) = 2.d0*atan2(fac*alfar(i),fac*beta(i))
          sum = sum + abs(alfai(i))/nsys
13    CONTINUE

      IF (sum.GT.x02aaf(0.d0)) THEN
          WRITE (6,FMT=*) 'Warning: imaginary eigenvalues detected'
          write(6,*) 'size of sum:',sum
      END IF

      RETURN

      END

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

      SUBROUTINE acomp(z,dz,iz,r,ir,sclrow,alpha,work,n,back,sgn,ifail)

C     .. Scalar Arguments ..
      DOUBLE PRECISION alpha,sgn
      INTEGER ifail,ir,iz,n,ndim
      LOGICAL back
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION dz(iz,n),r(ir,n),sclrow(n),work(8*n*n+n),z(iz,n)
C     ..
C     .. External Subroutines ..
      EXTERNAL compa
C     ..
      CALL compa(z,dz,iz,r,ir,sclrow,alpha,work(1),work(n*n+1),
     &           work(3*n*n+1),work(5*n*n+1),n,back,sgn,ifail)

      RETURN

      END

      SUBROUTINE compa(z,dz,iz,r,ir,sclrow,alpha,astore,bstore,cstore,
     &                 work,n,back,sgn,ifail)
C     .. Parameters ..
      DOUBLE PRECISION zero,half,one
      PARAMETER (zero=0.d0,half=5.d-1,one=1.d0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION alpha,sgn
      INTEGER ifail,ir,iz,n,ndim
      LOGICAL back
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION astore(n,n),bstore(2*n*n),
     &                 cstore(2*n*n),dz(iz,n),r(ir,n),work(n*(3*n+1)),
     &                 z(iz,n),sclrow(n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION ca2,dummy,d1,eps,sa2,arg1,arg2,arg3
      INTEGER i,id,j,k,p,ipjm1n,jpim1n,lwk(10)
      LOGICAL doy
C     ..
C     .. External Subroutines ..
      EXTERNAL f06qff,f06qhf,f03aff,f04ahf,getzdz,f04ayf,f01btf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC cos,sin
C     ..
C     .. External Functions ..
      DOUBLE PRECISION spexp,x01aaf,x02ajf
      EXTERNAL spexp,x01aaf,x02ajf
C     ..
      eps = x02ajf(one)
      sa2 = sin(half*alpha)
      ca2 = cos(half*alpha)

C Transpose z and dz:
      CALL transp(z,iz,n)
      CALL transp(dz,iz,n)
C Form matrix A:
      CALL lncomb(astore,n,z,iz,dz,iz,n,ca2,-sa2)
      CALL f06qff('G',n,n,astore,n,work(1),n)
C Now multiply z and dz by the inverse of a:
C Step 1: copy z and dz into a single unit of storage:

      doy = .true.
      IF (ca2**2.GE.half) doy = .false.

      IF (doy) THEN
         CALL f06qff('G',n,n,z,iz,bstore(1),n)
      ELSE
         CALL f06qff('G',n,n,dz,iz,bstore(1),n)
      END IF

C Step 2: solve an appropriate linear system with multiple RHS:
      ifail = 1
C Get PLU factorisation:
      CALL f01btf(n,work(1),n,work(n*n+1),d1,ifail)
      IF (ifail.NE.0) THEN
          IFAIL = 10
          RETURN
      END IF
C      CALL f03aff(n,eps,work(1),n,d1,id,work(n*n+1),ifail)
      sgn = sign(one,d1)
      index = 1
      DO 17 i = 1,n
         sgn = sgn*sign(one,work(index))
         index = index + n + 1
17    CONTINUE
C Solve system:
      ifail = 1
      CALL f06qff('G',n,n,bstore(1),n,cstore(1),n)
      CALL f04ayf(n,n,work(1),n,work(n*n+1),cstore(1),n,ifail)
C      CALL f04ahf(n,n,astore,n,work(1),n,work(n*n+1),bstore,n,eps,
C     &            cstore,n,work(n*(n+1)+1),n,k,ifail)
      IF (ifail.NE.0) THEN
          IFAIL = 10
          RETURN
      END IF
C Recover the required solutions:
      IF (doy) THEN
C Recover z:
         CALL f06qff('G',n,n,cstore(1),n,z,iz)
C Store an identity matrix:
         CALL f06qhf('G',n,n,zero,one,cstore(1),n)
C Get dz from z:
         CALL getzdz(z,dz,iz,cstore(1),n,n,ca2,sa2,1)
      ELSE
C Recover dz:
         CALL f06qff('G',n,n,cstore(1),n,dz,iz)
C Store an identity matrix:
         CALL f06qhf('G',n,n,zero,one,cstore(1),n)
C Get z from dz:
         CALL getzdz(z,dz,iz,cstore(1),n,n,ca2,sa2,0)
      END IF
C Take account to sclrow scaling:
      CALL recovr(astore,n,z,iz,sclrow,n,1)
      CALL recovr(bstore(1),n,dz,iz,sclrow,n,1)

C Transpose z and dz again:
      CALL transp(z,iz,n)
      CALL transp(dz,iz,n)

      CALL f06qff('G',n,n,astore,n,z,iz)
      CALL f06qff('G',n,n,bstore(1),n,dz,iz)

C We can set sclrow to zero now.
      CALL f06qhf('G',n,1,zero,zero,sclrow,n)

      IF (back) THEN
C Apply the transformation z = r.z, dz = r.dz
       IF (doy) THEN
C set astore = 0
          CALL f06qhf('G',n,n,zero,zero,astore,n)
          CALL f06yaf('N','N',n,n,n,one,r,ir,z,iz,one,astore,n)
          CALL f06qff('G',n,n,astore,n,z,iz)
          iopt = 1
       ELSE
          CALL f06qhf('G',n,n,zero,zero,astore,n)
          CALL f06yaf('N','N',n,n,n,one,r,ir,dz,iz,one,astore,n)
          CALL f06qff('G',n,n,astore,n,dz,iz)
          iopt = 0
       END IF

       CALL getzdz(z,dz,iz,r,ir,n,ca2,sa2,iopt)

      END IF

      RETURN

      END

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

      SUBROUTINE getzdz(z,dz,iz,r,ir,n,ca2,sa2,iopt)
      DOUBLE PRECISION zero,one
      PARAMETER (zero=0.d0,one=1.d0)
      INTEGER iz,n,iopt
      DOUBLE PRECISION z(iz,n),dz(iz,n),r(ir,n),ca2,sa2,arg1,arg2
      EXTERNAL lncomb
      IF (iopt.EQ.1) THEN
C Get dz from z
         arg1 = -one/sa2
         arg2 = ca2/sa2
         CALL f06qhf('G',n,n,zero,zero,dz,iz)
         CALL lncomb(dz,iz,r,ir,z,iz,n,arg1,arg2)
      ELSE
C Get z from dz
         arg1 = one/ca2
         arg2 = sa2/ca2
         CALL f06qhf('G',n,n,zero,zero,z,iz)
         CALL lncomb(z,iz,r,ir,dz,iz,n,arg1,arg2)
      END IF
      RETURN
      END

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

C Subroutine to compute the sum of the elements of a vector
      SUBROUTINE sumvec(vec,nvec,sum)
C     .. Scalar Arguments ..
      DOUBLE PRECISION sum
      INTEGER nvec
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION vec(nvec)
C     ..
C     .. Local Scalars ..
      INTEGER i
C     ..
      sum = 0.d0
      DO 10 i = 1,nvec
          sum = sum + vec(i)
10    CONTINUE
      RETURN

      END

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

      DOUBLE PRECISION FUNCTION armax(v,n)
C     .. Scalar Arguments ..
      INTEGER n
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION v(1:n)
C     ..
C     .. Local Scalars ..
      INTEGER i
C     ..
      armax = v(1)
      DO 10 i = 2,n
          IF (armax.LT.v(i)) armax = v(i)
10    CONTINUE
      RETURN

      END

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

      DOUBLE PRECISION FUNCTION armin(v,n)
C     .. Scalar Arguments ..
      INTEGER n
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION v(1:n)
C     ..
C     .. Local Scalars ..
      INTEGER i
C     ..
      armin = v(1)
      DO 10 i = 2,n
          IF (armin.GT.v(i)) armin = v(i)
10    CONTINUE
      RETURN

      END

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

      SUBROUTINE lindep(v1,v2,n,tf)
C This routine examines two vectors v1 and v2 for linear dependence
C     .. Parameters ..
      DOUBLE PRECISION zero,one
      PARAMETER (zero=0.d0,one=1.d0)
C     ..
C     .. Scalar Arguments ..
      INTEGER n
      LOGICAL tf
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION v1(n),v2(n)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION v1norm,v1v2ip,v2norm
      INTEGER i
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x02ajf
      EXTERNAL x02ajf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC abs,sqrt
C     ..
      v1norm = zero
      v2norm = zero
      v1v2ip = zero
      DO 10 i = 1,n
          v1norm = v1norm + v1(i)**2
          v2norm = v2norm + v2(i)**2
          v1v2ip = v1v2ip + v1(i)*v2(i)
10    CONTINUE
      IF (v1norm.LE.x02ajf(one) .OR. v2norm.LE.x02ajf(one)) THEN
          tf = .true.
          RETURN

      END IF

      IF (abs(abs(v1v2ip)/sqrt(v1norm*v2norm)-one).LT.
     &    sqrt(x02ajf(one))) THEN
          tf = .true.

      ELSE
          tf = .false.
      END IF

      RETURN

      END

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

      SUBROUTINE recovr(y,iy,a,ia,d,n,iopt)
      DOUBLE PRECISION zero
      PARAMETER (zero=0.d0)
      INTEGER iy,ia,n,iopt
      DOUBLE PRECISION y(iy,n),a(ia,n),d(n)
      DOUBLE PRECISION diff
      DOUBLE PRECISION spexp
      EXTERNAL spexp
      GOTO (100,200,300,400) iopt
100   DO 10 j = 1,n
         DO 20 i = 1,n
            diff = d(i)-d(j)
            IF (diff.LE.zero) THEN
               y(i,j) = spexp(diff)*a(j,i)
               y(j,i) = y(i,j)
            ELSE
               y(j,i) = spexp(-diff)*a(i,j)
               y(i,j) = y(j,i)
            END IF
20       CONTINUE
10    CONTINUE

      RETURN

200   DO 30 j = 1,n
         DO 40 i = 1,n
            y(i,j) = a(j,i)
40       CONTINUE
30    CONTINUE

      RETURN

300   DO 50 j = 1,n
         DO 60 i = 1,n
            diff = -d(i)-d(j)
            y(i,j) = spexp(diff)*a(j,i)
60       CONTINUE
50    CONTINUE

      RETURN

400   DO 70 j = 1,n
            diff = spexp(-d(j))
         DO 80 i = 1,n
            y(i,j) = diff*a(j,i)
80       CONTINUE
70    CONTINUE

      RETURN

      END

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

      SUBROUTINE transp(a,ia,n)
      INTEGER ia,n
      DOUBLE PRECISION a(ia,n)
      DOUBLE PRECISION dummy
      DO 10 i = 1,n
         DO 20 j = i+1,n
            dummy = a(i,j)
            a(i,j) = a(j,i)
            a(j,i) = dummy
20       CONTINUE
10    CONTINUE
      RETURN
      END

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

      SUBROUTINE lncomb(a,ia,b,ib,c,ic,n,facb,facc)
      INTEGER ia,ib,ic,n
      DOUBLE PRECISION a(ia,n),b(ib,n),c(ic,n)
      DOUBLE PRECISION facb,facc
      DO 10 i = 1,n
         DO 20 j = 1,n
            a(j,i) = facb*b(j,i)+facc*c(j,i)
20       CONTINUE
10    CONTINUE
      RETURN
      END

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

      DOUBLE PRECISION FUNCTION spexp(x)
C     .. Scalar Arguments ..
      DOUBLE PRECISION x
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION eneg
      LOGICAL first
C     ..
C     .. External Functions ..
      DOUBLE PRECISION x02amf
      EXTERNAL x02amf
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC exp,log
C     ..
C     .. Save statement ..
      SAVE first,eneg
C     ..
C     .. Data statements ..
      DATA first/.true./
C     ..
      IF (first) THEN
          eneg = log(x02amf(x))
          first = .false.
      END IF

      spexp = x02amf(x)
      IF (x.GT.eneg) spexp = exp(x)
      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 psi(v)
C     .. Parameters ..
      DOUBLE PRECISION b0,b1,b2,b3,b4,b5
      PARAMETER (b0=1.D0,b1=b0/2.D0,b2=b1/12.D0,b3=b2/30.D0,
     &           b4=b3/56.D0,b5=b4/90.D0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION v
C     ..
      psi = b0 + v* (b1+v* (b2+v* (b3+v* (b4+v*b5))))
      RETURN

      END

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

      DOUBLE PRECISION FUNCTION scl(th,u)
C     .. Parameters ..
      DOUBLE PRECISION one
      PARAMETER (one=1.d0)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION th,u
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC atan2,cos,sin
C     ..
      scl = th + atan2((u-one)*sin(th)*cos(th),one+ (u-one)*sin(th)**2)
      RETURN

      END