initialize_nlesolver_variables Subroutine

private subroutine initialize_nlesolver_variables(me, n, m, max_iter, tol, alpha, alpha_min, alpha_max, tolx, fmin_tol, backtrack_c, backtrack_tau, use_broyden, broyden_update_n, step_mode, func, grad, grad_sparse, export_iteration, user_input_check, verbose, iunit, n_uphill_max, n_intervals, sparsity_mode, irow, icol, atol, btol, conlim, damp, itnlim, nout, lusol_method, localSize, custom_solver_sparse)

Constructor for the class.

Type Bound

nlesolver_type

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(inout) :: me
integer, intent(in) :: n

number of optimization variables

integer, intent(in) :: m

number of constraints

integer, intent(in) :: max_iter

maximum number of iterations

real(kind=wp), intent(in) :: tol

function convergence tolerance

real(kind=wp), intent(in), optional :: alpha

constant step length for step_mode=1 (0,1]

real(kind=wp), intent(in), optional :: alpha_min

minimum step length (0,1]

real(kind=wp), intent(in), optional :: alpha_max

maximum step length (0,1]

real(kind=wp), intent(in), optional :: tolx

convergence tolerance for changes in x

real(kind=wp), intent(in), optional :: fmin_tol

convergence tolerance for fmin (used when step_mode=3)

real(kind=wp), intent(in), optional :: backtrack_c

backtracking linesearch parameter (0,1)

real(kind=wp), intent(in), optional :: backtrack_tau

backtracking linesearch parameter (0,1)

logical, intent(in), optional :: use_broyden

use a Broyden update (default is False)

integer, intent(in), optional :: broyden_update_n

For Broyden mode, update the full Jacobian at most every this many iterations (must be >1) If <=1 then Jacobian is only computed on the first iteration.

integer, intent(in), optional :: step_mode

step mode:

  • 1 = use the specified alpha (0,1]
  • 2 = backtracking linesearch (between alpha_min and alpha_max)
  • 3 = exact linesearch (between alpha_min and alpha_max)
  • 4 = evaluate function at specified fixed points (between alpha_min and alpha_max)
procedure(func_func) :: func

computes the function vector

procedure(grad_func), optional :: grad

computes the jacobian [required for dense mode: sparsity_mode=1]

procedure(grad_func_sparse), optional :: grad_sparse

computes the jacobian [required for sparse mode: sparsity_mode>1]

procedure(export_func), optional :: export_iteration

function to export each iteration

procedure(wait_func), optional :: user_input_check

check for user input (to stop solver if necessary)

logical, intent(in), optional :: verbose

for verbose status printing

integer, intent(in), optional :: iunit

unit for verbose printing (assumed to be open). by default this is output_unit.

integer, intent(in), optional :: n_uphill_max

maximum number of consecutive steps to allow where the value of f increases

integer, intent(in), optional :: n_intervals

number of intervals for fixed point linesearch

integer, intent(in), optional :: sparsity_mode

sparsity mode:

  • 1 - assume dense (use dense solver) Must specify grad for this mode.
  • 2 - assume sparse (use LSQR sparse solver).
  • 3 - assume sparse (use LUSOL sparse solver).
  • 4 - assume sparse (use LSMR sparse solver).
  • 5 - assume sparse (use the user provided sparse solver custom_solver_sparse). Must also specify grad_sparse and the irow and icol sparsity pattern for mode>1.
integer, intent(in), optional, dimension(:) :: irow

sparsity pattern nonzero elements row indices. must be specified with icol and be the same length (n_nonzeros).

integer, intent(in), optional, dimension(:) :: icol

sparsity pattern nonzero elements column indices must be specified with icol and be the same length (n_nonzeros).

real(kind=wp), intent(in), optional :: atol

LSQR: relative error in definition of A

real(kind=wp), intent(in), optional :: btol

LSQR: relative error in definition of b

real(kind=wp), intent(in), optional :: conlim

LSQR: An upper limit on cond(Abar), the apparent condition number of the matrix Abar.

real(kind=wp), intent(in), optional :: damp

LSQR: damp factor

integer, intent(in), optional :: itnlim

LSQR: max iterations

integer, intent(in), optional :: nout

LSQR: output unit for printing

integer, intent(in), optional :: lusol_method

LUSOL method:

  • 0 => TPP: Threshold Partial Pivoting.
  • 1 => TRP: Threshold Rook Pivoting.
  • 2 => TCP: Threshold Complete Pivoting.
integer, intent(in), optional :: localSize

LSMR: Number of vectors for local reorthogonalization:

  • 0 No reorthogonalization is performed.
  • 0 This many n-vectors "v" (the most recent ones) are saved for reorthogonalizing the next v.

localSize need not be more than min(m,n). At most min(m,n) vectors will be allocated.

procedure(sparse_solver_func), optional :: custom_solver_sparse

for sparsity_mode=5, this is the user-provided linear solver.


Calls

proc~~initialize_nlesolver_variables~~CallsGraph proc~initialize_nlesolver_variables nlesolver_module::nlesolver_type%initialize_nlesolver_variables proc~set_status nlesolver_module::nlesolver_type%set_status proc~initialize_nlesolver_variables->proc~set_status

Source Code

    subroutine initialize_nlesolver_variables(me,&
                    n,m,max_iter,tol,alpha,alpha_min,alpha_max,tolx,fmin_tol,&
                    backtrack_c,backtrack_tau,&
                    use_broyden,broyden_update_n,step_mode,&
                    func,grad,grad_sparse,&
                    export_iteration,user_input_check,&
                    verbose,iunit,n_uphill_max,n_intervals,&
                    sparsity_mode,irow,icol,&
                    atol,btol,conlim,damp,itnlim,nout,&
                    lusol_method,localSize,custom_solver_sparse )

    implicit none

    class(nlesolver_type),intent(inout) :: me
    integer,intent(in)                :: n                  !! number of optimization variables
    integer,intent(in)                :: m                  !! number of constraints
    integer,intent(in)                :: max_iter           !! maximum number of iterations
    real(wp),intent(in)               :: tol                !! function convergence tolerance
    procedure(func_func)              :: func               !! computes the function vector
    procedure(grad_func),optional        :: grad            !! computes the jacobian [required for dense mode: `sparsity_mode=1`]
    procedure(grad_func_sparse),optional :: grad_sparse     !! computes the jacobian [required for sparse mode: `sparsity_mode>1`]
    integer,intent(in),optional       :: step_mode          !! step mode:
                                                            !!
                                                            !!  * 1 = use the specified `alpha` (0,1]
                                                            !!  * 2 = backtracking linesearch (between `alpha_min` and `alpha_max`)
                                                            !!  * 3 = exact linesearch (between `alpha_min` and `alpha_max`)
                                                            !!  * 4 = evaluate function at specified fixed points  (between `alpha_min` and `alpha_max`)
    real(wp),intent(in),optional      :: alpha              !! constant step length for `step_mode=1` (0,1]
    real(wp),intent(in),optional      :: alpha_min          !! minimum step length (0,1]
    real(wp),intent(in),optional      :: alpha_max          !! maximum step length (0,1]
    real(wp),intent(in),optional      :: tolx               !! convergence tolerance for changes in `x`
    real(wp),intent(in),optional      :: fmin_tol           !! convergence tolerance for [[fmin]] (used when `step_mode=3`)
    real(wp),intent(in),optional      :: backtrack_c        !! backtracking linesearch parameter (0,1)
    real(wp),intent(in),optional      :: backtrack_tau      !! backtracking linesearch parameter (0,1)
    logical,intent(in),optional       :: use_broyden        !! use a Broyden update (default is False)
    integer,intent(in),optional       :: broyden_update_n   !! For Broyden mode, update the full Jacobian
                                                            !! at most every this many iterations (must be >1)
                                                            !! If <=1 then Jacobian is only computed on the
                                                            !! first iteration.
    procedure(export_func),optional   :: export_iteration   !! function to export each iteration
    procedure(wait_func),optional     :: user_input_check   !! check for user input (to stop solver if necessary)
    logical,intent(in),optional       :: verbose            !! for verbose status printing
    integer,intent(in),optional       :: iunit              !! unit for verbose printing (assumed to be open).
                                                            !! by default this is `output_unit`.
    integer,intent(in),optional       :: n_uphill_max       !! maximum number of consecutive steps
                                                            !! to allow where the value of `f` increases
    integer,intent(in),optional       :: n_intervals        !! number of intervals for fixed point linesearch

    integer,intent(in),optional :: sparsity_mode !! sparsity mode:
                                                 !!
                                                 !! * 1 - assume dense (use dense solver)
                                                 !!   Must specify `grad` for this mode.
                                                 !! * 2 - assume sparse (use LSQR sparse solver).
                                                 !! * 3 - assume sparse (use LUSOL sparse solver).
                                                 !! * 4 - assume sparse (use LSMR sparse solver).
                                                 !! * 5 - assume sparse (use the user provided sparse
                                                 !!   solver `custom_solver_sparse`).
                                                 !! Must also specify `grad_sparse` and the `irow` and `icol`
                                                 !! sparsity pattern for `mode>1`.
    integer,dimension(:),intent(in),optional :: irow !! sparsity pattern nonzero elements row indices.
                                                     !! must be specified with `icol` and be the same length (`n_nonzeros`).
    integer,dimension(:),intent(in),optional :: icol !! sparsity pattern nonzero elements column indices
                                                     !! must be specified with `icol` and be the same length (`n_nonzeros`).
    real(wp),intent(in),optional      :: atol    !! `LSQR`: relative error in definition of `A`
    real(wp),intent(in),optional      :: btol    !! `LSQR`: relative error in definition of `b`
    real(wp),intent(in),optional      :: conlim  !! `LSQR`: An upper limit on `cond(Abar)`, the apparent
                                                 !! condition number of the matrix `Abar`.
    real(wp),intent(in),optional      :: damp    !! `LSQR`: damp factor
    integer,intent(in),optional       :: itnlim  !! `LSQR`: max iterations
    integer,intent(in),optional       :: nout    !! `LSQR`: output unit for printing
    integer,intent(in),optional       :: lusol_method !! `LUSOL` method:
                                                      !!
                                                      !! * 0 => TPP: Threshold Partial   Pivoting.
                                                      !! * 1 => TRP: Threshold Rook      Pivoting.
                                                      !! * 2 => TCP: Threshold Complete  Pivoting.
    integer,intent(in),optional :: localSize !! `LSMR`: Number of vectors for local reorthogonalization:
                                             !!
                                             !!  *  0    No reorthogonalization is performed.
                                             !!  * >0    This many n-vectors "v" (the most recent ones)
                                             !!          are saved for reorthogonalizing the next v.
                                             !!
                                             !! localSize need not be more than `min(m,n)`.
                                             !! At most `min(m,n)` vectors will be allocated.
    procedure(sparse_solver_func),optional :: custom_solver_sparse !! for `sparsity_mode=5`, this is the
                                                                   !! user-provided linear solver.

    logical :: status_ok !! true if there were no errors

    status_ok = .true.

    !required:
    me%n        = abs(n)
    me%m        = abs(m)
    me%max_iter = abs(max_iter)
    me%tol      = abs(tol)
    me%func     => func
    if (present(grad)) me%grad               => grad
    if (present(grad_sparse)) me%grad_sparse => grad_sparse

    !optional:

    if (present(step_mode)) then
        select case (step_mode)
        case(1) ! = use the specified `alpha` (0,1]
            me%linesearch => simple_step
        case(2) ! = backtracking linesearch (between `alpha_min` and `alpha_max`)
            me%linesearch => backtracking_linesearch
        case(3) ! = exact linesearch (between `alpha_min` and `alpha_max`)
            me%linesearch => exact_linesearch
        case(4) ! = evaluate function at specified fixed points (between `alpha_min` and `alpha_max`)
            me%linesearch => fixed_point_linesearch
        case default
            status_ok = .false.
            call me%set_status(istat = -5, string = 'Error: invalid step_mode:',i=step_mode)
            return
        end select
    else
        me%linesearch => simple_step
    end if

    if (present(alpha))             me%alpha             = abs(alpha)
    if (present(alpha_min))         me%alpha_min         = abs(alpha_min)
    if (present(alpha_max))         me%alpha_max         = abs(alpha_max)
    if (present(tolx))              me%tolx              = abs(tolx)
    if (present(backtrack_c))       me%c                 = abs(backtrack_c)
    if (present(backtrack_tau))     me%tau               = abs(backtrack_tau)
    if (present(use_broyden))       me%use_broyden       = use_broyden
    if (present(broyden_update_n))  me%broyden_update_n  = abs(broyden_update_n)
    if (present(verbose))           me%verbose           = verbose
    if (present(iunit))             me%iunit             = iunit
    if (present(n_uphill_max))      me%n_uphill_max      = abs(n_uphill_max)
    if (present(n_intervals))       me%n_intervals       = max(abs(n_intervals),1)
    if (present(fmin_tol))          me%fmin_tol          = abs(fmin_tol)

    if (present(export_iteration)) me%export_iteration  => export_iteration
    if (present(user_input_check)) me%user_input_check  => user_input_check

    ! error checks:
    if (me%alpha<zero .or. me%alpha>one) then
        status_ok = .false.
        call me%set_status(istat = -1, string = 'Error: invalid alpha:',r=me%alpha)
        return
    end if
    if (me%alpha_min<zero .or. me%alpha_min>one) then
        status_ok = .false.
        call me%set_status(istat = -2, string = 'Error: invalid alpha_min:',r=me%alpha_min)
        return
    end if
    if (me%alpha_max<zero .or. me%alpha_max>one) then
        status_ok = .false.
        call me%set_status(istat = -3, string = 'Error: invalid alpha_max:',r=me%alpha_max)
        return
    end if
    if (me%alpha_max<=me%alpha_min) then
        status_ok = .false.
        call me%set_status(istat = -4, string = 'Error: alpha_min must be < alpha_max')
        return
    end if
    if (me%c<zero .or. me%c>one) then
        status_ok = .false.
        call me%set_status(istat = -12, string = 'Error: backtracking linesearch c must be in range (0, 1):',r=me%c)
        return
    end if
    if (me%tau<zero .or. me%tau>one) then
        status_ok = .false.
        call me%set_status(istat = -13, string = 'Error: backtracking linesearch tau must be in range (0, 1):',r=me%tau)
        return
    end if

    ! initialize:
    me%n_nonzeros = -1 ! not used
    if (allocated(me%irow)) deallocate(me%irow)
    if (allocated(me%icol)) deallocate(me%icol)

    if (present(custom_solver_sparse)) me%custom_solver_sparse => custom_solver_sparse

    if (present(sparsity_mode)) then
        me%sparsity_mode = sparsity_mode
        if (sparsity_mode>NLESOLVER_SPARSITY_DENSE) then ! sparse solver method
            if (present(irow) .and. present(icol) .and. present(grad_sparse)) then
                if (size(irow) == size(icol)) then
                    me%n_nonzeros = size(irow)
                    me%irow = irow
                    me%icol = icol
                else
                    call me%set_status(istat = -15, string = 'Error: irow and icol must be the same length')
                    return
                end if
            else
                call me%set_status(istat = -14, string = 'Error: must specify grad_sparse, irow, and icol for sparsity_mode > 1')
                return
            end if
            ! LSQR optional inputs:
            if (present(atol))   me%atol   = atol
            if (present(btol))   me%btol   = btol
            if (present(conlim)) me%conlim = conlim
            if (present(damp))   me%damp   = damp
            if (present(itnlim)) me%itnlim = itnlim
            if (present(nout))   me%nout   = nout

            ! LUSOL method
            if (present(nout)) me%lusol_method = lusol_method

            ! LSMR optional inputs:
            if (present(localSize)) me%localSize = localSize

        end if
        if (sparsity_mode==NLESOLVER_SPARSITY_CUSTOM_SPARSE) then
            if (.not. associated(me%custom_solver_sparse)) then
                call me%set_status(istat = -16, string = 'Error: must specify custom_solver_sparse for sparsity_mode = 5')
                return
            end if
        end if
    end if

    if (status_ok) then
        call me%set_status(istat = 0, string = 'Class successfully initialized')
    end if

    end subroutine initialize_nlesolver_variables