faccur Subroutine

private subroutine faccur(me, h0, h1, facc, x0, twoinf, f0, f1, ifail)

This procedure attempts to estimate the level of rounding errors in the calculated function values near the point x0+h0 by fitting a least-squares straight-line approximation to the function at the six points x0+h0-j*h1, (j = 0,1,3,5,7,9), and then setting facc to twice the largest deviation of the function values from this line. hi is adjusted if necessary so that it is approximately 8 times the smallest spacing at which the function values are unequal near x0+h0.

Type Bound

diff_func

Arguments

Type IntentOptional Attributes Name
class(diff_func), intent(inout) :: me
real(kind=wp), intent(inout) :: h0
real(kind=wp), intent(inout) :: h1
real(kind=wp), intent(out) :: facc
real(kind=wp), intent(in) :: x0
real(kind=wp), intent(in) :: twoinf
real(kind=wp), intent(in) :: f0
real(kind=wp), intent(in) :: f1
integer, intent(out) :: ifail

0 if no error, -1 if user termination.

Called by

proc~~faccur~~CalledByGraph proc~faccur diff_module::diff_func%faccur proc~diff diff_module::diff_func%diff proc~diff->proc~faccur proc~compute_jacobian_with_diff numerical_differentiation_module::compute_jacobian_with_diff proc~compute_jacobian_with_diff->proc~diff

Source Code

    subroutine faccur(me,h0,h1,facc,x0,twoinf,f0,f1,ifail)

    implicit none

    class(diff_func),intent(inout) :: me
    real(wp), intent(inout)  :: h0
    real(wp), intent(inout)  :: h1
    real(wp), intent(out)    :: facc
    real(wp), intent(in)     :: x0
    real(wp), intent(in)     :: twoinf
    real(wp), intent(in)     :: f0
    real(wp), intent(in)     :: f1
    integer, intent(out)     :: ifail  !! 0 if no error, -1 if user termination.

    real(wp) :: a0,a1,f00,f001,f2,deltaf,t0,t1,df(5)
    integer :: j

    ifail = 0
    t0 = 0.0_wp
    t1 = 0.0_wp
    if (h0 /= 0.0_wp) then
      if (x0+h0 /= 0.0_wp) then
        f00 = f1
      else
        h0 = 0.875_wp*h0
        f00 = me%f(x0+h0)
        if (me%stop) then
          ifail = -1
          return
        end if
      end if
      if (abs(h1) >= 32.0_wp*twoinf) h1 = h1/8.0_wp
      if (16.0_wp*abs(h1) > abs(h0)) h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
      f001 = me%f(x0+h0-h1)
      if (me%stop) then
        ifail = -1
        return
      end if
      if (f001 == f00) then
        if (256.0_wp*abs(h1) <= abs(h0)) then
          h1 = 2.0_wp*h1
          do
              f001 = me%f(x0+h0-h1)
              if (me%stop) then
                ifail = -1
                return
              end if
              if (f001 /= f00 .or. 256.0_wp*abs(h1) > abs(h0)) exit
              h1 = 2.0_wp*h1
          end do
          h1 = 8.0_wp*h1
        else
          h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
        end if
      else
        if (256.0_wp*twoinf <= abs(h0)) then
          do
              f001 = me%f(x0+h0-h1/2.0_wp)
              if (me%stop) then
                ifail = -1
                return
              end if
              if (f001 == f00 .or. abs(h1) < 4.0_wp*twoinf) exit
              h1 = h1/2.0_wp
          end do
          h1 = 8.0_wp*h1
          if (16.0_wp*abs(h1) > abs(h0)) h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
        else
          h1 = sign(h1,1.0_wp)*abs(h0)/16.0_wp
        end if
      end if
    else
      f00 = f0
    end if

    do j = 1,5
      f2 = me%f(x0+h0-real(2*j-1,wp)*h1)
      if (me%stop) then
        ifail = -1
        return
      end if
      df(j) = f2 - f00
      t0 = t0+df(j)
      t1 = t1+real(2*j-1,wp)*df(j)
    end do
    a0 = (33.0_wp*t0-5.0_wp*t1)/73.0_wp
    a1 = (-5.0_wp*t0+1.2_wp*t1)/73.0_wp
    facc = abs(a0)
    do j = 1,5
      deltaf = abs(df(j)-(a0+real(2*j-1,wp)*a1))
      if (facc < deltaf) facc = deltaf
    end do
    facc = 2.0_wp*facc

    end subroutine faccur