FIXME : Add documentation
eps is declared as real on purpose. We don't require higher order precision for convergence test because of a performance reason. The following definition means eps is a geometric mean of the machine-epsilons in real and double precisions. (note that real/double are promoted to double/double or double/quad depending on the compiler.)
[same applies to eps in find_zero below.]
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_bisect_newton_radau (n, xold, xnew, pold, pnew, zz)
use mp, only: mp_abort
use file_utils, only: error_unit
implicit none
integer, intent (in) :: n
double precision, intent (in) :: xold, xnew, pold, pnew
! real, intent (in) :: eps
double precision, intent (out) :: zz
integer :: i, maxit=100
real :: eps
double precision :: x1, x2, p1, p2, pz
eps = sqrt(epsilon(eps)*epsilon(x1)) * 2.0
i=0
x1 = xold
x2 = xnew
p1 = pold
p2 = pnew
if (debug) write(*,'(4f10.5)') x1, p1, x2, p2
! bisection
do i=1, 8
zz = (x1+x2) * dble(.5)
pz = radau_p(n,zz)
if (abs(pz) <= epsilon(pz)) return
if (pz*p1 < 0.0) then
p2=pz ; x2=zz
else
p1=pz ; x1=zz
end if
if (debug) write(*,'(4f10.5)') x1, p1, x2, p2
end do
if (debug) print*, 'find_zero_bisect_newton_radau: finished bisection'
! newton-raphson
if (exactly_equal(zz, x1)) x1 = x2
do i=1, maxit
x1 = zz
p1 = radau_p(n,x1)
zz = x1 - p1 / radau_pp(n,x1)
pz = radau_p(n,zz)
if (debug) write (*,'(4f10.5)') 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 find_zero_bisect_newton_radau'
if (debug) call mp_abort("Abort in gauss_quad:find_zero_bisect_newton_radau as debug=.true.")
end subroutine find_zero_bisect_newton_radau