module integrate ! http://www.mymathlib.com/c_source/diffeq/other/bulirsch_stoer.c ! updated to f90 and vectors use types, only: & real64, int64 use errors, only: & IERR_OK, & IERR_NONCONV, & IERR_DIVZERO, & IERR_TAB, & IERR_EPS 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 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 types, 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 types, 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 print*,'[extrat] ERROR-A: division by zero in col=', col 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 print*,'[extrat] ERROR-B: division by zero in col =', col 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 types, 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 print*,'[extpol] ERROR-A: division by zero for i= ', i 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 types, 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