FIXME : Add documentation
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer, | intent(in) | :: | n | |||
double precision, | intent(in) | :: | xold | |||
double precision, | intent(in) | :: | xnew | |||
double precision, | intent(in) | :: | pold | |||
double precision, | intent(in) | :: | pnew | |||
double precision, | intent(out) | :: | zz |
subroutine find_zero (n, xold, xnew, pold, pnew, zz)
use file_utils, only: error_unit
implicit none
integer, intent (in) :: n
double precision, intent (in) :: xold, xnew, pold, pnew
double precision, intent (out) :: zz
integer :: i, maxit=100
real :: eps
double precision :: x1, x2, p1, p2, pz
! <doc>
! eps is declared as real on purpose.
! [see comment in find_zero_bisect_newton above.]
! </doc>
eps = sqrt(epsilon(eps)*epsilon(x1)) * 2.0
x1 = xold
x2 = xnew
p1 = pold
p2 = pnew
if (debug) write (*,'(a,4es15.5e3)') 'initial ', x1, p1, x2, p2
! bisection
do i=1, 8
zz = (x1+x2) * dble(.5)
pz = laguerre_l(n,zz)
if (abs(pz) <= epsilon(dble(0.0))) return
if (pz*p1 < 0.) then
p2=pz ; x2=zz
else
p1=pz ; x1=zz
end if
if (debug) write (*,'(a,4es15.5e3)') 'bisection ', x1, p1, x2, p2
end do
! newton-raphson
if (exactly_equal(zz, x1)) x1 = x2
do i=1, maxit
x1 = zz
p1 = laguerre_l(n,x1)
zz = x1 - p1 / laguerre_lp(n,x1)
pz = laguerre_l(n,zz)
if (debug) write (*,'(a,4es15.5e3)') &
'newton ', zz, pz, x1, p1
if (min(abs(zz/x1-1.0), abs(pz)) < eps) exit
end do
if (i==maxit+1) write (error_unit(),*) &
& 'WARNING: too many iterations in get_laguerre_grids'
end subroutine find_zero