      subroutine eptrk(n,f,t,u,te,rtol,atol,mf,
     $     xsteps,work,lwork,info)
c-----!----------------------------------------------------------------- 
c
c  EPTRK ... Explicit pseudo two-step RK methods for parallel computers 
c            with shared memory.
c
c  EPTRK solves the initial value problems for nonstiff systems of
c  first order ODEs
c
c      du
c      -- = f(t,u), u=(u(1),...,u(n)), f=(f(1),...,f(n))
c      dt
c
c  Version of August 25, 1997.
c
c  We have implemented a 5-order and a 8-order method (see "mf" below) 
c  with stepsize control by embedding. 
c
c  The code uses the BLAS and the LAPACK libraries.
c
c  The code contains compiler directives "loop_parallel" for automatic 
c  parallelization for Convex Exemplar machines. Compilers for other
c  machines will ignore these lines.
c 
c
c  Authors:
c 
c     N. H. Cong,
c     Faculty of Mathematics, Mechanics and Informatics
c     Hanoi University of Sciences
c     90 Nguyen Trai, Dong Da, Hanoi, Vietnam
c 
c     H. Podhaisky, R. Weiner,
c     Fachbereich Mathematik und Informatik, Universitaet Halle
c     06099 Halle, Germany
c     e-mail: podhaisky@mathematik.uni-halle.de
c             rw@mail.mathematik.uni-halle.de
c 
c-----!----------------------------------------------------------------- 
c
c  References:
c
c  [1]  N.H. Cong: 
c         Explicit pseudo two-step Runge-Kutta methods for parallel
c         computers. Submitted for publication.
c
c  [2]  N.H. Cong, H. Podhaisky, R. Weiner:
c         Numerical experiments with some explicit pseudo two-step
c         RK methods on a shared memory computer. Technical Report.
c         FB Mathematik und Informatik, Universitaet Halle
c     
c For retrieving the latest version visit URL: 
c  http://www.mathematik.uni-halle.de/institute/numerik/software 
c
c-----!----------------------------------------------------------------- 
c
c  Input parameters:
c
c     n   = Number of equations of the ODE system.
c
c     f   = Name (external) of user-supplied subroutine computing
c           the value f(t,u):
c         
c               subroutine f(n,t,u,uprime)
c               integer n
c               real*8 t,u(n),uprime(n)
c           
c               uprime(1)=...
c               ...
c
c      t   = Initial t-value. Changed on exit.
c
c      u   = Initial u-value. Changed on exit.
c
c     te   = Final t-value.
c
c     rtol = Relative tolerance.
c
c     atol = Absolute tolerance.
c
c     mf   = Method flag: 
c               mf = 5 .. use 5-order method EPTRK5
c               mf = 8 .. use 8-order method EPTRK8
c       
c    xsteps = Maximal number of allowed steps.
c     
c     work = Array of working space of length "lwork".
c
c    lwork = Length of "work". "lwork" must be at least 45 * n. 
c
c
c-----!----------------------------------------------------------------- 
c
c  Output parameters:
c
c     info  = Reports on success upon return:
c
c       0    computation successful
c       1    working array too small
c       2    stepsize too small
c       3    NaN occured
c       4    too many steps
c       5    PIRK method for the first step does not converges.
c       6    input is not consistent
c
c  
c    For  statistics the common block "ceptrk1" can be used:
c
c       integer nfeval,nstep,naccpt,nrejct
c       common  /ceptrk1/nfeval,nstep,naccpt,nrejct
c
c    with
c
c      nfeval = total Number of function evaluations.
c      nstep  = total Number of computed steps.
c      naccpt = total Number of accepted steps.
c      nrejct = total Number of rejected steps.
c
c
c-----!----------------------------------------------------------------- 
      implicit none

      integer n,mf,xsteps,info,lwork
      real*8 t,u(n),te,rtol,atol,work(lwork)
      external f

      integer pmstart
      integer pmfm,pmfu0,pmfun,pmun,pmyn,pmfyn,pmy
      integer pmscal,pmuu
      integer pmend


      logical lpirk

      integer sx
      parameter (sx=8)

c --- Set Pointers in array "work".
      pmstart=1

      pmfm=pmstart
      pmfu0=pmfm+n
      pmfun=pmfu0+n
      pmfyn=pmfun+2*n*sx
      pmun=pmfyn+n*sx
      pmyn=pmun+n*sx
      pmy=pmyn+n*sx
      pmscal=pmy+n
      pmuu=pmscal+n

      pmend=pmuu+n-1

      if (pmend.gt.lwork) then
        write (*,*) 'lwork.lt.',pmend
        info = 1
        return
      end if

      if (mf.ne.5.and.mf.ne.8) then
        write (*,*) 'Illegal value for parameter mf=',mf
        write (*,*) 'The method flag "mf" must be 5 or 8'
        info = 6
        return
       end if


c --- we use the pirk method for the first step only
      lpirk=.false.

      call eptrkc(n,f,t,u,te,rtol,atol,mf,xsteps,sx,lpirk,
     $ work(pmfm),work(pmfu0),work(pmfun),work(pmfyn),work(pmun),
     $ work(pmyn),work(pmy),work(pmscal),work(pmuu),info )
      return
      end 

c-----!----------------------------------------------------------------- 
c --- S U B R O U T I N E   E P T R K C 
c-----!----------------------------------------------------------------- 
c 
c
      subroutine eptrkc(n,f,t,u,tend,rtol,atol,mf,xsteps,
     $ sxx,lpirk,fm,fu0,fun,fyn,un,
     $ yn,y,scal,uu,info)

      implicit none
      integer n,info,sxx,mf,order,xsteps
      external f
      real*8 t,u(n),tend,rtol,atol,h,hold

      logical lpirk

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


      integer update,all
      data update, all /0,1/

      integer i,j,m
      real*8 alpha,err
      
      real*8 dnrm2,hnew,told

      real*8 fm(n),fu0(n),fun(n,sxx,2),un(n,sxx),yn(n,sxx)
      real*8 fyn(n,sxx),y(n),scal(n),uu(n)

      integer sx,sc
      integer iu,iy
      parameter (sx=8)
      real*8 c(sx),a(sx,sx),q(sx,sx),r(sx,sx),ark(sx,sx)
      real*8 g(sx),p(sx,sx),p1(sx,sx),b(sx),b1(sx),eye(sx,sx)
      real*8 rr(sx,sx),qinv(sx,sx)

      logical reached

c iembd =
c    0 ... use p(p-2)-order pair for embedding
c    1 ... use p(p-1)-order pair for embedding

	  integer iembd
	  common /ceptrk2/ iembd

	  iembd=0

      nfeval=0
      nstep=0
      naccpt=0
      nrejct=0
      info=0

      reached=.false.
      iu=1
      iy=2

      alpha=1d0
c --- compute all coefficients of the method
      call coeff(sc,all,mf,order,alpha,c,q,r,rr,a,
     $           ark,g,p,p1,b,b1,eye,qinv)

c --- for embedding
      call daxpy(sc,-1d0,b,1,b1,1)
c
c inital step size
c
      h=1d-6
      hold=h

c
c Begin Integration
c
 1000 continue
      if (t.ge.tend.or.reached) return

      if (nstep.ge.xsteps) then
        info=4
        return
      end if

      if (t+h*1.01d0.gt.tend) then
        h=tend-t
        reached=.true.
      end if

      nstep=nstep+1
      do i=1,n
        scal(i)=1d0/(atol+rtol*dabs(u(i)))
      end do
      if (naccpt.eq.0.or.lpirk) then
c
c PIRK-method for the first step
c
 90   continue
      m=0
      do j=1,sc
        call dcopy(n,u,1,yn(1,j),1)
      end do
c ---
 50   continue
      do j=1,sc
        call f(n,t+c(j)*h,yn(1,j),fun(1,j,iu))
        nfeval=nfeval+1
      end do
      do j=1,sc
        call dcopy(n,u,1,un(1,j),1)
        do i=1,sc
          call daxpy(n,h*ark(j,i),fun(1,i,iu),1,un(1,j),1)
        end do
      end do
      m=m+1
c --- PIRK-method does not converges. 
c     Possible remedy: a smaller inital stepsize.
      if (m.ge.20) then
	   h=0.5d0*h
	   if (h.gt.1d-10*(tend-t)) then 
	   	 goto 90
	   end if
       info=5
       return
      end if
      call daxpy(n*sc,-1d0,un(1,1),1,yn,1)

      err=0d0
	  do j=1,sc
      do i=1,n
        err=err+(yn(i,j)*scal(i))**2
      end do
	  end do
      err=dsqrt(err/dble(sc*n))

      call dcopy(n*sc,un(1,1),1,yn,1)

      if (err.gt.1d-8) goto 50
c ---
      else
c
c EPTRK-method
c
      alpha=h/hold
c --- recompute the RK-coefficients
      call coeff(sc,update,mf,order,alpha,c,q,r,rr,
     $   a,ark,g,p,p1,b,b1,eye,qinv)

c --- stage values
c$dir loop_parallel
      do j=1,sc
        call dcopy(n,u,1,un(1,j),1)
        do i=1,sc
          call daxpy(n,h*a(j,i),fun(1,i,iy),1,un(1,j),1)
        end do
      end do

c --- function evaluations
c$dir loop_parallel

      do i=1,sc
        call f(n,t+c(i)*h,un(1,i),fun(1,i,iu))
      end do
      nfeval=nfeval+sc

      end if
     
c
c new solution uu
c
      call dcopy(n,u,1,uu,1)
      do i=1,sc
        call daxpy(n,h*b(i),fun(1,i,iu),1,uu,1)
      end do
       
c
c embedded solution fu0
c
      call dcopy(n,0d0,0,fu0,1)
      do i=1,sc
        call daxpy(n,h*b1(i),fun(1,i,iu),1,fu0,1)
      end do

c
c step size control
c
      err=0d0
      do i=1,n
        err=err+(fu0(i)*scal(i))**2
      end do
      err=dsqrt(err/dble(n))

      if (err.ne.err) then
      write (*,9000) t,naccpt
 9000 format ('NaN occured at t=',f8.2,1x,'naccpt=',i5)
        info=3
        return
      end if

      hnew=(1d0/err)**(1d0/(order-1+iembd)) * 0.8d0

      if (hnew.lt.0.3d0) hnew=0.3d0
      if (hnew.gt.3d0) hnew=3d0
      hnew=h*hnew

      if (1d-1*hnew.le.dabs(t)*1d-14) then
      write (*,*) 'Stepsize too small'
        info=2
        return
      end if

      if (err.lt.1) then
c
c Step is accepted
c
      naccpt=naccpt+1
      told=t
      hold=h
      t=t+h
      call dcopy(n,uu,1,u,1)
      iu=3-iu
      iy=3-iy
      h=hnew
      goto 1000
      else
c
c Step is rejected
c
      nrejct=nrejct+1
      reached=.false.
      h=hnew
      goto 1000
      end if
        
      end

c-----!----------------------------------------------------------------- 
c --- S U B R O U T I N E   C O E F F
c-----!----------------------------------------------------------------- 
c
c The subroutine computes all RK-coefficients if "task" is equal 1.
c Otherwise the matrix A which depends on the stepsize ratio 
c alpha=h/hold is recomputed only. 
c 
c 
c
      subroutine coeff(sc,task,mf,order,alpha,c,q,r,rr,a,
     $                ark,g,p,p1,b,b1,eye,qinv)
      implicit none 
      
      integer sx,sc,order,mf
      parameter (sx=8)
      real*8 c(sx),a(sx,sx),q(sx,sx),r(sx,sx),ark(sx,sx)
      real*8 g(sx),p(sx,sx),p1(sx,sx),b(sx),b1(sx),eye(sx,sx)
      real*8 rr(sx,sx),qinv(sx,sx)
      real*8 alpha

      integer task,all
      data  all /1/

      integer i,j,info,iwork(sx)

	  integer iembd
	  common /ceptrk2/ iembd


      if (task.eq.all) then
      
      if (mf.eq.5) then
        sc=5
        order=5
        c(1)= 0.0885879595127039473955461d0
        c(2)= 0.409466864440734710864926d0
        c(3)= 0.787659461760847056025242d0
        c(4)=1.0d0
        c(5)=1.0d0+c(2)
       else 
        sc=8
        order=8
        c(1)= 0.0571041961145176821931212d0
        c(2)= 0.276843013638123827680046d0
        c(3)= 0.583590432368916820056698d0
        c(4)= 0.860240135656219447847913d0
        c(5)=1.0d0
        c(6)=1.0d0+c(2)
        c(7)=1.0d0+c(3)
        c(8)=1.0d0+c(4)
      end if

      do i=1,sc
        g(i)=1d0/dble(i)
      end do

      do i=1,sc
        r(i,1)=1d0
        q(i,1)=1d0
        p1(i,1)=c(i)**1*g(1)
        do j=2,sc
          p1(i,j)=c(i)**j*g(j)
          r(i,j)=c(i)**dble(j-1)
          q(i,j)=(c(i)-1)**dble(j-1)
        end do
       end do

c
c A
c
      call dcopy(sx*sx,0d0,0,qinv,1)
      call dcopy(sc,1d0,0,qinv(1,1),sx+1)
      call dgesv(sc,sc,q,sx,iwork,qinv,sx,info)

      end if

      do i=1,sc
        p(i,1)=p1(i,1)
        do j=2,sc
          p(i,j)=alpha**(j-1)*p1(i,j)
        end do
      end do

      call dgemm('n','n',sc,sc,sc,1d0,p,sx,qinv,sx,0d0,a,sx)

      if (task.eq.all) then

c
c b
c
      call dcopy(sx*sx,0d0,0,eye,1)
      call dcopy(sc,1d0,0,eye(1,1),sx+1)
      call dcopy(sx*sx,r,1,rr,1)

      call dgesv(sc,sc,r,sx,iwork,eye,sx,info)
      if (info.ne.0) then
        write (*,*) 'eptrk::coeff info.ne.0'
        stop
      end if
      call dgemv('t',sc,sc,1d0,eye,sx,g,1,0d0,b,1)
c
c ARK
c
      call dgemm('n','n',sc,sc,sc,1d0,p1,sx,eye,sx,0d0,ark,sx)
c
c for embedding: b1
c
      call dcopy (sx,0d0,0,b1,1)

      do i=1,sc-2+iembd
        r(i,1)=1d0
        do j=2,sc-2+iembd
          r(i,j)=c(i+3-1-iembd)**dble(j-1)
        end do
       end do

      call dcopy(sx*sx,0d0,0,eye,1)
      call dcopy(sc,1d0,0,eye(1,1),sx+1)
      call dgesv(sc-2+iembd,sc-2+iembd,
     $           r,sx,iwork,eye,sx,info)
      if (info.ne.0) then
        write (*,*) 'eptrk::coeff info.ne.0'
        stop
      end if
      call dgemv('t',sc-2+iembd,sc-2+iembd,
     $           1d0,eye,sx,g,1,0d0,b1(3-iembd),1)

      end if

      return 
      end   

