02_root_finding.f90

Go to the documentation of this file.

!----------------------------------------------------------------------------------------
!
! This file is part of EFTCAMB.
!
! Copyright (C) 2013-2016 by the EFTCAMB authors
!
! The EFTCAMB code is free software;
! You can use it, redistribute it, and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation;
! either version 3 of the License, or (at your option) any later version.
! The full text of the license can be found in the file eftcamb/LICENSE at
! the top level of the EFTCAMB distribution.
!
!----------------------------------------------------------------------------------------



!----------------------------------------------------------------------------------------


module eftcamb_rootfind

    use precision

    implicit none

    private

    public zbrac, zbrent

contains

    ! ---------------------------------------------------------------------------------------------
    subroutine zbrac( func, x1, x2, succes, funcZero )

        implicit none

        real(dl) func
        real(dl), intent(inout)  :: x1
        real(dl), intent(inout)  :: x2
        logical , intent(out)    :: succes
        real(dl), intent(in)     :: funcZero

        integer , parameter      :: ntry = 1000
        real(dl), parameter      :: FACTOR = 5._dl  

        real(dl) delta, temp, f1,f2
        integer  j
        external func

        ! initial check:
        if ( x1.eq.x2 ) stop 'you have to guess an initial range in zbrac'
        ! get an ordered interval:
        if (x2<x1) then
            temp=x2
            x2=x1
            x1=temp
        end if
        ! compute the function:
        f1 = func(x1) -funczero
        f2 = func(x2) -funczero

        succes=.true.
        do j=1, ntry
            ! check wether the interval is bracketing:
            if ( f1*f2 .lt. 0._dl ) return
            ! compute the interval and enlarge it:
            delta=abs(x2-x1)
            x1 = x1 -factor*delta
            x2 = x2 +factor*delta
            ! recompute the functions:
            f1 = func(x1) -funczero
            f2 = func(x2) -funczero
        end do
        succes=.false.

        return

    end subroutine zbrac

    ! ---------------------------------------------------------------------------------------------
    function zbrent(func,x1,x2,tol,funcZero,succes)

        implicit none

        real(dl) :: func
        real(dl) :: x1
        real(dl) :: x2
        real(dl) :: tol
        real(dl) :: funcZero
        logical  :: succes

        real(dl) :: zbrent

        external func

        integer ,parameter :: ITMAX = 2000
        real(dl),parameter :: EPS   = 3.d-15 
        integer iter
        real(dl) a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm

        succes = .true.
        a=x1
        b=x2
        fa=func(a)-funczero
        fb=func(b)-funczero
        if((fa.gt.0..and.fb.gt.0.).or.(fa.lt.0..and.fb.lt.0.)) then
            succes=.false.
            return
        end if
        c=b
        fc=fb
        do 11 iter=1,itmax
            if((fb.gt.0..and.fc.gt.0.).or.(fb.lt.0..and.fc.lt.0.))then
                c=a
                fc=fa
                d=b-a
                e=d
            endif
            if(abs(fc).lt.abs(fb)) then
                a=b
                b=c
                c=a
                fa=fb
                fb=fc
                fc=fa
            endif
            tol1=2.*eps*abs(b)+0.5*tol
            xm=.5*(c-b)
            if(abs(xm).le.tol1 .or. fb.eq.0.)then
                zbrent=b
                return
            endif
            if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
                s=fb/fa
                if(a.eq.c) then
                    p=2.*xm*s
                    q=1.-s
                else
                    q=fa/fc
                    r=fb/fc
                    p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
                    q=(q-1.)*(r-1.)*(s-1.)
                endif
                if(p.gt.0.) q=-q
                p=abs(p)
                if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
                    e=d
                    d=p/q
                else
                    d=xm
                    e=d
                endif
            else
                d=xm
                e=d
            endif
            a=b
            fa=fb
            if(abs(d) .gt. tol1) then
                b=b+d
            else
                b=b+sign(tol1,xm)
            endif
            fb=func(b)-funczero
11      continue
        succes = .false.
        zbrent=b
        return

    end function zbrent

end module eftcamb_rootfind

!----------------------------------------------------------------------------------------