nonlin_linesearch.f90

! nonlin_linesearch.f90
 
! REFERENCES
! - https://scicomp.stackexchange.com/questions/26330/backtracking-armijo-line-search-algorithm
! - https://ctk.math.ncsu.edu/
 
module nonlin_linesearch
    use, intrinsic :: iso_fortran_env, only : int32, real64
    use nonlin_constants
    use nonlin_core
    use ferror, only : errors
    implicit none
    private
    public :: line_search
    public :: limit_search_vector
 
! ******************************************************************************
! TYPES
! ------------------------------------------------------------------------------
    type line_search
        private
        integer(int32) :: m_maxeval = 100
        real(real64) :: m_alpha = 1.0d-4
        real(real64) :: m_factor = 0.1d0
    contains
        procedure, public :: get_max_fcn_evals => ls_get_max_eval
        procedure, public :: set_max_fcn_evals => ls_set_max_eval
        procedure, public :: get_scaling_factor => ls_get_scale
        procedure, public :: set_scaling_factor => ls_set_scale
        procedure, public :: get_distance_factor => ls_get_dist
        procedure, public :: set_distance_factor => ls_set_dist
        generic, public :: search => ls_search_mimo, ls_search_miso
 
        procedure :: ls_search_mimo
        procedure :: ls_search_miso
    end type
 
contains
! ******************************************************************************
! LINE SEARCH MEMBERS
! ------------------------------------------------------------------------------
    pure function ls_get_max_eval(this) result(n)
        class(line_search), intent(in) :: this
        integer(int32) :: n
        n = this%m_maxEval
    end function
 
! --------------------
    subroutine ls_set_max_eval(this, x)
        class(line_search), intent(inout) :: this
        integer(int32), intent(in) :: x
        this%m_maxEval = x
    end subroutine
 
! ------------------------------------------------------------------------------
    pure function ls_get_scale(this) result(x)
        class(line_search), intent(in) :: this
        real(real64) :: x
        x = this%m_alpha
    end function
 
! --------------------
    subroutine ls_set_scale(this, x)
        class(line_search), intent(inout) :: this
        real(real64), intent(in) :: x
        this%m_alpha = x
    end subroutine
 
! ------------------------------------------------------------------------------
    pure function ls_get_dist(this) result(x)
        class(line_search), intent(in) :: this
        real(real64) :: x
        x = this%m_factor
    end function
 
! --------------------
    subroutine ls_set_dist(this, x)
        class(line_search), intent(inout) :: this
        real(real64), intent(in) :: x
        if (x <= 0.0d0) then
            this%m_factor = 0.1d0
        else if (x >= 1.0d0) then
            this%m_factor = 0.99d0
        else
            this%m_factor = x
        end if
    end subroutine
 
! ------------------------------------------------------------------------------
    subroutine ls_search_mimo(this, fcn, xold, grad, dir, x, fvec, fold, fx, &
            ib, err)
        ! Arguments
        class(line_search), intent(in) :: this
        class(vecfcn_helper), intent(in) :: fcn
        real(real64), intent(in), dimension(:) :: xold, dir, grad
        real(real64), intent(out), dimension(:) :: x, fvec
        real(real64), intent(in), optional :: fold
        real(real64), intent(out), optional :: fx
        type(iteration_behavior), optional :: ib
        class(errors), intent(inout), optional, target :: err
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
        real(real64), parameter :: p5 = 0.5d0
        real(real64), parameter :: one = 1.0d0
        real(real64), parameter :: two = 2.0d0
        real(real64), parameter :: three = 3.0d0
 
        ! Local Variables
        logical :: xcnvrg, fcnvrg
        integer(int32) :: i, m, n, neval, niter, flag, maxeval
        real(real64) :: alam, alam1, alamin, f1, slope, temp, test, tmplam, &
            alpha, tolx, lambdamin, f, fo
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
        character(len = 128) :: errmsg
 
        ! Initialization
        xcnvrg = .false.
        fcnvrg = .false.
        neval = 0
        niter = 0
        m = fcn%get_equation_count()
        n = fcn%get_variable_count()
        tolx = two * epsilon(tolx)
        alpha = this%m_alpha
        lambdamin = this%m_factor
        maxeval = this%m_maxEval
        if (present(fx)) fx = zero
        if (present(ib)) then
            ib%iter_count = niter
            ib%fcn_count = neval
            ib%converge_on_fcn = fcnvrg
            ib%converge_on_chng = xcnvrg
            ib%converge_on_zero_diff = .false.
        end if
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Input Checking
        if (.not.fcn%is_fcn_defined()) then
            ! ERROR: No function is defined
            call errmgr%report_error("ls_search_mimo", &
                "No function has been defined.", &
                nl_invalid_operation_error)
            return
        end if
        flag = 0
        if (size(xold) /= n) then
            flag = 3
        else if (size(grad) /= n) then
            flag = 4
        else if (size(dir) /= n) then
            flag = 5
        else if (size(x) /= n) then
            flag = 6
        else if (size(fvec) /= m) then
            flag = 7
        end if
        if (flag /= 0) then
            ! One of the input arrays is not sized correctly
            write(errmsg, 100) "Input number ", flag, &
                " is not sized correctly."
            call errmgr%report_error("ls_search_mimo", trim(errmsg), &
                nl_array_size_error)
            return
        end if
 
        ! Compute 1/2 F * F (* = dot product) if not provided
        if (present(fold)) then
            fo = fold
        else
            ! Evaluate the function, and compute the dot product
            call fcn%fcn(xold, fvec)
            fo = p5 * dot_product(fvec, fvec)
            neval = neval + 1
        end if
 
        ! Compute the slope parameter
        slope = dot_product(grad, dir)
        if (slope >= zero) then
            ! ERROR: The slope should not be pointing uphill - invalid direction
            call errmgr%report_error("ls_search_mimo", &
                "The search direction vector appears to be pointing in " // &
                "an uphill direction -  away from a minimum.", &
                nl_divergent_behavior_error)
            return
        end if
 
        ! Compute the minimum lambda value (length along the search direction)
        test = zero
        do i = 1, n
            temp = abs(dir(i)) / max(abs(xold(i)), one)
            if (temp > test) test = temp
        end do
        alamin = tolx / test
        alam = one
 
        ! Iteration Loop
        flag = 0 ! Used to check for convergence errors
        do
            ! Step along the specified direction by the amount ALAM
            x = xold + alam * dir
            call fcn%fcn(x, fvec)
            f = p5 * dot_product(fvec, fvec)
            neval = neval + 1
            niter = niter + 1
 
            ! Check the step
            if (alam < alamin) then
                ! Either the solution has converged due to negligible change in
                ! the root values, or the line search may have fully
                ! backtracked.  In the event the solution fully backtracked,
                ! we'll issue a warning to inform the user of the potential
                ! issue.
                if (norm2(x - xold) == zero) then
                    ! The line search fully backtracked
                    call errmgr%report_warning("ls_search_mimo", &
                        "The line search appears to have fully " // &
                        "backtracked.  As such, check results carefully, " // &
                        "and/or consider attempting the solve without " // &
                        "the line search.", &
                        nl_convergence_error)
                end if
                x = xold
                xcnvrg = .true.
                exit
            else if (f <= fo + alpha * alam * slope) then
                ! The function has converged
                fcnvrg = .true.
                exit
            else
                ! Convergence has not yet occurred, continue backtracking
                tmplam = min_backtrack_search(niter, fo, f, f1, alam, &
                    alam1, slope)
            end if
 
            ! Set up parameters for the cubic model as we've already been
            ! through once with the quadratic model without success.
            alam1 = alam
            f1 = f
            alam = max(tmplam, lambdamin * alam)
 
            ! Ensure we haven't performed too many function evaluations
            if (neval >= maxeval) then
                ! ERROR: Too many function evaluations
                flag = 1
                exit
            end if
        end do
        if (present(fx)) fx = f
 
        ! Report out iteration statistics
        if (present(ib)) then
            ib%iter_count = niter
            ib%fcn_count = neval
            ib%converge_on_fcn = fcnvrg
            ib%converge_on_chng = xcnvrg
            ib%converge_on_zero_diff = .false.
        end if
 
        ! Check for convergence issues
        if (flag /= 0) then
            write(errmsg, 100) "The line search failed to " // &
                "converge.  Function evaluations performed: ", neval, "."
            call errmgr%report_error("ls_search_mimo", trim(errmsg), &
                nl_convergence_error)
        end if
 
        ! Formatting
100     format(a, i0, a)
    end subroutine
 
! ------------------------------------------------------------------------------
    subroutine ls_search_miso(this, fcn, xold, grad, dir, x, fold, fx, &
            ib, err)
        ! Arguments
        class(line_search), intent(in) :: this
        class(fcnnvar_helper), intent(in) :: fcn
        real(real64), intent(in), dimension(:) :: xold, dir, grad
        real(real64), intent(out), dimension(:) :: x
        real(real64), intent(in), optional :: fold
        real(real64), intent(out), optional :: fx
        type(iteration_behavior), optional :: ib
        class(errors), intent(inout), optional, target :: err
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
        real(real64), parameter :: p5 = 0.5d0
        real(real64), parameter :: one = 1.0d0
        real(real64), parameter :: two = 2.0d0
        real(real64), parameter :: three = 3.0d0
 
        ! Local Variables
        logical :: xcnvrg, fcnvrg
        integer(int32) :: i, n, neval, niter, flag, maxeval
        real(real64) :: alam, alam1, alamin, f1, slope, temp, test, tmplam, &
            alpha, tolx, lambdamin, f, fo
        class(errors), pointer :: errmgr
        type(errors), target :: deferr
        character(len = 128) :: errmsg
 
        ! Initialization
        xcnvrg = .false.
        fcnvrg = .false.
        neval = 0
        niter = 0
        n = fcn%get_variable_count()
        tolx = two * epsilon(tolx)
        alpha = this%m_alpha
        lambdamin = this%m_factor
        maxeval = this%m_maxEval
        if (present(fx)) fx = zero
        if (present(ib)) then
            ib%iter_count = niter
            ib%fcn_count = neval
            ib%converge_on_fcn = fcnvrg
            ib%converge_on_chng = xcnvrg
            ib%converge_on_zero_diff = .false.
        end if
        if (present(err)) then
            errmgr => err
        else
            errmgr => deferr
        end if
 
        ! Input Checking
        if (.not.fcn%is_fcn_defined()) then
            ! ERROR: No function is defined
            call errmgr%report_error("ls_search_miso", &
                "No function has been defined.", &
                nl_invalid_operation_error)
            return
        end if
        flag = 0
        if (size(xold) /= n) then
            flag = 3
        else if (size(grad) /= n) then
            flag = 4
        else if (size(dir) /= n) then
            flag = 5
        else if (size(x) /= n) then
            flag = 6
        end if
        if (flag /= 0) then
            ! One of the input arrays is not sized correctly
            write(errmsg, 100) "Input number ", flag, &
                " is not sized correctly."
            call errmgr%report_error("ls_search_miso", trim(errmsg), &
                nl_array_size_error)
            return
        end if
 
        ! Establish the "old" function value
        if (present(fold)) then
            fo = fold
        else
            ! Evaluate the function
            fo = fcn%fcn(xold)
            neval = neval + 1
        end if
 
        ! Compute the slope parameter
        slope = dot_product(grad, dir)
        if (slope >= zero) then
            ! ERROR: The slope should not be pointing uphill - invalid direction
            call errmgr%report_error("ls_search_miso", &
                "The search direction vector appears to be pointing in " // &
                "an uphill direction -  away from a minimum.", &
                nl_divergent_behavior_error)
            return
        end if
 
        ! Compute the minimum lambda value (length along the search direction)
        test = zero
        do i = 1, n
            temp = abs(dir(i)) / max(abs(xold(i)), one)
            if (temp > test) test = temp
        end do
        alamin = tolx / test
        alam = one
 
        ! Iteration Loop
        flag = 0 ! Used to check for convergence errors
        do
            ! Step along the specified direction by the amount ALAM
            x = xold + alam * dir
            f = fcn%fcn(x)
            neval = neval + 1
            niter = niter + 1
 
            ! Check the step
            if (alam < alamin) then
                ! Either the solution has converged due to negligible change in
                ! the root values, or the line search may have fully
                ! backtracked.  In the event the solution fully backtracked,
                ! we'll issue a warning to inform the user of the potential
                ! issue.
                if (norm2(x - xold) == zero) then
                    ! The line search fully backtracked
                    call errmgr%report_warning("ls_search_miso", &
                        "The line search appears to have fully " // &
                        "backtracked.  As such, check results carefully, " // &
                        "and/or consider attempting the solve without " // &
                        "the line search.", &
                        nl_convergence_error)
                end if
                x = xold
                xcnvrg = .true.
                exit
            else if (f <= fo + alpha * alam * slope) then
                ! The function has converged
                fcnvrg = .true.
                exit
            else
                ! Convergence has not yet occurred, continue backtracking
                tmplam = min_backtrack_search(niter, fo, f, f1, alam, &
                    alam1, slope)
            end if
 
            ! Set up parameters for the cubic model as we've already been
            ! through once with the quadratic model without success.
            alam1 = alam
            f1 = f
            alam = max(tmplam, lambdamin * alam)
 
            ! Ensure we haven't performed too many function evaluations
            if (neval >= maxeval) then
                ! ERROR: Too many function evaluations
                flag = 1
                exit
            end if
        end do
        if (present(fx)) fx = f
 
        ! Report out iteration statistics
        if (present(ib)) then
            ib%iter_count = niter
            ib%fcn_count = neval
            ib%converge_on_fcn = fcnvrg
            ib%converge_on_chng = xcnvrg
            ib%converge_on_zero_diff = .false.
        end if
 
        ! Check for convergence issues
        if (flag /= 0) then
            write(errmsg, 100) "The line search failed to " // &
                "converge.  Function evaluations performed: ", neval, "."
            call errmgr%report_error("ls_search_miso", trim(errmsg), &
                nl_convergence_error)
        end if
 
        ! Formatting
100     format(a, i0, a)
    end subroutine
 
! ------------------------------------------------------------------------------
    pure function min_backtrack_search(mode, f0, f, f1, alam, alam1, &
            slope) result(lam)
        ! Arguments
        integer(int32), intent(in) :: mode
        real(real64), intent(in) :: f0, f, f1, alam, alam1, slope
        real(real64) :: lam
 
        ! Parameters
        real(real64), parameter :: zero = 0.0d0
        real(real64), parameter :: p5 = 0.5d0
        real(real64), parameter :: two = 2.0d0
        real(real64), parameter :: three = 3.0d0
 
        ! Local Variables
        real(real64) :: rhs1, rhs2, a, b, disc
 
        ! Process
        if (mode == 1) then
            ! Use a quadratic function approximation
            lam = -slope / (two * (f - f0 - slope))
        else
            ! Use a cubic function
            rhs1 = f - f0 - alam * slope
            rhs2 = f1 - f0 - alam1 * slope
            a = (rhs1 / alam**2 - rhs2 / alam1**2) / (alam - alam1)
            b = (-alam1 * rhs1 / alam**2 + alam * rhs2 / alam1**2) / &
                (alam - alam1)
            if (a == zero) then
                lam = -slope / (two * b)
            else
                disc = b**2 - three * a * slope
                if (disc < zero) then
                    lam = p5 * alam
                else if (b <= zero) then
                    lam = (-b + sqrt(disc)) / (three * a)
                else
                    lam = -slope / (b + sqrt(disc))
                end if
            end if
            if (lam > p5 * alam) lam = p5 * alam
        end if
    end function
 
! ------------------------------------------------------------------------------
    subroutine limit_search_vector(x, lim)
        ! Arguments
        real(real64), intent(inout), dimension(:) :: x
        real(real64), intent(in) :: lim
 
        ! Local Variables
        real(real64) :: mag
 
        ! Process
        mag = norm2(x)
        if (mag == 0.0d0) return
        if (mag > lim) x = (lim / mag) * x
    end subroutine
 
! ------------------------------------------------------------------------------
end module