!*********************************************************************
!  n		dimension of system, leading dimension of sv,q,b
!  k	  	set k==0 for initialisation, 
!	 	on exit, number of completed calls/iterations.
!  kmax		max number of iterations, leading dim of h
!  ncgd		number of converged eigenvalues requested
!  dp           dot product: d = dp(n,a,b)
!  tol		tolerence of converged eigvals from H_{k-1} and H_k
!  sv		on entry, A.q_k ,
!	 	on exit, latest basis vector q_{k+1}.
!  h		Hessenburg matrix H
!  q		Krylov subspace Q
!  b		converged eigenvectors
!  wr2, wi2	converged eigenvalues.
!  ifail == -1  converged!!!
!            0  ok, continue
!            1  reached max iterations
!            2  invariant subspace
!          > 2  error from Lapack
!	 					A. P. Willis 2002
!*********************************************************************
      subroutine arnold(n,k,kmax,ncgd,dp,tol,sv,h,q,b,wr2,wi2,ifail)
      implicit none

! arguments
      integer		n,k,kmax,ncgd
      double precision  dp
      external          dp
      double precision  tol, sv(n), q(n,kmax), b(n,ncgd*2),  &
      			wr2(ncgd), wi2(ncgd), h(kmax,kmax)
      integer		ifail

! local variables
      integer		i,j,m,p
      double precision	norm,nm,small
      parameter		(small=1d-80)

! lapack stuff
      double precision	scale(kmax), tau(kmax),  &
     			wr(kmax), wi(kmax), work((kmax+2)*kmax),  &
     			h2(kmax,kmax), h3(kmax,kmax), vr(kmax,ncgd*2)
      logical		select(kmax)
      integer		ilo, ihi, ifailr(ncgd*2)


!---------------------------------------------------------------------
      ifail = 0

!     too many iterations
      if(k.ge.kmax) then
         ifail = 1
         return
      end if

!  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Arnoldi Process
!
!     copy soln vect to krylov subspace at k+1
      do i = 1, n
         q(i,k+1) = sv(i)
      end do

!     modified gram-schmidt
      do j = 1, k
!        get inner product --> k^th col of kxk matrix h
         h(j,k) = dp(n,q(1,k+1),q(1,j))
!        orthogonalise against previous k vectors
         do i = 1, n
            q(i,k+1) = q(i,k+1) - h(j,k)*q(i,j)
         end do
      end do
!     reorthogonolise
      do j = 1, k
         norm = dp(n,q(1,k+1),q(1,j))
         do i = 1, n
            q(i,k+1) = q(i,k+1) - norm*q(i,j)
         end do
         h(j,k) = h(j,k) + norm
      end do

!     normalise the k+1 basis vector
      norm = dp(n,q(1,k+1),q(1,k+1))
      norm = dsqrt(norm)
      if(norm.lt.small) then
         ifail = 2
         return
      end if
      do i = 1, n
         q(i,k+1) = q(i,k+1) / norm
         sv(i) = q(i,k+1)
      end do

!     for starting vector just copy to subspace and normalise,
!     zero the hessenburg matrix
      if(k.eq.0) then
         do j = 1, kmax
            do i = 1, kmax
               h(i,j) = 0d0
            end do
         end do
         goto 199
      end if

!     subdiagonal of hessenburg
      h(k+1,k) = norm
!  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Eigenvalues / Eigenvectors
!
!     need large enough matrix to get eigvecs
      if(k.lt.ncgd) goto 199

! get eigvals and eigvecs of kxk hessenberg matrix h by QR
      do j = 1, k
         do i = 1, kmax
            h2(i,j) = h(i,j)
         end do
      end do

!     balance the matrix
      call dgebal('B',k,h2,kmax,ilo,ihi,scale,ifail)      
      if(ifail.ne.0) then
         ifail = ifail + 1000
         return
      end if

!     reduce
      m = (kmax+2)*kmax
      call dgehrd(k,ilo,ihi,h2,kmax,tau,work,m,ifail)
      if(ifail.ne.0) then
         ifail = ifail + 2000
         return
      end if

!     QR
      do j = 1, k
         do i = 1, kmax
            h3(i,j) = h2(i,j)
         end do
      end do
      call dhseqr('E','N',k,ilo,ihi,h3,kmax,wr,wi,0d0,1,work,m,ifail)      
      if(ifail.ne.0) then
         ifail = ifail + 3000
         return
      end if
      
!     order the eigvals, largest magnitude
!     compare with eigvals from k-1
!     copy new values to old
      do i = 1, k
         select(i) = .false.
      end do
      p = 0
      do i = 1, ncgd
         norm = -1d0
         do j = 1, k
            if(.not.select(j)) then
               nm = wr(j)*wr(j) + wi(j)*wi(j)
               if(nm.gt.norm) then
                  norm = nm
                  m = j
               end if            
            end if
         end do
         select(m) = .true.
         if(dabs(wi(m)).gt.small) select(m+1)=.true.
         nm = dsqrt( (wr(m)-wr2(i))**2 + (wi(m)-wi2(i))**2 )
         if(nm .lt. tol) p = p+1
         wr2(i) = wr(m)
         wi2(i) = wi(m)
      end do

!     check if reqiured num converged...
      if(p.lt.ncgd) goto 199
      
!     Enough converged!! 
!     Get eigvecs of hessenburg by inverse iteration.
      do i = 1, k
         if(select(i) .and. dabs(wi(i)).gt.small) select(i+1)=.false.
      end do
      call dhsein('R','Q','N',select,k,h2,kmax,wr,wi,  &
                  0d0,1,vr,kmax,ncgd*2,m, work, 1,ifailr,ifail )
      if(ifail.ne.0) then
         ifail = ifail + 4000
         return
      end if
      
      call dgebak('B','R',k,ilo,ihi,scale,m,vr,kmax,ifail)
      if(ifail.ne.0) then
         ifail = ifail + 5000
         return
      end if

!     ? check correct number of cols
      p = 0
      do i = 1, k
         if(select(i)) then
            p = p+1
            if(dabs(wi(i)).gt.small) p = p+1
         endif
      end do
      if(p.ne.m) then
         ifail = m + p*100 + 6000
      end if

!     order the cols of vr
      m = 0
      do j = 1, k
         if(select(j)) then
            i = 0
            p = 0
 101        continue
            i = i+1
            if(wr2(i).ne.wr(j) .or. wi2(i).ne.wi(j)) then
               p = p+1
               if(dabs(wi2(i)).gt.small) p = p+1
               goto 101
            end if
            p = p+1
            m = m+1
            do i = 1, k
               work((p-1)*k+i) = vr(i,m)
            end do
            if(dabs(wi(j)).gt.small) then
               p = p+1
               m = m+1
               do i = 1, k
                  work((p-1)*k+i) = vr(i,m)
               end do            
            end if
         end if
      end do
      do p = 1, m
         do i = 1, k
            vr(i,p) = work((p-1)*k+i)
         end do
      end do

!     Get eigvecs, m cols
      do p = 1, m
         do i = 1, n
            b(i,p) = 0d0
            do j = 1, k
               b(i,p) = b(i,p) + q(i,j)*vr(j,p)
            end do
         end do
      end do

!     converged !!!
      ifail = -1
!  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 199  continue
      k = k+1

      end
!*********************************************************************
!  wr + i wi = e^{sigma T}
!  overwrite wr,wi with sigma.
!  do nothing to eigvecs.
!  if  w v = A v  with w,v complex, then equating re and imag parts
!     wr vr - wi vi = A vr  and  wi vr + wr vi = A vi
!
!*********************************************************************
      subroutine logeig(n,ncgd,b,wr,wi,t,ifail)
      implicit none

! arguments
      integer		n, ncgd
      double precision  b(n,ncgd*2), wr(ncgd), wi(ncgd), t
      integer		ifail

! local variables
      integer		i,j,p
      double precision	norm,small,re,im,th,pi
      parameter		(small=1d-80)
      
!---------------------------------------------------------------------
      ifail = 0
      pi = 4d0 * datan(1d0)

!     log of the eigvals
      do i = 1, ncgd
         re = wr(i)
         im = wi(i)
         if( dabs(re).lt.small ) then
            if( dabs(im).lt.small ) then
               ifail = i
               return
            else if( im.gt.0d0 ) then
               th =  0.5d0*pi
            else 
               th = -0.5d0*pi
            end if
         else if( re.gt.0d0 ) then
            th = datan(im/re)
         else
            if( im.gt.0d0 ) then
               th =  pi + datan(im/re)
            else
               th = -pi + datan(im/re)
            end if
         end if
         norm = dsqrt( re*re + im*im )
         wr(i) = dlog( norm ) / t
         wi(i) = th / t
      end do
      
!     eigenvectors
!      p = 1
!      do i = 1, ncgd
!         if(dabs(wi(i)).lt.small) then
!            p = p+1
!         else
!            do j = 1, n
!               re = b(j,p)
!               im = b(j,p+1)
!               b(j,p)   = wr(i)*re - wi(i)*im
!               b(j,p+1) = wr(i)*im + wi(i)*re
!            end do
!            p = p+2
!         end if
!      end do

      end
!*********************************************************************