!*******************************************************
!*  Program to demonstrate the Newton root subroutine  *
!*             in the complex domain                   *
!* --------------------------------------------------- *
!*  Reference; BASIC Scientific Subroutines, Vol. II   *
!*  By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
!*                                                     *
!*                 F90 version by J-P Moreau, Paris.   *
!*                         (www.jpmoreau.fr)           *
!* --------------------------------------------------- *
!* SAMPLE RUN:                                         *
!*                                                     *
!* (Example: find a complex root of  z^2 + 1 = 0)      *
!*                                                     *
!* What is the initial guess:                          *
!*                                                     *
!*    X0 = .2                                          *
!*    Y0 = 1.2                                         *
!*                                                     *
!* Convergence criterion: 1e-8                         *
!* Maximum number of iterations: 10                    *
!*                                                     *
!* The root estimate is:                               *
!*                                                     *
!*    X0 =  0.000000                                   *
!*    Y0 =  1.000000                                   *
!*                                                     *
!* The number of iterations performed was: 5           *
!*                                                     *
!*******************************************************
PROGRAM DEMO_ZNEWTON

integer k,n
real*8  a,e,x0,y0,x,y

  print *,' '
  print *,'What is the initial guess:'
  print *,' '
  write(*,"('    X0 = ')",advance='no')
  read *, x0
  write(*,"('    Y0 = ')",advance='no')
  read *, y0
  print *,' '
  write(*,"(' Convergence criterion: ')",advance='no')
  read *, e
  write(*,"(' Maximum number of iterations: ')",advance='no')
  read *, n

  call ZNewton(n,k,a,e,x0,y0,x,y)      ! Call ZNewton subroutine

  if (a.ne.0.d0) then
    print *,' '
    print *,'The root estimate is:'
    print *,' '
    write(*,50) x
    write(*,60) y
    write(*,70) k
  end if

  stop

50 format('    X0 = ',f9.6)
60 format('    Y0 = ',f9.6)
70 format(/' The number of iterations performed was: ',I2//)

end

!*********************************************************
!  Functions subroutine
! --------------------------------------------------------
Subroutine Eval(x,y,u,v,u1,v1,u2,v2)
real*8 x,y,u,v,u1,v1,u2,v2
  u=x*x-y*y+1.d0; v=2.d0*x*y
  u1=2.d0*x; u2=-2.d0*y
  v1=2.d0*y; v2=2.d0*x
  return
end
!*********************************************************


!*************************************************
!*  Complex root seeking using Newton's method   *
!* --------------------------------------------- *
!* This routine uses the complex domain form of  *
!* Newton's method for iteratively searching     *
!* for roots. The complex function and its first *
!* partial derivatives must be available in the  *
!* form: F(Z) = U(X,Y) + I V(X,Y). The required  *
!* derivatives are DU/DX and DU/DY.              *
!* --------------------------------------------- *
!* INPUTS; initial guess, x0 and y0, convergence *
!* criteria, e, maximum number of iterations, n. *
!* OUTPUTS; approximation to the root, x and y,  *
!* number of performed iterations, k.            *
!*************************************************
Subroutine ZNewton(n,k,a,e,x0,y0,x,y)
  !Label: 100
  parameter(TINY=1.d-12)  !small number
  integer n,k
  real*8 e,x0,y0,x,y
  real*8 u,v,u1,v1,u2,v2
  k=0 
100 k=k+1
  ! Get u,v and the derivatives u1,v1,u2,v2
  call Eval(x0,y0,u,v,u1,v1,u2,v2)
  a=u1*u1+u2*u2
  ! Guard against a=0
  if (a < TINY) then
    print *,' '
    print *,'ZERO DIVIDE ERROR - u1*u1+u2*u2 must be <> 0.'
    print *,' '
    return
  end if
  x=x0+(v*u2-u*u1)/a
  y=y0-(v*u1+u*u2)/a
  ! Check for convergence in euclidean space
  if ((x0-x)*(x0-x)+(y0-y)*(y0-y)<=e*e) return
  if (k>=n) return
  x0=x ; y0=y
  goto 100
end ! ZNewton()


! End of file znewton.f90