fzero Subroutine

subroutine fzero(f,b,c,iflag,rguess,tol_rel,tol_abs)
  !
  !    Input, real ( kind = 8 ) R, a (better) guess of a zero of F which could
  !    help in speeding up convergence.  If F(B) and F(R) have opposite signs, 
  !    a root will be found in the interval (B,R); if not, but F(R) and F(C) 
  !    have opposite signs, a root will be found in the interval (R,C); 
  !    otherwise, the interval (B,C) will be searched for a possible root.
  !    When no better guess is known, it is recommended that R be set to B or C;
  !    because if R is not interior to the interval (B,C), it will be ignored.
  !
  !    Input, real ( kind = 8 ) RE, the relative error used for RW in the 
  !    stopping criterion.  If the input RE is less than machine precision,
  !    then RW is set to approximately machine precision.
  !
  !    Input, real ( kind = 8 ) AE, the absolute error used in the stopping
  !    criterion.  If the given interval (B,C) contains the origin, then a
  !    nonzero value should be chosen for AE.
  !
  !    Output, integer IFLAG, a status code.  The user must check IFLAG 
  !    after each call.  Control returns to the user in all cases.
  !
  !    1, B is within the requested tolerance of a zero.
  !      the interval (b,c) collapsed to the requested
  !      tolerance, the function changes sign in (b,c), and
  !      f(x) decreased in magnitude as (b,c) collapsed.
  !
  !    2, F(B) = 0.  however, the interval (b,c) may not have
  !      collapsed to the requested tolerance.
  !
  !    3, B may be near a singular point of f(x).
  !      the interval (b,c) collapsed to the requested tolerance and 
  !      the function changes sign in (b,c), but
  !      f(x) increased in magnitude as (b,c) collapsed,i.e.
  !      max ( ABS ( f(b in) ), ABS ( f(c in) ) ) < ABS ( f(b out) ).
  !
  !    4, no change in sign of f(x) was found although the
  !      interval (b,c) collapsed to the requested tolerance.
  !      the user must examine this case and decide whether
  !      b is near a local minimum of f(x), or B is near a
  !      zero of even multiplicity, or neither of these.
  !
  !    5, too many (more than 500) function evaluations used.
  !
  implicit none
  interface                                                             
     function f(x)                                                   
       real(8),intent(in) :: x
       real(8)            :: f
     end function f
  end interface
  real(8),optional :: rguess
  real(8),optional :: tol_rel
  real(8),optional :: tol_abs
  real(8) :: re
  real(8) :: ae
  real(8) :: a
  real(8) :: acbs
  real(8) :: acmb
  real(8) :: aw
  real(8) :: b
  real(8) :: c
  real(8) :: cmb
  real(8) :: er
  ! real(8) ::, external :: f
  real(8) :: fa
  real(8) :: fb
  real(8) :: fc
  real(8) :: fx
  real(8) :: fz
  integer :: ic
  integer :: iflag
  integer :: kount
  real(8) :: p
  real(8) :: q
  real(8) :: r
  real(8) :: rw
  real(8) :: t
  real(8) :: tol
  real(8) :: z
  !
  re=1d-12;if(present(tol_rel))re=tol_rel
  ae=1d-9;if(present(tol_abs))ae=tol_abs
  r = c;if(present(rguess))r=rguess
  !   
  er = 2.0D+00 * epsilon ( er )
  !
  !  Initialize.
  !
  z = r
  if ( r <= min ( b, c ) .or. max ( b, c ) <= r ) then
     z = c
  end if
  rw = max ( re, er )
  aw = max ( ae, 0.0D+00 )
  ic = 0
  t = z
  fz = f(t)
  fc = fz
  t = b
  fb = f(t)
  kount = 2
  if ( sign ( 1.0D+00, fz ) /= sign ( 1.0D+00, fb ) ) then
     c = z
  else if ( z /= c ) then
     t = c
     fc = f(t)
     kount = 3
     if ( sign ( 1.0D+00, fz ) /= sign ( 1.0D+00, fc ) ) then
        b = z
        fb = fz
     end if
  end if
  a = c
  fa = fc
  acbs = abs ( b - c )
  fx = max ( abs ( fb ), abs ( fc ) )
  do
     !
     !  Interchange
     !
     if ( abs ( fc ) < abs ( fb ) ) then
        a = b
        fa = fb
        b = c
        fb = fc
        c = a
        fc = fa
     end if
     cmb = 0.5D+00 * ( c - b )
     acmb = abs ( cmb )
     tol = rw * abs ( b ) + aw
     !
     !  Test stopping criterion and function count.
     !
     if ( acmb <= tol ) then
        exit
     end if
     if ( fb == 0.0D+00 ) then
        iflag = 2
        return
     end if
     if ( 500 <= kount ) then
        iflag = 5
        return
     end if
     !
     !  Calculate new iterate implicitly as b+p/q
     !  where we arrange 0 <= p.
     !  The implicit form is used to prevent overflow.
     !
     p = ( b - a ) * fb
     q = fa - fb
     if ( p < 0.0D+00 ) then
        p = -p
        q = -q
     end if
     !
     !  Update A and check for satisfactory reduction
     !  in the size of the bracketing interval.
     !  If not, perform bisection.
     !
5    continue
     a = b
     fa = fb
     ic = ic + 1
     if ( ic < 4 ) then
        go to 6
     end if
     if ( acbs <= 8.0D+00 * acmb ) then
        b = b + cmb
        go to 9
     end if
     ic = 0
     acbs = acmb
     !
     !  Test for too small a change
     !
6    continue
     if ( abs ( q ) * tol < p ) then
        go to 7
     end if
     !
     !  Increment by tolerance
     !
     b = b + sign ( tol, cmb )
     go to 9
     !
     !  Root ought to be between b and (c+b)/2.
     !
7    continue
     !
     !  Use the secant rule or bisection.
     !
     if ( p < cmb * q ) then
        b = b + p / q
     else
        b = b + cmb
     end if
     !
     !  Have completed computation for new iterate B.
     !
9    continue
     t = b
     fb = f(t)
     kount = kount + 1
     !
     !  Decide whether next step is interpolation or extrapolation.
     !
     if ( sign ( 1.0D+00, fb ) == sign ( 1.0D+00, fc ) ) then
        c = a
        fc = fa
     end if
  end do
  !
  !  Finished.  Process results for proper setting of IFlAG.
  !
  if ( sign ( 1.0D+00, fb ) == sign ( 1.0D+00, fc ) ) then
     iflag = 4
     return
  end if
  if ( fx < abs ( fb ) ) then
     iflag = 3
     return
  end if
  iflag = 1    
  return
end subroutine fzero