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

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

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

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

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

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

assume sparse (use a user provided sparse solver).

Abstract Interfaces

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

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

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)

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.

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