module IVPRKP ! version 20120824
! This file ivprkp.f90 is part of the EPPEER-package for explicit
! parallel peer methods by B.A.Schmitt and R.Weiner available at
! www.mathematik.uni-marburg.de\~schmitt\peer
! It contains implementations of well-known Runge-Kutta methods of other
! authors from literature
!  -- variable definitions
  integer, save :: rstages, rorder, rmxsteps, rkfcn,rksteps,rkrej
  real(8), save :: ratol, rrtol, rhmin, rhmax
  ! Runge-Kutta tables, function value first stage
  real(8), dimension(:), allocatable, save :: rkb, rkc, rke
  real(8), dimension(:,:), allocatable, save :: rka
  logical :: rfsal

  contains
  !-- RKSETPARAM sets integration parameters
  subroutine rksetparam(rtol,atol)
    implicit none
    real(8), intent(in) :: rtol, atol
    rrtol = rtol
    ratol = atol
  end subroutine rksetparam

  !-- RKMETH performs Runge-Kutta integration
  subroutine rkmeth(n,fcn,t,tend,y,hs,jobrep,par)
    implicit none
    interface
      pure subroutine fcn(t,u,udot,par)
      real(8), intent(in) :: t
      real(8), dimension(:), intent(in) :: u
      real(8), dimension(:), intent(out) :: udot
      real(8), dimension(:), intent(in) :: par
      end subroutine fcn
    end interface
    integer, intent(in) :: n
    integer :: jobrep ! control and status report
    real(8) :: t,tend,hs
    real(8), dimension(n) :: y
    real(8), dimension(:), intent(in) :: par
    integer :: stg,ic,id  ! local variables
    real(8) tc, rer, rc, hn
    real(8), dimension(n) :: fv1, yc, ye
    real(8), dimension(n,2:rstages) :: ks
    logical :: o1step
!   --
    rksteps = 0;  rkrej = 0
    rhmax = tend-t
    o1step = (jobrep.ne.0) ! only one single step
    jobrep = 0
    if (.not. allocated(rka)) jobrep=-1  ! rksetcoeff not called
    call fcn(t,y,fv1,par); rkfcn = 1
    if (hs==0D0) & ! compute initial stepsize
      call rkinsteps(n,fcn,t,y,fv1,hs,par)
    do while (abs(t-tend)>1D-14 .and. (jobrep==0)) ! repeat until tend
      rksteps = rksteps+1
      if ((t-tend)*(t+hs-tend)<0D0) hs = tend-t ! allows tend<t0
      do stg = 2,rstages ! compute stages
        tc = t+hs*rkc(stg)
        yc = y+hs*rka(stg,1)*fv1
        do ic = 2, stg-1
          yc = yc+hs*rka(stg,ic)*ks(:,ic)
        end do
        call fcn(tc,yc,ks(:,stg),par); rkfcn = rkfcn+1
      end do
      ye = hs*rke(1)*fv1 ! compute error estimate:
      do stg =2, rstages
        ye = ye + hs*rke(stg)*ks(:,stg)
      end do
      rer = 0.D0  ! compute error quotient
      do id = 1,n
        rc = ratol+rrtol*abs(y(id))
        rer = rer +(ye(id)/rc)**2
      end do
      rer = sqrt(rer/n)
      hn = 0.9D0*exp(-log(1.D-8+rer)/rorder) ! new stepsize estimate
      hn = hs*max(0.2D0,min(hn,2.D0)) ! safety interval
      if (rer > 1.D0) then ! reject step
        rkrej = rkrej+1
        if (hn<rhmin) jobrep=1! report failure, stepsize too small
      else ! accept step and compute new solution
        t = t+hs
        y = y+hs*rkb(1)*fv1
        do stg = 2,rstages
         y = y+ hs*rkb(stg)*ks(:,stg)
        end do
        if (o1step) then ! cancel after first step
          jobrep = -2
        else
          if (rfsal) then
            fv1 = ks(:,rstages)
          else
            call fcn(t,y,fv1,par); rkfcn=rkfcn+1
          end if
        end if
      end if
      if (rksteps > rmxsteps) jobrep=2 ! failure, too many steps
      hs = hn
    end do ! integration steps
  end subroutine rkmeth

  !-- RKSETCOEFF creates coefficient matrices, sets defaults
  subroutine rksetcoeff(mnr,mname)
    implicit none
    integer :: mnr
    character(len=16), intent(out), optional :: mname
    rfsal = .false.
    select case(mnr) !
    case(2) ! RK2(1)
      if (present(mname)) mname = "RK1(2)"
      rstages = 2;  rorder = 2
      allocate(rka(rstages,rstages),rkb(rstages),rkc(rstages),rke(rstages))
      rka(2,1)=1.D0/2.D0
      rkc = (/ 0.D0, 1.D0/2.D0 /)
      rkb = (/ 0.D0, 1.D0 /)
      rke = (/-1.D0, 1.D0 /)
    case(3) ! RK3(2)
      if (present(mname)) mname = "RK2(3)"
      rstages = 3;  rorder = 3
      allocate(rka(rstages,rstages),rkb(rstages),rkc(rstages),rke(rstages))
      rka(2,1)=1.D0
      rka(3,1)=1.D0/4.D0; rka(3,2)=1.D0/4.D0;
      rkc = (/ 0.D0, 1.D0, 0.5D0 /)
      rkb = (/ 1.D0/6.D0, 1.D0/6.D0, 2.D0/3.D0 /)
      rke = (/-1.D0/3.D0,-1.D0/3.D0, 2.D0/3.D0 /)
    case(4) ! RKFB4(5)
      if (present(mname)) mname = "RKFB4(5)"
      rstages = 6; rorder = 4
      allocate(rka(rstages,rstages),rkb(rstages),rkc(rstages),rke(rstages))
      rka(2,1) = 1D0/4
      rka(3,1) = 3D0/32; rka(3,2)=9D0/32
      rka(4,1)=1932.D0/2197; rka(4,2)=-72.D2/2197; rka(4,3)=7296.D0/2197
      rka(5,1)=439.D0/216; rka(5,2)=-8.D0; rka(5,3)=368.D1/513;
      rka(5,4)=-845.D0/4104;
      rka(6,1)=-8.D0/27; rka(6,2)=2.D0; rka(6,3)=-3544.D0/2565;
      rka(6,4)=1859.D0/4104; rka(6,5)=-11.D0/40;
      rkc = (/ 0.D0, 1.D0/4, 3.D0/8, 12.D0/13, 1.D0, 1.D0/2 /)
      rkb = (/ 25D0/216,0D0,1408D0/2565,2197D0/4104,-1D0/5,0D0 /)
      rke = (/ 16D0/135,0D0,6656D0/12825,28561D0/56430,-9D0/50,2D0/55 /)
      rke = rkb - rke ! for difference of solutions
    case(5) ! DoPri5(4)
      if (present(mname)) mname = "DoPri5(4)"
      rstages = 7; rorder = 5
      rfsal = .true. ! first stage as last
      allocate(rka(rstages,rstages),rkb(rstages),rkc(rstages),rke(rstages))
      rka(2,1)=0.2D0
      rka(3,1)=3.D0/40.D0; rka(3,2)=9.D0/40.D0
      rka(4,1)=44.D0/45.D0; rka(4,2)=-56.D0/15.D0; rka(4,3)=32.D0/9.D0
      rka(5,1)=19372.D0/6561.D0; rka(5,2)=-25360.D0/2187.D0;
      rka(5,3)=64448.D0/6561.D0; rka(5,4)=-212.D0/729.D0
      rka(6,1)=9017.D0/3168.D0; rka(6,2)=-355.D0/33.D0
      rka(6,3)=46732.D0/5247.D0; rka(6,4)=49.D0/176.D0
      rka(6,5)=-5103.D0/18656.D0
      rka(7,1)=35.D0/384.D0; rka(7,2)=0.D0; rka(7,3)=500.D0/1113.D0
      rka(7,4)=125.D0/192.D0; rka(7,5)=-2187.D0/6784.D0
      rka(7,6)=11.D0/84.D0
      rkb = (/35.D0/384.D0,0.D0,500.D0/1113.D0,125.D0/192.D0,-2187.D0/6784.D0, &
        11.D0/84.D0, 0D0 /)
      rkc = (/0D0, 0.2D0, 0.3D0,0.8D0,8.D0/9.D0,1.D0,1.D0 /)
      rke = (/71.D0/57600.D0,0.D0,-71.D0/16695.D0,71.D0/1920.D0,-17253.D0/339200.D0, &
       22.D0/525.D0,-1.D0/40.D0 /)
    case default !
      if (present(mname)) mname = "error"
      mnr = -1
    end select
    ! set default parameters
    rmxsteps = 1000000
    rhmin = 1.D-10
    ratol = 1.D-5; rrtol=ratol
  end subroutine rksetcoeff

!  -- inisteps computes initial stepsize [dopri/Hairer-Wanner]
  subroutine rkinsteps(n,fcn,t,y,fv1,hs,par)
    implicit none
    interface
      pure subroutine fcn(t,u,udot,par)
      real(8), intent(in) :: t
      real(8), dimension(:), intent(in) :: u
      real(8), dimension(:), intent(out) :: udot
      real(8), dimension(:), intent(in) :: par
      end subroutine fcn
    end interface
    integer, intent(in) :: n
    real(8), dimension(:), intent(in) :: y, fv1
    real(8), dimension(:), intent(in) :: par
    real(8) :: t,hs, h,h1,rc,dnf,dny,der2,der12
    real(8), dimension(n) :: y2, f2
    integer :: id
    dnf = 0.D0; dny = 0.D0
    do id = 1,n
      rc = ratol+rrtol*abs(y(id))
      dnf = dnf +(fv1(id)/rc)**2
      dny = dny +(y(id)/rc)**2
    end do
    if (min(dnf,dny)<1D-10) then
      h = 1.D-6
    else
      h =1.D-2*sqrt(dny/dnf)
    end if
    h = min(h,rhmax)
    y2 = y+h*fv1 ! explicit Euler step
    call fcn(t+h,y2,f2,par)
    der2 = 0.D0
    do id = 1,n
      rc = ratol+rrtol*abs(y(id))
      der2 = der2 +((f2(id)-fv1(id))/rc)**2
    end do
    der2 = sqrt(der2)/h ! estimate for second derivative
    der12 = max(der2,sqrt(dnf))
    if (der12<=1.D-15) then
      h1 = max(1.D-6,h*1.D-3)
    else
      h1 = (1.D-2/der12)**(1.D0/rorder)
    end if
    hs = min(1D2*h,h1,rhmax)
  end subroutine rkinsteps

  !-- RKGETSTATS returns integration statistics
  subroutine rkgetstats(nsteps,nrej,nfcn)
    implicit none
    integer, intent(out) :: nsteps, nrej, nfcn
    nsteps = rksteps
    nrej = rkrej
    nfcn = rkfcn
  end subroutine rkgetstats

!  -- FREECOEFF releases all memory
  subroutine rkfreecoeff
    deallocate(rka,rkb,rkc,rke)
  end subroutine rkfreecoeff

end module IVPRKP