module integrate

  ! http://www.mymathlib.com/c_source/diffeq/other/bulirsch_stoer.c
  ! updated to f90 and vectors

  use typedef, only: &
       real64, int64

  use errors, only: &
       IERR_OK, &
       IERR_NONCONV, &
       IERR_DIVZERO, &
       IERR_TAB, &
       IERR_EPS

  use parameters, only: &
       verb_integrate

  implicit none

  ! integer(kind=int64), parameter :: &
  !      nnnsteps = 12
  ! integer(kind=int64), dimension(nnnsteps), parameter :: &
  !      nnsteps = (/ 2,4,6,8,12,16,24,32,48,64,96,128 /)
  integer(kind=int64), parameter :: &
       nnnsteps = 9
  integer(kind=int64), dimension(nnnsteps), parameter :: &
       nnsteps = (/ 2,4,6,8,10,12,14,16,18 /)

  real(kind=real64), dimension(nnnsteps), parameter :: &
       rnsteps = 1.d0 / nnsteps

  abstract interface
     function template(t, y) result(yp)
       import real64
       real(kind=real64), intent(in) :: &
            t
       real(kind=real64), intent(in), dimension(:) :: &
            y
       real(kind=real64), dimension(size(y, 1)) :: &
            yp
     end function template
  end interface


  abstract interface
     subroutine regularize(y)
       import real64
       real(kind=real64), intent(inout), dimension(:) :: &
            y
     end subroutine regularize
  end interface


  abstract interface
     subroutine extrapolate(fzero, tableau, x, f, ntab, n, ierr)
       import real64, int64
       integer(kind=int64), intent(in) :: &
            n, ntab
       real(kind=real64), intent(out), dimension(n) :: &
            fzero
       real(kind=real64), intent(inout), dimension(n, ntab+1) :: &
            tableau
       real(kind=real64), intent(in), dimension(ntab) :: &
            x
       real(kind=real64), intent(in), dimension(n) :: &
            f
       integer(kind=int64), intent(out) :: &
            ierr
     end subroutine extrapolate
  end interface

contains

  function gragg(f, t0, y0, f0, h0, nstep, rnstep, n) result(y)

    use typedef, only: &
         real64, int64

    implicit none

    procedure(template) :: &
         f

    real(kind=real64), intent(in) :: &
         t0, h0
    real(kind=real64), intent(in), dimension(n) :: &
         y0, f0
    integer(kind=int64), intent(in) :: &
         n, nstep
    real(kind=real64), intent(in) :: &
         rnstep
    real(kind=real64), dimension(n) :: &
         y

    real(kind=real64) :: &
         h, h2, t
    real(kind=real64), dimension(n) :: &
         y1, y2
    integer(kind=int64) :: &
         i

   h = h0 * rnstep
   h2 =  h + h

   t = t0
   y = y0
   y1 = y + h * f0

   do i = 2, nstep
      t = t + h
      y2 = y + h2 * f(t, y1)
      y = y1
      y1 = y2
   end do
   y = 0.5d0 * (y + y1 + h * f(t0 + h0, y1))

 end function gragg


 subroutine ExtRat(fzero, tableau, x, f, ntab, n, ierr)

   use typedef, only: &
        real64, int64

   implicit none

   real(kind=real64), intent(out), dimension(n) :: &
        fzero
   real(kind=real64), intent(inout), dimension(n, ntab+1) :: &
        tableau
   real(kind=real64), intent(in), dimension(ntab) :: &
        x
   real(kind=real64), intent(in), dimension(n) :: &
        f
   integer(kind=int64), intent(in) :: &
        n, ntab
   integer(kind=int64), intent(out) :: &
        ierr

   ! local variables

   real(kind=real64), dimension(n) :: &
        up, across, denominator, dum, t
   integer(kind=int64) :: &
        col

   if (ntab == 1) then
      fzero(:) = f(:)
      tableau(:,1) = f(:)
      ierr = IERR_OK
      return
   endif
   if (x(ntab) == 0.d0) then
      fzero(:) = f(:)
      ierr = IERR_TAB
      return
   end if

   across(:) = 0.d0
   up(:) = tableau(:,1)
   tableau(:,1) = f(:)

   do col = 2,ntab
      denominator(:) = tableau(:,col-1) - across(:)
      if (all(denominator(:)==0.d0)) then
         if (verb_integrate) then
            write(6, "('[extrat] ERROR-A: division by zero in col =',i5)") &
                 col
         endif
         ierr = IERR_DIVZERO
         return
      end if
      where (denominator(:) == 0.d0)
         dum(:) = 1.d0
      elsewhere
         dum(:) = 1.d0 - (tableau(:,col-1) - up(:)) / denominator(:)
      end where
      denominator(:) = (x(ntab+1-col) / x(ntab)) * dum(:) - 1.d0
      if (all(denominator(:)==0.d0)) then
         if (verb_integrate) then
            write(6, "('[extrat] ERROR-B: division by zero in col =', i5)") &
                 col
         endif
         ierr = IERR_DIVZERO
         return
      endif
      where (denominator(:)==0.d0)
         t(:) = tableau(:,col-1)
      elsewhere
         t(:) = tableau(:,col-1)+(tableau(:,col-1)-up(:))/denominator(:)
      end where
      across(:) = up(:)
      up(:) = tableau(:,col)
      tableau(:,col) = t(:)
   end do
   ierr = IERR_OK
   fzero(:) = t(:)
   return
 end subroutine ExtRat


 subroutine ExtPol(fzero, tableau, x, f, ntab, n, ierr)

   use typedef, only: &
        real64, int64

   implicit none

   real(kind=real64), intent(out), dimension(n) :: &
        fzero
   real(kind=real64), intent(inout), dimension(n, ntab+1) :: &
        tableau
   real(kind=real64), intent(in), dimension(ntab) :: &
        x
   real(kind=real64), intent(in), dimension(n) :: &
        f
   integer(kind=int64), intent(in) :: &
        n, ntab
   integer(kind=int64), intent(out) :: &
        ierr

   ! local variables

   real(kind=real64), dimension(n) :: &
        back_two_columns, old_aux, vertical_diff, &
        backslant_diff, forwardslant_diff, denominator
   integer(kind=int64) :: &
        i

   if (ntab==1) then
      tableau(:,1) = f(:)
      ierr = IERR_OK
      return
   endif

   if (x(ntab)==0.d0) then
      fzero(:) = f
      ierr = IERR_TAB
      return
   endif

   back_two_columns(:) = 0.d0
   old_aux(:) = tableau(:,1)
   tableau(:,1) = f(:)
   do i = 1, ntab-1
      vertical_diff(:) = tableau(:,i) - old_aux(:)
      backslant_diff(:) = tableau(:,i) - back_two_columns(:)
      forwardslant_diff(:) = backslant_diff(:) - vertical_diff(:)
      denominator(:) = (x(ntab-i)/x(ntab)) * forwardslant_diff(:) - backslant_diff(:)
      if (all(denominator(:)==0.d0)) then
         if (verb_integrate) then
            write(6, "('[extpol] ERROR-A: division by zero for i =',i5)") &
                 i
         endif
         ierr = IERR_DIVZERO
         return
      endif
      back_two_columns(:) = old_aux(:)
      old_aux(:) = tableau(:,i+1)
      where (denominator(:)==0.d0)
         tableau(:,i+1) = tableau(:,i)
      elsewhere
         tableau(:,i+1) = tableau(:,i)+vertical_diff(:)*backslant_diff(:)/denominator(:)
      end where
   end do
   fzero(:) = tableau(:,ntab)
   ierr = IERR_OK
   return
 end subroutine ExtPol


 subroutine bs(f, t0, y0, h0, eps, ryscale, iext, n, y, hnew, neval, ierr)

   use typedef, only: &
        int64, real64

   implicit none

   procedure(template) :: &
        f
   real(kind=real64), intent(in) :: &
        t0, h0, eps
   real(kind=real64), intent(in), dimension(n) :: &
        y0, ryscale
   integer(kind=int64), intent(in) :: &
        n, iext

   real(kind=real64), intent(out), dimension(n) :: &
        y
   integer(kind=int64), intent(out) :: &
        neval, ierr
   real(kind=real64), intent(out) :: &
        hnew

   ! local variables

   procedure(extrapolate), pointer :: &
        extrapol
   integer(kind=int64) :: &
        i, jerr
   real(kind=real64), dimension(nnnsteps) :: &
        step_size2
   real(kind=real64), dimension(n,nnnsteps+1) :: &
        tableau
   real(kind=real64), dimension(n) :: &
        est, old_est, f0
   real(kind=real64) :: &
        eps2

   neval = 0

   if (eps <= 0.d0) then
      ierr = IERR_EPS
      return
   endif
   eps2 = eps**2

   if (iext==1) then
      extrapol => ExtRat
   else
      extrapol => ExtPol
   endif

   f0 = f(t0, y0(:))
   est(:) = gragg(f, t0, y0(:), f0(:), h0, nnsteps(1), rnsteps(1), n)
   neval = neval + nnsteps(1) + 2

   step_size2(1) = (h0 * rnsteps(1))**2
   y(:) = est(:)

   i = 1
   call extrapol(y, tableau, step_size2, est, i, n, jerr)

   if (jerr < 0) then
      ierr = jerr
      return
   endif

   do i = 2, nnnsteps
      old_est(:) = y(:)
      est(:) = gragg(f, t0, y0(:), f0(:), h0, nnsteps(i), rnsteps(i), n)
      neval = neval + nnsteps(i) + 1
      step_size2(i) = (h0  * rnsteps(i))**2
      call extrapol(y, tableau, step_size2, est, i, n, jerr)
      if (jerr < 0) then
         ierr = jerr
         return
      endif

      if (sum(((y(:) - old_est(:)) * ryscale(:))**2) < eps2) then
         if (i > 2) then
            hnew = 8.d0 * h0 * rnsteps(i-1)
         else
            hnew = h0
         endif
         ierr = 0
         return
      end if
   end do
   ierr = IERR_NONCONV

 end subroutine bs

end module integrate