c
c Sample driver for EPTRK. 
c 
c link dr_eptrk eptrk diffu2 lapack blas
c
c Version of August 25, 1997.
c
c-----!----------------------------------------------------------------- 
      program dr_eptrk
      implicit none
      integer n,lwork,nx
      parameter (nx=2500,lwork=45*nx)
      integer info,tol,i
      integer mf
      integer xsteps

      integer nfeval,nstep,naccpt,nrejct
      common  /ceptrk1/nfeval,nstep,naccpt,nrejct

      real*8 u(nx),t,te,rtol,atol,work(lwork)
      real*8 u2(nx),errgem,errlog

      external f
      character*8 sample,method

c --- maximal number of steps
      xsteps=5000

c --- first use the 5-order method, next the 8-order method
      do mf=5,8,3

      write (method,9100) mf
 9100 format ('eptrk',i1)

c --- use tolerances 10^{-2} ... 10^{-6}
      do tol=-2,-6,-1
      atol=dexp(dble(tol)*dlog(1d1))
      rtol=atol
      
c --- get the initial values
      call problem(n,t,te,u,sample)

      if (n.gt.nx) then
        write (*,*) 'dr_eptrk: nx .lt. ',n
        stop
      end if

c
c Call to the integrator EPTRK.
c

      call eptrk(n,f,t,u,te,rtol,atol,mf,
     $  xsteps, work,lwork, info)


c 
c Compare with the reference solution.
c
      call loesung (n,te,u2)
      call daxpy (n,-1d0,u2,1,u,1)
      errgem=0d0
      do i=1,n
        errgem=errgem+(u(i)/(1+dabs(u2(i))))**2 
      end do
      errgem=dsqrt(errgem/n)
      errlog=dlog10(errgem)

      call summary(tol,naccpt,nrejct,nfeval,errlog,
     $     sample,method,info)

      end do
      end do
      end

c-----!----------------------------------------------------------------- 
c --- S U B R O U T I N E  S U M M A R Y
c-----!----------------------------------------------------------------- 
c This subroutine shows the integration statistics. If EPTRK
c gives a return value "info" not equal zero we put a indicating 
c character in the first column.
c
      subroutine summary(tol,naccpt,nrejct,nfeval,err,
     $          sample,method,info)
      implicit none
      integer tol,naccpt,nrejct,nfeval,info
      real*8 err
      character*8 sample,method
      character  code(0:6)
      data code /' ','I','S','N','X','C','?'/
c
c       0    computation successful
c       1    working array too small (input not valid)
c       2    stepsize too small
c       3    NaN occured
c       4    to many steps
c       5    convergence
c       6    input not consistent

   
      if(tol.eq.-2) then 
        write (*,2001) method,sample 
        write (*,2000)
 2000 format ( 
     $'  TOL  STEPS      FEVAL   ERRLOG ' )
 2001 format (/ 'Method ',a10,'on ',a10)
      end if

      write(*,1000)code(info),tol,naccpt,nrejct,nfeval,err

 1000 format (a1,1x,i2,i6,'+',i5,i7,1x,f6.2,f8.2)

      return
      end
