root_finding_m.f90 Source File

module root_finding_m
    !! Routines for finding zeros
    use precision_m, only: FP, DP
    implicit none  
    
    type, public, abstract :: func_1D_t
        !! Defines function type for find_zero
     contains
       procedure(func_if), deferred, public :: func
    end type func_1D_t

    abstract interface
        real(DP) function func_if(this, t, iuser, ruser)
            !! Defines interface of function for find_zero
             import DP, func_1D_t
             class(func_1D_t), intent(in) :: this
             !! Instance of type
             real(DP), intent(in) :: t
             !! Abscissa
             integer, intent(in) :: iuser(:)
             !! Integer parameters of function
             real(DP), intent(in) :: ruser(:)
             !! Real parameters of function f
        end function
    end interface

    private :: faux

    ! Pointers for passing function and its parameters 
    ! to auxiliary function faux. These are threadprivate to
    ! make the usage of the find_zero routine threadsafe
    class(func_1D_t), private, pointer :: func_mw => null()
    integer, contiguous, private, pointer :: iuser_mw(:) => null()
    real(DP), contiguous, private, pointer :: ruser_mw(:) => null()
    !$omp threadprivate(func_mw, iuser_mw, ruser_mw)
        
contains

    subroutine find_zero(a, b, tol, f, iuser, ruser, xzero, ifail)
        !! Finds the zero of function f within interval [a,b]
        !! Wrapper for Fortran 77 routine in zeroin.f
        real(DP), intent(in) :: a
        !! Left limit of interval
        real(DP), intent(in) :: b
        !! Right limit of interval
        real(DP), intent(in) :: tol
        !! Absolute tolerance of result  
        class(func_1D_t), target, intent(in) :: f
        !! Function f
        integer, contiguous, target, intent(in) :: iuser(:)
        !! Integer parameters of function f
        real(DP), contiguous, target, intent(in) :: ruser(:)
        !! Real parameters of function f
        real(DP), intent(out) :: xzero
        !! On success: zero of function f
        integer, intent(out) :: ifail
        !! On success = 0, otherwise -1

        real(DP) :: zeroin, sgn

        func_mw => f
        iuser_mw => iuser
        ruser_mw => ruser

        ! check if f(1) has different sign as f(b)
        sgn = faux(a) * faux(b)
        if (sgn >= 0.0_DP) then
            ifail = -1
            return
        endif

        ! Call of Fortran 77 routine in zeroin.f
        xzero = zeroin(a, b, faux, 0.5_DP * tol)
        ifail = 0

        func_mw => null()
        iuser_mw => null()
        ruser_mw => null()

    end subroutine

    real(DP) function faux(t)
        !! Auxiliary function with only abscissa in interface
        real(DP), intent(in) :: t
        !! Abscissa
        
        faux = func_mw%func(t, iuser_mw, ruser_mw)
        
    end function
    
end module