solve Subroutine

private subroutine solve(me, ax, bx, xzero, fzero, iflag, fax, fbx, bisect_on_failure)

Main wrapper routine for all the methods.

Type Bound

root_solver

Arguments

Type IntentOptional Attributes Name
class(root_solver), intent(inout) :: me
real(kind=wp), intent(in) :: ax

left endpoint of initial interval
real(kind=wp), intent(in) :: bx

right endpoint of initial interval
real(kind=wp), intent(out) :: xzero

abscissa approximating a zero of f in the interval ax,bx
real(kind=wp), intent(out) :: fzero

value of f at the root (f(xzero))
integer, intent(out) :: iflag

status flag (-1=error, 0=root found, -4=ax must be /= bx)
real(kind=wp), intent(in), optional :: fax

if f(ax) is already known, it can be input here
real(kind=wp), intent(in), optional :: fbx

if f(bx) is already known, it can be input here
logical, intent(in), optional :: bisect_on_failure

if true, then if the specified method fails, it will be retried using the bisection method. (default is False). Note that this can use up to maxiter additional function evaluations.

Calls

proc~~solve~~CallsGraph proc~solve root_module::root_solver%solve find_root find_root proc~solve->find_root interface~root_scalar root_module::root_scalar proc~solve->interface~root_scalar proc~get_fa_fb root_module::root_solver%get_fa_fb proc~solve->proc~get_fa_fb proc~solution root_module::root_solver%solution proc~solve->proc~solution proc~root_scalar_by_name root_module::root_scalar_by_name interface~root_scalar->proc~root_scalar_by_name proc~root_scalar_by_type root_module::root_scalar_by_type interface~root_scalar->proc~root_scalar_by_type proc~root_scalar_by_name->interface~root_scalar proc~root_name_to_method root_module::root_name_to_method proc~root_scalar_by_name->proc~root_name_to_method proc~root_scalar_by_type->proc~solve proc~initialize_root_solver root_module::root_solver%initialize_root_solver proc~root_scalar_by_type->proc~initialize_root_solver proc~lowercase root_module::lowercase proc~root_name_to_method->proc~lowercase

Called by

proc~~solve~~CalledByGraph proc~solve root_module::root_solver%solve interface~root_scalar root_module::root_scalar proc~solve->interface~root_scalar proc~root_scalar_by_type root_module::root_scalar_by_type proc~root_scalar_by_type->proc~solve interface~root_scalar->proc~root_scalar_by_type proc~root_scalar_by_name root_module::root_scalar_by_name interface~root_scalar->proc~root_scalar_by_name proc~root_scalar_by_name->interface~root_scalar

Source Code

    subroutine solve(me,ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure)

    implicit none

    class(root_solver),intent(inout) :: me
    real(wp),intent(in)              :: ax      !! left endpoint of initial interval
    real(wp),intent(in)              :: bx      !! right endpoint of initial interval
    real(wp),intent(out)             :: xzero   !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
    real(wp),intent(out)             :: fzero   !! value of `f` at the root (`f(xzero)`)
    integer,intent(out)              :: iflag   !! status flag (`-1`=error, `0`=root found, `-4`=ax must be /= bx)
    real(wp),intent(in),optional     :: fax     !! if `f(ax)` is already known, it can be input here
    real(wp),intent(in),optional     :: fbx     !! if `f(bx)` is already known, it can be input here
    logical,intent(in),optional      :: bisect_on_failure !! if true, then if the specified method fails,
                                                          !! it will be retried using the bisection method.
                                                          !! (default is False). Note that this can use up
                                                        !! to `maxiter` additional function evaluations.

    real(wp) :: fa !! `f(ax)` passed to the lower level routine
    real(wp) :: fb !! `f(bx)` passed to the lower level routine

    if (ax==bx) then
        ! ax must be /= bx
        iflag = -4
        xzero = ax  ! just to return something
        fzero = fa  !
    else

        call me%get_fa_fb(ax,bx,fax,fbx,fa,fb)

        ! check trivial cases first:
        if (me%solution(ax,fa,xzero,fzero)) then

            iflag = 0

        else if (me%solution(bx,fb,xzero,fzero)) then

            iflag = 0

        else if (fa*fb>0.0_wp) then

            ! f(ax) and f(bx) do not have different signs
            iflag = -1
            xzero = ax  ! just to return something
            fzero = fa  !

        else

            ! call the root solver.
            ! make sure order is correct.
            if (ax<bx) then
                call me%find_root(ax,bx,fa,fb,xzero,fzero,iflag)
            else
                call me%find_root(bx,ax,fb,fa,xzero,fzero,iflag)
            end if

            ! if it failed, then we have the option to then try bisection
            if (iflag /= 0) then
                if (present(bisect_on_failure)) then
                    if (bisect_on_failure) then
                        ! use the wrapper routine for that with the input class
                        call root_scalar(root_method_bisection,&
                                         func_wrapper,ax,bx,xzero,fzero,iflag,&
                                         me%ftol,me%rtol,me%atol,me%maxiter,fa,fb,&
                                         bisect_on_failure = .false.)
                    end if
                end if
            end if

        end if

    end if

    contains

    function func_wrapper(x) result(f)
        !! wrapper function to use bisection
        implicit none
        real(wp),intent(in) :: x
        real(wp) :: f
        f = me%f(x)
    end function func_wrapper

    end subroutine solve