nlesolver_type – nlesolver-fortran

type, public :: nlesolver_type

Nonlinear equations solver class.

Components

Type	Visibility	Attributes		Name		Initial
integer,	private		::	n	=	0	number of opt vars
integer,	private		::	m	=	0	number of constraints
integer,	private		::	max_iter	=	100	maximum number of iterations
real(kind=wp),	private		::	tol	=	1.0e-6_wp	convergence tolerance for function values
real(kind=wp),	private		::	alpha	=	1.0_wp	step length (when specified constant)
real(kind=wp),	private		::	alpha_min	=	0.1_wp	minimum step length (when allowed to vary)
real(kind=wp),	private		::	alpha_max	=	1.0_wp	maximum step length (when allowed to vary)
real(kind=wp),	private		::	tolx	=	1.0e-8_wp	convergence tolerance for `x`
real(kind=wp),	private		::	c	=	0.5_wp	backtracking linesearch parameter (0,1)
real(kind=wp),	private		::	tau	=	0.5_wp	backtracking linesearch parameter (0,1)
real(kind=wp),	private		::	fmin_tol	=	1.0e-5_wp	tolerance for "exact" linesearch
integer,	private		::	n_intervals	=	2	number of intervals for fixed point linesearch
logical,	private		::	use_broyden	=	.false.	if true, a Broyden update is used rather than computing the Jacobian at every step. The `grad` function is only called for the initial evaluation.
integer,	private		::	broyden_update_n	=	4	if this value is `>0`, the Broyden update is computed at most this many times before the full Jacobian is recomputed.
integer,	private		::	n_uphill_max	=	5	maximum number of consecutive steps to allow where the value of `f` increases
logical,	private		::	verbose	=	.false.	verbose output printing
integer,	private		::	iunit	=	output_unit	output unit for printing (assumed to be open).
character(len=:),	private,	allocatable	::	message			latest status message
integer,	private		::	istat	=	-999	latest status message
integer,	private		::	norm_mode	=	NLESOLVER_2_NORM	how to compute the norm of the function vector: 1 = 2-norm (default) 2 = infinity-norm 3 = 1-norm
integer,	private		::	bounds_mode	=	NLESOLVER_IGNORE_BOUNDS	how to handle the `x` variable bounds: 0 = ignore bounds (default). 1 = scalar mode : adjust the step vector (`xnew-x`) by moving each individual scalar component of `xnew` to be within the bounds. Note that this can change the direction and magnitude of the search vector. 2 = vector mode: backtrack the search vector so that no bounds are violated. note that this does not change the direction of the vector, only the magnitude. 3 = wall mode: adjust the step vector (`xnew-x`) by moving each individual scalar component of `xnew` to be within the bounds. And then holding the ones changed fixed during the line search. Note: if `bounds_mode>0`, then the initial `x` vector is also adjusted to be within the bounds before the first step, if necessary. This is just done by adjusting each scalar component of `x` to be within bounds.
real(kind=wp),	private,	dimension(:), allocatable	::	xlow			lower bounds for `x` (size is `n`). only used if `bounds_mode>0`.
real(kind=wp),	private,	dimension(:), allocatable	::	xupp			upper bounds for `x` (size is `n`). only used if `bounds_mode>0`.
procedure(func_func),	private,	pointer	::	func	=>	null()	user-supplied routine to compute the function
procedure(export_func),	private,	pointer	::	export_iteration	=>	null()	user-supplied routine to export iterations
procedure(wait_func),	private,	pointer	::	user_input_check	=>	null()	user-supplied routine to enable user to stop iterations
procedure(linesearch_func),	private,	pointer	::	linesearch	=>	null()	line search method (determined by `step_mode` user input in initialize)
integer,	private		::	sparsity_mode	=	NLESOLVER_SPARSITY_DENSE	sparsity mode: 1 - assume dense (use dense solver). 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 a user provided sparse solver).
integer,	private		::	n_nonzeros	=	-1	number of nonzero Jacobian elements (used for `sparsity_mode > 1`)
integer,	private,	dimension(:), allocatable	::	irow			sparsity pattern nonzero elements row indices.
integer,	private,	dimension(:), allocatable	::	icol			sparsity pattern nonzero elements column indices
real(kind=wp),	private		::	atol	=	zero	relative error in definition of `A`
real(kind=wp),	private		::	btol	=	zero	relative error in definition of `b`
real(kind=wp),	private		::	conlim	=	zero	An upper limit on `cond(Abar)`, the apparent condition number of the matrix `Abar`.
integer,	private		::	itnlim	=	100	max iterations
integer,	private		::	nout	=	0	output unit for printing
real(kind=wp),	private		::	damp	=	zero	damp parameter for LSQR
integer,	private		::	localSize	=	0	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.
integer,	private		::	lusol_method	=	0	0 => TPP: Threshold Partial Pivoting. 1 => TRP: Threshold Rook Pivoting. 2 => TCP: Threshold Complete Pivoting.
procedure(grad_func),	private,	pointer	::	grad	=>	null()	user-supplied routine to compute the gradient of the function (dense version)
procedure(grad_func_sparse),	private,	pointer	::	grad_sparse	=>	null()	user-supplied routine to compute the gradient of the function (sparse version)
procedure(sparse_solver_func),	private,	pointer	::	custom_solver_sparse	=>	null()	user-supplied sparse linear solver routine (used for `sparsity_mode=5`)
procedure(norm_func),	private,	pointer	::	norm	=>	null()	function for computing the norm of the `f` vector

Type-Bound Procedures

procedure, public :: initialize => initialize_nlesolver_variables

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, bounds_mode, xlow, xupp, norm_mode)

Constructor for the class.

Arguments

Type	Intent	Optional	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: Read more…
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: Read more…
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: Read more…
integer,	intent(in),	optional		::	localSize	`LSMR`: Number of vectors for local reorthogonalization: Read more…
procedure(sparse_solver_func),		optional		::	custom_solver_sparse	for `sparsity_mode=5`, this is the user-provided linear solver.
integer,	intent(in),	optional		::	bounds_mode	how to handle the `x` variable bounds: Read more…
real(kind=wp),	intent(in),	optional,	dimension(n)	::	xlow	lower bounds for `x` (size is `n`). only used if `bounds_mode>0` and both `xlow` and `xupp` are specified.
real(kind=wp),	intent(in),	optional,	dimension(n)	::	xupp	upper bounds for `x` (size is `n`). only used if `bounds_mode>0` and both `xlow` and `xupp` are specified.
integer,	intent(in),	optional		::	norm_mode	how to compute the norm of the function vector: Read more…

procedure, public :: solve => nlesolver_solver

private subroutine nlesolver_solver(me, x)

Main solver.

Arguments

Type	Intent	Optional	Attributes		Name
class(nlesolver_type),	intent(inout)			::	me
real(kind=wp),	intent(inout),		dimension(:)	::	x

procedure, public :: destroy => destroy_nlesolver_variables

private subroutine destroy_nlesolver_variables(me)

Destructor

Arguments

Type	Intent	Optional	Attributes		Name
class(nlesolver_type),	intent(out)			::	me

procedure, public :: status => get_status

private subroutine get_status(me, istat, message)

Return the status code and message from the nlesolver_type class.

Arguments

Type	Intent	Optional	Attributes		Name
class(nlesolver_type),	intent(inout)			::	me
integer,	intent(out),	optional		::	istat	Integer status code.
character(len=:),	intent(out),	optional,	allocatable	::	message	Text status message

procedure, private :: set_status

private subroutine set_status(me, istat, string, i, r)

Set status flag and message.

Arguments

Type	Intent	Optional		Name
class(nlesolver_type),	intent(inout)		::	me
integer,	intent(in)		::	istat	status code
character(len=*),	intent(in)		::	string	status message
integer,	intent(in),	optional	::	i	an integer value to append
real(kind=wp),	intent(in),	optional	::	r	a real value to append

procedure, private :: adjust_x_for_bounds

private subroutine adjust_x_for_bounds(me, x)

if necessary, adjust the x vector to be within the bounds. used for the initial guess.

Arguments

Type	Intent	Optional	Attributes		Name
class(nlesolver_type),	intent(inout)			::	me
real(kind=wp),	intent(inout),		dimension(me%n)	::	x	the `x` vector to adjust

procedure, private :: adjust_search_direction

private subroutine adjust_search_direction(me, x, p, pnew, modified)

if necessary, adjust the search direction vector p so that bounds are not violated.

Arguments

Type	Intent	Attributes		Name
class(nlesolver_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(me%n)	::	x	initial `x`. it is assumed that this does not violate any bounds.
real(kind=wp),	intent(in),	dimension(me%n)	::	p	search direction `p = xnew - x`
real(kind=wp),	intent(out),	dimension(me%n)	::	pnew	new search direction
logical,	intent(out),	dimension(me%n)	::	modified	indicates the elements of `p` that were modified

procedure, private :: compute_next_step

private subroutine compute_next_step(me, xold, search_direction, alpha, modified, xnew)

Compute the next step.

Arguments

Type	Intent	Attributes		Name
class(nlesolver_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(me%n)	::	xold	initial `x`
real(kind=wp),	intent(in),	dimension(me%n)	::	search_direction	search direction vector
real(kind=wp),	intent(in)		::	alpha	step length to take in the search direction
logical,	intent(in),	dimension(me%n)	::	modified	indicates the elements of `p` that were modified because they violated the bounds. Output of adjust_search_direction.
real(kind=wp),	intent(out),	dimension(me%n)	::	xnew	final `x`

Source Code

        type,public :: nlesolver_type

        !! Nonlinear equations solver class.

        private

        integer     :: n                = 0             !! number of opt vars
        integer     :: m                = 0             !! number of constraints
        integer     :: max_iter         = 100           !! maximum number of iterations
        real(wp)    :: tol              = 1.0e-6_wp     !! convergence tolerance for function values
        real(wp)    :: alpha            = 1.0_wp        !! step length (when specified constant)
        real(wp)    :: alpha_min        = 0.1_wp        !! minimum step length (when allowed to vary)
        real(wp)    :: alpha_max        = 1.0_wp        !! maximum step length (when allowed to vary)
        real(wp)    :: tolx             = 1.0e-8_wp     !! convergence tolerance for `x`
        real(wp)    :: c                = 0.5_wp        !! backtracking linesearch parameter (0,1)
        real(wp)    :: tau              = 0.5_wp        !! backtracking linesearch parameter (0,1)
        real(wp)    :: fmin_tol         = 1.0e-5_wp     !! tolerance for "exact" linesearch
        integer     :: n_intervals      = 2             !! number of intervals for fixed point linesearch
        logical     :: use_broyden      = .false.       !! if true, a Broyden update is used
                                                        !! rather than computing the Jacobian
                                                        !! at every step. The `grad` function is
                                                        !! only called for the initial evaluation.
        integer     :: broyden_update_n = 4             !! if this value is `>0`, the Broyden update
                                                        !! is computed at most this many times before
                                                        !! the full Jacobian is recomputed.
        integer     :: n_uphill_max     = 5             !! maximum number of consecutive steps
                                                        !! to allow where the value of `f` increases
        logical     :: verbose          = .false.       !! verbose output printing
        integer     :: iunit            = output_unit   !! output unit for printing (assumed to be open).
        character(len=:),allocatable :: message         !! latest status message
        integer     :: istat            = -999          !! latest status message

        integer :: norm_mode = NLESOLVER_2_NORM  !! how to compute the norm of the function vector:
                                                 !!
                                                 !! * 1 = 2-norm (default)
                                                 !! * 2 = infinity-norm
                                                 !! * 3 = 1-norm

        integer :: bounds_mode = NLESOLVER_IGNORE_BOUNDS  !! how to handle the `x` variable bounds:
                                    !!
                                    !!  * 0 = ignore bounds (default).
                                    !!  * 1 = scalar mode : adjust the step vector (`xnew-x`) by moving each individual scalar component
                                    !!    of `xnew` to be within the bounds. Note that this can change the direction and magnitude of
                                    !!    the search vector.
                                    !!  * 2 = vector mode: backtrack the search vector so that no bounds are violated. note that
                                    !!    this does not change the direction of the vector, only the magnitude.
                                    !!  * 3 = wall mode: adjust the step vector (`xnew-x`) by moving each individual scalar component
                                    !!    of `xnew` to be within the bounds. And then holding the ones changed fixed during the line search.
                                    !!
                                    !! Note: if `bounds_mode>0`, then the initial `x` vector is also adjusted
                                    !! to be within the bounds before the first step, if necessary. This is just done
                                    !! by adjusting each scalar component of `x` to be within bounds.
        real(wp),dimension(:),allocatable :: xlow !! lower bounds for `x` (size is `n`). only used if `bounds_mode>0`.
        real(wp),dimension(:),allocatable :: xupp !! upper bounds for `x` (size is `n`). only used if `bounds_mode>0`.

        procedure(func_func),pointer       :: func             => null() !! user-supplied routine to compute the function
        procedure(export_func),pointer     :: export_iteration => null() !! user-supplied routine to export iterations
        procedure(wait_func),pointer       :: user_input_check => null() !! user-supplied routine to enable user to stop iterations
        procedure(linesearch_func),pointer :: linesearch       => null() !! line search method (determined by `step_mode` user input in [[nlesolver_type:initialize]])

        ! sparsity options:
        integer :: sparsity_mode = NLESOLVER_SPARSITY_DENSE  !! sparsity mode:
                                                             !!
                                                             !! * 1 - assume dense (use dense solver).
                                                             !! * 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 a user provided sparse solver).
        integer :: n_nonzeros = -1 !! number of nonzero Jacobian elements (used for `sparsity_mode > 1`)
        integer,dimension(:),allocatable :: irow !! sparsity pattern nonzero elements row indices.
        integer,dimension(:),allocatable :: icol !! sparsity pattern nonzero elements column indices

        ! LSQR parameters:
        real(wp) :: atol    = zero  !! relative error in definition of `A`
        real(wp) :: btol    = zero  !! relative error in definition of `b`
        real(wp) :: conlim  = zero  !! An upper limit on `cond(Abar)`, the apparent
                                    !! condition number of the matrix `Abar`.
        integer  :: itnlim  = 100   !! max iterations
        integer  :: nout    = 0     !! output unit for printing
        real(wp) :: damp    = zero  !! damp parameter for LSQR

        integer :: localSize = 0 !! 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.

        ! LUSOL parameters:
        integer :: lusol_method = 0 !! * 0 => TPP: Threshold Partial   Pivoting.
                                    !! * 1 => TRP: Threshold Rook      Pivoting.
                                    !! * 2 => TCP: Threshold Complete  Pivoting.

        ! dense version:
        procedure(grad_func),pointer :: grad => null() !! user-supplied routine to compute the gradient of the function (dense version)

        ! sparse version:
        procedure(grad_func_sparse),pointer :: grad_sparse => null() !! user-supplied routine to compute the gradient of the function (sparse version)

        ! custom sparse solver:
        procedure(sparse_solver_func),pointer :: custom_solver_sparse => null() !! user-supplied sparse linear solver routine (used for `sparsity_mode=5`)

        procedure(norm_func),pointer :: norm => null() !! function for computing the norm of the `f` vector

        contains

        private

        procedure,public :: initialize => initialize_nlesolver_variables
        procedure,public :: solve      => nlesolver_solver
        procedure,public :: destroy    => destroy_nlesolver_variables
        procedure,public :: status     => get_status

        procedure :: set_status
        procedure :: adjust_x_for_bounds
        procedure :: adjust_search_direction
        procedure :: compute_next_step

        end type nlesolver_type

nlesolver_type Derived Type

Contents

Variables

Type-Bound Procedures

Source Code

type, public :: nlesolver_type

Components

Type-Bound Procedures

procedure, public :: initialize => initialize_nlesolver_variables

Arguments

procedure, public :: solve => nlesolver_solver

private subroutine nlesolver_solver(me, x)

Arguments

procedure, public :: destroy => destroy_nlesolver_variables

private subroutine destroy_nlesolver_variables(me)

Arguments

procedure, public :: status => get_status

private subroutine get_status(me, istat, message)

Arguments

procedure, private :: set_status

private subroutine set_status(me, istat, string, i, r)

Arguments

procedure, private :: adjust_x_for_bounds

private subroutine adjust_x_for_bounds(me, x)

Arguments

procedure, private :: adjust_search_direction

private subroutine adjust_search_direction(me, x, p, pnew, modified)

Arguments

procedure, private :: compute_next_step

private subroutine compute_next_step(me, xold, search_direction, alpha, modified, xnew)

Arguments

Source Code