!=======================================================================
! ch11_bisection.f90
!
! Find one root of f(x) = 0 by bisection on an initial interval.
!
! Origin:
!   Based on the bisection program by Alexander Godunov, January 2010.
!   The numerical algorithm is retained. Comments and failure reporting
!   have been revised for the Chapter 11 companion website.
!
! Example:
!   f(x) = sqrt(x) - cos(x) on [0,1].
!   Reference root = 0.64171437087288... .
!
! The program reports the iteration count, root estimate, and f(root).
! The requested tolerance applies to the width of the final bracket.
! A small bracket by itself does not establish continuity or exclude a
! singularity, so inspect the residual and the function's domain.
!=======================================================================
program main
implicit none
double precision f,x1,x2,root,eps
integer flag
external f

! Initial bracket and requested maximum bracket width
x1 = 0.0
x2 = 1.0
eps = 1.0e-6

write (*,100)

  call bisection(f,x1,x2,eps,root,flag)

! Report a missing bracket or invalid input before printing a root
if(flag == 0) then
  write(*,104)
  stop
end if

write(*,101) abs(flag), eps
if(flag > 0) then
    write(*,102) Root
    else
    write(*,103)
    write(*,102) Root
end if

! Check the residual of the reported estimate
write(*,105) f(Root)

100 format(' A single root of f(x) ',/, &
           ' Method - Bisection')
101 format(' iterations  = ',i3,/,&
           ' tolerance   = ',1pe12.5)
102 format(' root        = ',1pe12.5)
103 format(' iteration limit reached before the requested width')
104 format(' invalid interval or no initial sign change')
105 format(' f(root)     = ',1pe12.5)

end program main


   Subroutine bisection(f,x1,x2,eps,Root,flag)
!=======================================================================
! Locate one root of f(x) on [x1,x2] by repeatedly halving the bracket.
!
! Input:
!   f     - function to evaluate, continuous on the initial bracket
!   x1    - left endpoint of the initial bracket
!   x2    - right endpoint of the initial bracket
!   eps   - requested maximum width of the final bracket (positive)
!
! Output:
!   Root  - midpoint of the final bracket when flag is nonzero
!   flag  - convergence indicator and iteration count
!           flag > 0: width <= eps; flag is the iteration count
!           flag = 0: invalid interval/tolerance or no sign change
!           flag < 0: iteration limit reached; abs(flag) is the count
!
! Notes:
!   A sign change guarantees a root only when f is continuous.
!   The maximum number of iterations is 200. If the requested width
!   is not reached by then, the estimate is not marked converged.
!   The algorithm uses bracket width, not |f(Root)|, to stop.
!=======================================================================
implicit none
double precision f, x1, x2, eps, Root
double precision a, b, c
integer i, flag
integer, parameter:: iter=200

! Check that the inputs describe a usable initial bracket
if(x1 >= x2 .or. eps <= 0.0d0) then
  flag = 0
  return
end if
if(f(x1)*f(x2)>0.0) then
  flag = 0
  return
end if

! Initialize the current bracket
a=x1
b=x2

! Repeatedly retain the half whose endpoints have opposite signs
do i=1,iter
  ! This midpoint form avoids overflow in a+b when a and b are large.
  c=a+(b-a)/2.0
  if(f(a)*f(c).le.0.0) then
      b = c
    else
      a=c
  end if
  ! Stop when the current bracket is sufficiently narrow.
  if(abs(b-a)<= eps) exit  
end do
Root=(b+a)/2.0

! Report whether the requested bracket width was reached
if (i > iter) then
  flag = -iter
else
  flag = i
end if
end subroutine bisection

!-----------------------------------------------------------------------
  Function f(x)
implicit none
double precision f, x
! Example function. The chosen bracket must stay in the domain x >= 0.
! To use the alternative below, change the bracket to [-1,0].
!  f = x + cos(x)
  f = sqrt(x)-cos(x)
end function f