nlesolver_module Module

A basic multidimensional nonlinear equation solver, using a Newton-Raphson type direct method.

Features

  • Works with square, under-determined, or over-determined systems.
  • Uses LAPACK routines (dgesv or dgels) to solve the linear system.
  • Has a Broyden update option.
  • Has various line search options.

References

  • https://en.wikipedia.org/wiki/Backtracking_line_search
  • http://projecteuclid.org/download/pdf_1/euclid.pjm/1102995080
  • http://help.agi.com/stk/index.htm#gator/eq-diffcorr.htm
  • http://gmat.sourceforge.net/doc/R2015a/html/DifferentialCorrector.html
  • https://openmdao.org/newdocs/versions/latest/features/building_blocks/solvers/bounds_enforce.html

Author

  • Jacob Williams

License

  • BSD-3

Todo

add an istat output to func and grad, for user stopping or to take a smaller stop (if istat>0 take a smaller step, if istat<0 abort)

Note

The default real kind (wp) can be changed using optional preprocessor flags. This library was built with real kind: real(kind=real64) [8 bytes]


Uses

  • module~~nlesolver_module~~UsesGraph module~nlesolver_module nlesolver_module fmin_module fmin_module module~nlesolver_module->fmin_module iso_fortran_env iso_fortran_env module~nlesolver_module->iso_fortran_env lsmrModule lsmrModule module~nlesolver_module->lsmrModule lsqr_module lsqr_module module~nlesolver_module->lsqr_module lusol_ez_module lusol_ez_module module~nlesolver_module->lusol_ez_module

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: nlesolver_rk = real64

real kind used by this module [8 bytes]
integer, private, parameter :: wp = nlesolver_rk

local copy of nlesolver_rk with a shorter name
real(kind=wp), private, parameter :: zero = 0.0_wp
real(kind=wp), private, parameter :: one = 1.0_wp
real(kind=wp), private, parameter :: two = 2.0_wp
real(kind=wp), private, parameter :: eps = epsilon(one)

machine

real(kind=wp), private, parameter :: big = huge(one)
integer, public, parameter :: NLESOLVER_SPARSITY_DENSE = 1

sparsity_mode : assume dense (use dense solver).
integer, public, parameter :: NLESOLVER_SPARSITY_LSQR = 2

sparsity_mode : assume sparse (use LSQR sparse solver).
integer, public, parameter :: NLESOLVER_SPARSITY_LUSOL = 3

sparsity_mode : assume sparse (use LUSOL sparse solver).
integer, public, parameter :: NLESOLVER_SPARSITY_LSMR = 4

sparsity_mode : assume sparse (use LSMR sparse solver).
integer, public, parameter :: NLESOLVER_SPARSITY_CUSTOM_SPARSE = 5

sparsity_mode : assume sparse (use a user provided sparse solver).
integer, public, parameter :: NLESOLVER_IGNORE_BOUNDS = 0

bounds_mode : ignore bounds
integer, public, parameter :: NLESOLVER_SCALAR_BOUNDS = 1

bounds_mode : scalar mode
integer, public, parameter :: NLESOLVER_VECTOR_BOUNDS = 2

bounds_mode : vector mode
integer, public, parameter :: NLESOLVER_WALL_BOUNDS = 3

bounds_mode : wall mode
integer, public, parameter :: NLESOLVER_2_NORM = 1

norm_mode : 2-norm
integer, public, parameter :: NLESOLVER_INF_NORM = 2

norm_mode : infinity-norm
integer, public, parameter :: NLESOLVER_1_NORM = 3

norm_mode : 1-norm
integer, public, parameter :: NLESOLVER_LINESEARCH_SIMPLE = 1

linesearch : use the specified alpha (0,1]
integer, public, parameter :: NLESOLVER_LINESEARCH_BACKTRACKING = 2

linesearch : backtracking linesearch (between alpha_min and alpha_max)
integer, public, parameter :: NLESOLVER_LINESEARCH_EXACT = 3

linesearch : exact linesearch (between alpha_min and alpha_max)
integer, public, parameter :: NLESOLVER_LINESEARCH_FIXEDPOINT = 4

linesearch : evaluate function at specified fixed points (between alpha_min and alpha_max)

Abstract Interfaces

abstract interface

  • private pure function norm_func(me, fvec) result(f)

    function vector norm.

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(in) :: me
    real(kind=wp), intent(in), dimension(me%m) :: fvec

    the function vector

    Return Value real(kind=wp)

    norm of the vector

abstract interface

  • private subroutine sparse_solver_func(me, n_cols, n_rows, n_nonzero, irow, icol, a, b, x, istat)

    for a custom user-provided linear solver (sparse version) solve Ax = b for x, given A and b.

    Arguments

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

    n: number of columns in A.
    integer, intent(in) :: n_rows

    m: number of rows in A.
    integer, intent(in) :: n_nonzero

    number of nonzero elements of A.
    integer, intent(in), dimension(n_nonzero) :: irow

    sparsity pattern (size is n_nonzero)
    integer, intent(in), dimension(n_nonzero) :: icol

    sparsity pattern (size is n_nonzero)
    real(kind=wp), intent(in), dimension(n_nonzero) :: a

    matrix elements (size is n_nonzero)
    real(kind=wp), intent(in), dimension(n_rows) :: b

    right hand side (size is m)
    real(kind=wp), intent(out), dimension(n_cols) :: x

    solution (size is n)
    integer, intent(out) :: istat

    status code (=0 for success)

abstract interface

  • private subroutine func_func(me, x, f)

    compute the function

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(inout) :: me
    real(kind=wp), intent(in), dimension(:) :: x
    real(kind=wp), intent(out), dimension(:) :: f

abstract interface

  • private subroutine grad_func(me, x, g)

    compute the gradient of the function (Jacobian). Dense version.

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(inout) :: me
    real(kind=wp), intent(in), dimension(:) :: x
    real(kind=wp), intent(out), dimension(:,:) :: g

abstract interface

  • private subroutine grad_func_sparse(me, x, g)

    compute the gradient of the function (Jacobian). Sparse version.

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(inout) :: me
    real(kind=wp), intent(in), dimension(:) :: x
    real(kind=wp), intent(out), dimension(:) :: g

    sparse jacobian. length is n_nonzeros

abstract interface

  • private subroutine export_func(me, x, f, iter)

    export an iteration:

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(inout) :: me
    real(kind=wp), intent(in), dimension(:) :: x
    real(kind=wp), intent(in), dimension(:) :: f
    integer, intent(in) :: iter

    iteration number

abstract interface

  • private subroutine wait_func(me, user_stop)

    enable a user-triggered stop of the iterations:

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(inout) :: me
    logical, intent(out) :: user_stop

abstract interface

  • private subroutine linesearch_func(me, xold, p, x, f, fvec, fjac, fjac_sparse)

    line search method. Note that not all inputs/outputs are used by all methods.

    Arguments

    Type IntentOptional Attributes Name
    class(nlesolver_type), intent(inout) :: me
    real(kind=wp), intent(in), dimension(me%n) :: xold

    previous value of x
    real(kind=wp), intent(in), dimension(me%n) :: p

    search direction
    real(kind=wp), intent(out), dimension(me%n) :: x

    new x
    real(kind=wp), intent(inout) :: f
    • input: current magnitude of fvec,
    • output: new value of f
    real(kind=wp), intent(inout), dimension(me%m) :: fvec
    • input: current function vector,
    • output: new function vector
    real(kind=wp), intent(in), optional, dimension(:,:) :: fjac

    jacobian matrix [dense]
    real(kind=wp), intent(in), optional, dimension(:) :: fjac_sparse

    jacobian matrix [sparse]

Derived Types

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:

Read more…
integer, private :: bounds_mode = NLESOLVER_IGNORE_BOUNDS

how to handle the x variable bounds:

Read more…
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:

Read more…
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 This many n-vectors "v" (the most recent ones) are saved for reorthogonalizing the next v.

Read more…
integer, private :: lusol_method = 0 Read more…
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
procedure, public :: solve => nlesolver_solver
procedure, public :: destroy => destroy_nlesolver_variables
procedure, public :: status => get_status
procedure, private :: set_status
procedure, private :: adjust_x_for_bounds
procedure, private :: adjust_search_direction
procedure, private :: compute_next_step

Functions

private function int2str(i) result(s)

Convert an integer to a string. Works for up to 256 digits.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: i

integer to convert

Return Value character(len=:), allocatable

string result

private function real2str(r) result(s)

Convert a real to a string. Works for up to 256 digits.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: r

real to convert

Return Value character(len=:), allocatable

string result

private pure function norm_2(me, fvec) result(f)

2-norm function

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(in) :: me
real(kind=wp), intent(in), dimension(me%m) :: fvec

the function vector

Return Value real(kind=wp)

norm of the vector

private pure function norm_1(me, fvec) result(f)

1-norm function

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(in) :: me
real(kind=wp), intent(in), dimension(me%m) :: fvec

the function vector

Return Value real(kind=wp)

norm of the vector

private pure function norm_inf(me, fvec) result(f)

Infinity-norm function

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(in) :: me
real(kind=wp), intent(in), dimension(me%m) :: fvec

the function vector

Return Value real(kind=wp)

norm of the vector


Subroutines

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

Set status flag and message.

Arguments

Type IntentOptional Attributes 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

private subroutine get_status(me, istat, message)

Return the status code and message from the nlesolver_type class.

Read more…

Arguments

Type IntentOptional 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

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 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:

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…

private subroutine nlesolver_solver(me, x)

Main solver.

Arguments

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

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 IntentOptional Attributes Name
class(nlesolver_type), intent(inout) :: me
real(kind=wp), intent(inout), dimension(me%n) :: x

the x vector to adjust

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

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

Read more…

Arguments

Type IntentOptional 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

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

Compute the next step.

Arguments

Type IntentOptional 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

private subroutine destroy_nlesolver_variables(me)

Destructor

Arguments

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

private subroutine linear_solver(m, n, a, b, x, info)

Solve the linear system: , using a dense, direct method.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: m

number of rows in a
integer, intent(in) :: n

number of columns in a
real(kind=wp), intent(in), dimension(m,n) :: a

A matrix of the linear system
real(kind=wp), intent(in), dimension(m) :: b

RHS of the linear system
real(kind=wp), intent(out), dimension(n) :: x

the solution of the linear system.
integer, intent(out) :: info

output status flag (=0 if success)

private subroutine simple_step(me, xold, p, x, f, fvec, fjac, fjac_sparse)

Take a simple step in the search direction of p * alpha.

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(inout) :: me
real(kind=wp), intent(in), dimension(me%n) :: xold

previous value of x
real(kind=wp), intent(in), dimension(me%n) :: p

search direction
real(kind=wp), intent(out), dimension(me%n) :: x

new x
real(kind=wp), intent(inout) :: f

magnitude of fvec
real(kind=wp), intent(inout), dimension(me%m) :: fvec

function vector
real(kind=wp), intent(in), optional, dimension(:,:) :: fjac

jacobian matrix [dense]
real(kind=wp), intent(in), optional, dimension(:) :: fjac_sparse

jacobian matrix [sparse]

private subroutine backtracking_linesearch(me, xold, p, x, f, fvec, fjac, fjac_sparse)

Backtracking line search.

Read more…

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(inout) :: me
real(kind=wp), intent(in), dimension(me%n) :: xold

previous value of x
real(kind=wp), intent(in), dimension(me%n) :: p

search direction
real(kind=wp), intent(out), dimension(me%n) :: x

new x
real(kind=wp), intent(inout) :: f

magnitude of fvec
real(kind=wp), intent(inout), dimension(me%m) :: fvec

function vector
real(kind=wp), intent(in), optional, dimension(:,:) :: fjac

jacobian matrix [dense]
real(kind=wp), intent(in), optional, dimension(:) :: fjac_sparse

jacobian matrix [sparse]

private subroutine exact_linesearch(me, xold, p, x, f, fvec, fjac, fjac_sparse)

An exact linesearch that uses a derivative-free minimizer to find the minimum value of f(x) between x = xold + p * alpha_min and x = xold + p * alpha_max.

Read more…

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(inout) :: me
real(kind=wp), intent(in), dimension(me%n) :: xold

previous value of x
real(kind=wp), intent(in), dimension(me%n) :: p

search direction
real(kind=wp), intent(out), dimension(me%n) :: x

new x
real(kind=wp), intent(inout) :: f

magnitude of fvec
real(kind=wp), intent(inout), dimension(me%m) :: fvec

function vector
real(kind=wp), intent(in), optional, dimension(:,:) :: fjac

jacobian matrix [dense]
real(kind=wp), intent(in), optional, dimension(:) :: fjac_sparse

jacobian matrix [sparse]

private subroutine fixed_point_linesearch(me, xold, p, x, f, fvec, fjac, fjac_sparse)

A simple search that just evaluates the function at a specified number of points and picks the one with the minimum function value.

Arguments

Type IntentOptional Attributes Name
class(nlesolver_type), intent(inout) :: me
real(kind=wp), intent(in), dimension(me%n) :: xold

previous value of x
real(kind=wp), intent(in), dimension(me%n) :: p

search direction
real(kind=wp), intent(out), dimension(me%n) :: x

new x
real(kind=wp), intent(inout) :: f

magnitude of fvec
real(kind=wp), intent(inout), dimension(me%m) :: fvec

function vector
real(kind=wp), intent(in), optional, dimension(:,:) :: fjac

jacobian matrix [dense]
real(kind=wp), intent(in), optional, dimension(:) :: fjac_sparse

jacobian matrix [sparse]