find_zero Subroutine

private subroutine find_zero(n, xold, xnew, pold, pnew, zz)

Uses

FIXME : Add documentation

Arguments

Type IntentOptional 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

Contents

Source Code


Source Code

  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