subroutine root(f, fval, a, b, xerrbi, xerrsec, ier, soln)
implicit none
interface
real function f(x_)
implicit none
real, intent(in) :: x_
end function f
end interface
real, intent(in) :: fval, a, b, xerrbi, xerrsec
real, intent(out) :: soln
integer, intent(out) :: ier
real :: a1, b1, f1, f2, f3, trial, aold
integer :: i, niter
logical :: skip_bisection
ier = 0 ; a1 = a ; b1 = b
f1 = f(a1) - fval ; f2 = f(b1) - fval
skip_bisection = xerrbi <= 0 !Could be xerrbi <= abs(b1 - a1)?
if (.not. skip_bisection .and. f1 * f2 > 0.) then
write(*, *) 'f1 and f2 have same sign in bisection routine root'
write(*, *) 'a1,f1,b1,f2=', a1, f1, b1, f2
write(*, *) 'fval=', fval
ier = 1
skip_bisection = .true.
end if
if (.not. skip_bisection) then
! Presumably we require |b1 - a1| >= xerrbi here to avoid niter < 1?
niter = 1 + int(log(abs(b1 - a1) / xerrbi) / log(2.))
! Repeatedly bisect between limits, swapping one limit
! for the bisection point each iteration. Helps to narrow
! down the bracketing region substantially.
do i = 1, niter
trial = 0.5 * (a1 + b1)
f3 = f(trial) - fval
if (f3 * f1 > 0.) then
a1 = trial
f1 = f3
else
b1 = trial
f2 = f3
end if
end do
end if
!Swap so |f1| < |f2|
if( abs(f1) > abs(f2) ) then
f3 = f1 ; f1 = f2 ; f2 = f3
aold = a1 ; a1 = b1 ; b1 = aold
end if
! start secant method
do i = 1, 10
aold = a1 ; f3 = f2 - f1
if (abs(f3) < 1.e-11) f3 = 1.e-11 !Should we exit here? Saying f2~f1
a1 = a1 - f1 * (b1 - a1) / f3 !Update guess for root.
f2 = f1 ; b1 = aold ; f1 = f(a1) - fval
if (abs(a1 - b1) < xerrsec) exit !Root found closely enough so exit
end do
soln = a1
! Detect instances where we have not converged to a solution
if (abs(f1) > xerrsec) ier = ier + 2
end subroutine root