!=======================================================================
! ch11_newton.f90
!
! Find one root of f(x) = 0 by Newton's method.
!
! Origin:
!   Based on Newton1.f90 by Alexander Godunov, January 2010.
!   The Newton update and residual stopping test are retained.
!   An old large-step substitution has been removed because it was
!   not a Newton or bisection step.
!
! Example:
!   f(x) = x + cos(x), starting from x0 = -1.
!   Reference root = -0.73908513321516... .
!
! The program reports the iteration count, root estimate, and f(root).
! Newton's method is local: a different initial guess may converge to
! another root or fail. A small residual should be interpreted in the
! context of the function and the physical problem.
!=======================================================================
program main
implicit none
double precision f,x0,root,eps
integer flag
external f

! Initial guess and requested residual tolerance
x0  = -1.0d0
eps = 1.0d-7

write(*,100) x0
call Newton1(f,x0,eps,root,flag)

! Report invalid input or a zero derivative 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,0)

100 format(' A single root of f(x)',/,&
           ' Method - Newton',/,&
           ' Initial guess = ',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 or zero derivative at an iterate')
105 format(' f(root)            = ',1pe12.5)

end program main


subroutine Newton1(f,x0,eps,Root,flag)
!=======================================================================
! Locate one root near x0 using Newton's update x_new=x-f(x)/f'(x).
!
! Input:
!   f     - function evaluating f(x) for k=0 and f'(x) for k=1
!   x0    - initial guess
!   eps   - requested maximum residual |f(Root)| (positive)
!
! Output:
!   Root  - most recent Newton estimate when flag is nonzero
!   flag  - convergence indicator and iteration count
!           flag > 0: |f(Root)| <= eps; flag is the iteration count
!           flag = 0: invalid tolerance or zero derivative
!           flag < 0: iteration limit reached; abs(flag) is the count
!
! Notes:
!   This is the ordinary Newton method without damping or bracketing.
!   The maximum number of iterations is 200. If the initial guess
!   already satisfies the residual test, flag=1 denotes success even
!   though no Newton update was needed.
!=======================================================================
implicit none
double precision f,x0,eps,Root
integer flag,i
integer, parameter :: iter=200
double precision x,xnew,fx,dfx
external f

flag = 0
if(eps <= 0.0d0) return
x = x0

! A sufficiently accurate initial guess needs no Newton update.
fx = f(x,0)
if(abs(fx) <= eps) then
  Root = x
  flag = 1
  return
end if

! Repeatedly use the value and derivative at the current estimate.
do i=1,iter
  dfx = f(x,1)
  if(dfx == 0.0d0) return

  xnew = x - fx/dfx
  fx = f(xnew,0)

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

Root = xnew
flag = -iter
end subroutine Newton1


function f(x,k)
!-----------------------------------------------------------------------
! Evaluate the example function for k=0 and its derivative for k=1.
! To use another function, change both expressions consistently.
!-----------------------------------------------------------------------
implicit none
double precision f,x
integer k

if(k == 0) f = x + cos(x)
if(k == 1) f = 1.0d0 - sin(x)
end function f