nonlin_solve.f90

! nonlin_solve.f90
 
module nonlin_solve
    use, intrinsic :: iso_fortran_env, only : int32, real64
    use nonlin_constants
    use nonlin_core
    use nonlin_linesearch, only : line_search, limit_search_vector
    use ferror
    use linalg, only : qr_factor, form_qr, qr_rank1_update, lu_factor, &
        rank1_update, mtx_mult, recip_mult_array, solve_triangular_system, &
        solve_lu
    implicit none
    private
    public :: line_search_solver
    public :: quasi_newton_solver
    public :: newton_solver
    public :: brent_solver
    public :: newton_1var_solver
    public :: test_convergence
 
! ******************************************************************************
! NONLIN_SOLVE_LINE_SEARCH.F90
! ------------------------------------------------------------------------------
    type, abstract, extends(equation_solver) :: line_search_solver
        private
        class(line_search), allocatable :: m_linesearch
        logical :: m_uselinesearch = .true.
    contains
        procedure, public :: get_line_search => lss_get_line_search
        procedure, public :: set_line_search => lss_set_line_search
        procedure, public :: set_default_line_search => lss_set_default
        procedure, public :: is_line_search_defined => &
            lss_is_line_search_defined
        procedure, public :: get_use_line_search => lss_get_use_search
        procedure, public :: set_use_line_search => lss_set_use_search
    end type
 
! ------------------------------------------------------------------------------
    interface
        module subroutine lss_get_line_search(this, ls)
            class(line_search_solver), intent(in) :: this
            class(line_search), intent(out), allocatable :: ls
        end subroutine
 
        module subroutine lss_set_line_search(this, ls)
            class(line_search_solver), intent(inout) :: this
            class(line_search), intent(in) :: ls
        end subroutine
 
        module subroutine lss_set_default(this)
            class(line_search_solver), intent(inout) :: this
            type(line_search) :: ls
        end subroutine
 
        pure module function lss_is_line_search_defined(this) result(x)
            class(line_search_solver), intent(in) :: this
            logical :: x
        end function
 
        pure module function lss_get_use_search(this) result(x)
            class(line_search_solver), intent(in) :: this
            logical :: x
        end function
 
        module subroutine lss_set_use_search(this, x)
            class(line_search_solver), intent(inout) :: this
            logical, intent(in) :: x
        end subroutine
    end interface
 
! ******************************************************************************
! NONLIN_SOLVE_QUASI_NEWTON.F90
! ------------------------------------------------------------------------------
    type, extends(line_search_solver) :: quasi_newton_solver
        private
        integer(int32) :: m_jdelta = 5
    contains
        procedure, public :: solve => qns_solve
        procedure, public :: get_jacobian_interval => qns_get_jac_interval
        procedure, public :: set_jacobian_interval => qns_set_jac_interval
    end type
 
! ------------------------------------------------------------------------------
    interface
        module subroutine qns_solve(this, fcn, x, fvec, ib, err)
            class(quasi_newton_solver), intent(inout) :: this
            class(vecfcn_helper), intent(in) :: fcn
            real(real64), intent(inout), dimension(:) :: x
            real(real64), intent(out), dimension(:) :: fvec
            type(iteration_behavior), optional :: ib
            class(errors), intent(inout), optional, target :: err
        end subroutine
 
        pure module function qns_get_jac_interval(this) result(n)
            class(quasi_newton_solver), intent(in) :: this
            integer(int32) :: n
        end function
 
        module subroutine qns_set_jac_interval(this, n)
            class(quasi_newton_solver), intent(inout) :: this
            integer(int32), intent(in) :: n
        end subroutine
    end interface
 
! ******************************************************************************
! ------------------------------------------------------------------------------
    type, extends(line_search_solver) :: newton_solver
    contains
        procedure, public :: solve => ns_solve
    end type
 
! ------------------------------------------------------------------------------
    interface
        module subroutine ns_solve(this, fcn, x, fvec, ib, err)
            class(newton_solver), intent(inout) :: this
            class(vecfcn_helper), intent(in) :: fcn
            real(real64), intent(inout), dimension(:) :: x
            real(real64), intent(out), dimension(:) :: fvec
            type(iteration_behavior), optional :: ib
            class(errors), intent(inout), optional, target :: err
        end subroutine
    end interface
 
! ******************************************************************************
! NONLIN_SOLVE_BRENT.F90
! ------------------------------------------------------------------------------
    type, extends(equation_solver_1var) :: brent_solver
    contains
        procedure, public :: solve => brent_solve
    end type
 
! ------------------------------------------------------------------------------
    interface
        module subroutine brent_solve(this, fcn, x, lim, f, ib, err)
            class(brent_solver), intent(inout) :: this
            class(fcn1var_helper), intent(in) :: fcn
            real(real64), intent(inout) :: x
            type(value_pair), intent(in) :: lim
            real(real64), intent(out), optional :: f
            type(iteration_behavior), optional :: ib
            class(errors), intent(inout), optional, target :: err
        end subroutine
    end interface
 
! ******************************************************************************
! NONLIN_SOLVE_NEWTON1VAR.F90
! ------------------------------------------------------------------------------
    type, extends(equation_solver_1var) :: newton_1var_solver
    contains
        procedure, public :: solve => newt1var_solve
    end type
! ------------------------------------------------------------------------------
    interface
        module subroutine newt1var_solve(this, fcn, x, lim, f, ib, err)
            class(newton_1var_solver), intent(inout) :: this
            class(fcn1var_helper), intent(in) :: fcn
            real(real64), intent(inout) :: x
            type(value_pair), intent(in) :: lim
            real(real64), intent(out), optional :: f
            type(iteration_behavior), optional :: ib
            class(errors), intent(inout), optional, target :: err
        end subroutine
    end interface
 
contains
! ******************************************************************************
! GENERAL ROUTINES
! ------------------------------------------------------------------------------
    subroutine test_convergence(x, xo, f, g, lg, xtol, ftol, gtol, c, cx, cf, &
            cg, xnorm, fnorm)
        ! Arguments
        real(real64), intent(in), dimension(:) :: x, xo, f, g
        real(real64), intent(in) :: xtol, ftol, gtol
        logical, intent(in) :: lg
        logical, intent(out) :: c, cx, cf, cg
        real(real64), intent(out) :: xnorm, fnorm
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
        real(real64), parameter :: one = 1.0d0
        real(real64), parameter :: half = 0.5d0
 
        ! Local Variables
        integer(int32) :: i, nvar, neqn
        real(real64) :: test, dxmax, fc, den
 
        ! Initialization
        nvar = size(x)
        neqn = size(f)
        cx = .false.
        cf = .false.
        cg = .false.
        c = .false.
        fc = half * dot_product(f, f)
        fnorm = zero
        xnorm = zero
 
        ! Test for convergence on residual
        do i = 1, neqn
            fnorm = max(abs(f(i)), fnorm)
        end do
        if (fnorm < ftol) then
            cf = .true.
            c = .true.
            return
        end if
 
        ! Test the change in solution
        do i = 1, nvar
            test = abs(x(i) - xo(i)) / max(abs(x(i)), one)
            xnorm = max(test, xnorm)
        end do
        if (xnorm < xtol) then
            cx = .true.
            c = .true.
            return
        end if
 
        ! Test for zero gradient slope - do not set convergence flag
        if (lg) then
            test = zero
            den = max(fc, half * nvar)
            do i = 1, nvar
                dxmax = abs(g(i)) * max(abs(x(i)), one) / den
                test = max(test, dxmax)
            end do
            if (test < gtol) then
                cg = .true.
            end if
        end if
    end subroutine
 
! ------------------------------------------------------------------------------
end module