!=======================================================================
! ch11_secant.f90
!
! Find one root of f(x) = 0 by the secant method.
!
! Origin:
!   Based on Secant.f90 by Alexander Godunov, January 2010.
!   The secant update and residual stopping test are retained.
!   The two starting guesses are now explicit. An old large-step
!   midpoint substitution has been removed because it was not a
!   bracketed bisection step or an ordinary secant step.
!
! Example:
!   f(x) = x + cos(x), starting from x1 = -1 and x2 = -0.9.
!   Reference root = -0.73908513321516... .
!
! The program reports the iteration count, root estimate, and f(root).
! The secant method does not preserve a bracket. Starting values can
! affect convergence and may lead to a different root or failure.
!=======================================================================
program main
implicit none
double precision f,x1,x2,root,eps
integer flag
external f

! Two initial guesses and requested residual tolerance
x1  = -1.0d0
x2  = -0.9d0
eps = 1.0d-7

write(*,100) x1,x2
call Secant(f,x1,x2,eps,root,flag)

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

write(*,101) eps,abs(flag)
if(flag < 0) write(*,103)
write(*,102) root

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

100 format(' A single root of f(x)',/,&
           ' Method - Secant',/,&
           ' Initial guesses = ',f8.3,', ',f8.3)
101 format(' residual tolerance = ',1pe12.5,/,&
           ' iterations         = ',i3)
102 format(' root               = ',1pe12.5)
103 format(' iteration limit reached before the requested residual')
104 format(' invalid tolerance/guesses or zero secant slope')
105 format(' f(root)            = ',1pe12.5)

end program main


subroutine Secant(f,x1,x2,eps,Root,flag)
!=======================================================================
! Locate one root using secant lines through successive iterates.
!
! Input:
!   f     - function evaluating f(x)
!   x1    - first initial guess
!   x2    - second, distinct initial guess
!   eps   - requested maximum residual |f(Root)| (positive)
!
! Output:
!   Root  - most recent secant estimate when flag is nonzero
!   flag  - convergence indicator and iteration count
!           flag > 0: |f(Root)| <= eps; flag is the iteration count
!           flag = 0: invalid input or zero secant slope
!           flag < 0: iteration limit reached; abs(flag) is the count
!
! Notes:
!   This is the ordinary secant method without damping or bracketing.
!   The maximum number of iterations is 200. An initial guess that
!   already satisfies the residual test is returned with flag=1 to
!   follow the positive-success convention, without an update.
!=======================================================================
implicit none
double precision f,x1,x2,eps,Root
integer flag,i
integer, parameter :: iter=200
double precision a,b,c,fa,fb,fc,df
external f

flag = 0
if(eps <= 0.0d0 .or. x1 == x2) return

a = x1
b = x2
fa = f(a)
fb = f(b)

! An accurate starting guess needs no secant update.
if(abs(fa) <= eps) then
  Root = a
  flag = 1
  return
end if
if(abs(fb) <= eps) then
  Root = b
  flag = 1
  return
end if

! Form the secant slope from the two most recent function values.
do i=1,iter
  df = (fb-fa)/(b-a)
  if(df == 0.0d0) return
  c = b - fb/df
  fc = f(c)

  ! Stop on the residual, as in the original program.
  if(abs(fc) <= eps) then
    Root = c
    flag = i
    return
  end if

  a = b
  fa = fb
  b = c
  fb = fc
  if(a == b) return
end do

Root = c
flag = -iter
end subroutine Secant


function f(x)
!-----------------------------------------------------------------------
! Example function. To use another problem, change this expression.
!-----------------------------------------------------------------------
implicit none
double precision f,x
f = x + cos(x)
end function f