numdiff_type – NumDiff

type, public :: numdiff_type

base type for sparsity and Jacobian computations.

Inherits

Help

Components

Type	Visibility	Attributes		Name		Initial
logical,	private		::	exception_raised	=	.false.	if true, an exception has been raised
integer,	private		::	istat	=	0	a non-zero value will cause all routine to exit. this can be set to `-1` by calling terminate.
character(len=:),	private,	allocatable	::	error_msg			error message (if `istat/=0`)
integer,	private		::	n	=	0	number of `x` variables
integer,	private		::	m	=	0	number of `f` functions
real(kind=wp),	private,	dimension(:), allocatable	::	xlow			lower bounds on `x`
real(kind=wp),	private,	dimension(:), allocatable	::	xhigh			upper bounds on `x`
logical,	private		::	print_messages	=	.true.	if true, warning messages are printed to the `error_unit` for any errors.
integer,	private		::	chunk_size	=	100	chuck size for allocating the arrays (>0)
integer,	private		::	perturb_mode	=	1	perturbation mode: 1 - perturbation is `dx=dpert`, 2 - perturbation is `dx=dpertx`, 3 - perturbation is `dx=dpert(1+x)`
real(kind=wp),	private,	dimension(:), allocatable	::	dpert			perturbation vector for `x`
logical,	private		::	partition_sparsity_pattern	=	.false.	to partition the sparsity pattern using dsm
type(sparsity_pattern),	private		::	sparsity			the sparsity pattern
real(kind=wp),	private,	dimension(:), allocatable	::	xlow_for_sparsity			lower bounds on `x` for computing sparsity (optional)
real(kind=wp),	private,	dimension(:), allocatable	::	xhigh_for_sparsity			upper bounds on `x` for computing sparsity (optional)
real(kind=wp),	private,	dimension(:), allocatable	::	dpert_for_sparsity			perturbation vector for `x` when computing sparsity for `sparsity_mode=4`
integer,	private		::	num_sparsity_points	=	3	for `sparsity_mode=4`, the number of jacobian evaluations used to estimate the sparsity pattern. See compute_sparsity_random_2
integer,	private		::	sparsity_perturb_mode	=	1	perturbation mode (if `sparsity_mode=4`): 1 - perturbation is `dx=dpert`, 2 - perturbation is `dx=dpertx`, 3 - perturbation is `dx=dpert(1+x)`
logical,	private		::	compute_linear_sparsity_pattern	=	.false.	to also compute the linear sparsity pattern
real(kind=wp),	private		::	linear_sparsity_tol	=	epsilon(1.0_wp)	the equality tolerance for derivatives to indicate a constant jacobian element (linear sparsity)
real(kind=wp),	private		::	function_precision_tol	=	epsilon(1.0_wp)	the function precision. two functions values that are the within this tolerance are considered the same value. This is used when estimating the sparsity pattern when `sparsity_mode=2` in compute_sparsity_random
integer,	private		::	mode	=	1	1 = use `meth` (specified methods), 2 = use `class` (specified class, method is selected on-the-fly).
type(finite_diff_method),	private,	dimension(:), allocatable	::	meth			the finite difference method to use compute the `n`th column of the Jacobian `size(n)`. Either this or `class` is used
integer,	private,	dimension(:), allocatable	::	class			the class of method to use to compute the `n`th column of the Jacobian `size(n)`. Either this or `meth` is used
type(meth_array),	private,	dimension(:), allocatable	::	class_meths			array of methods for the specified classes. used with `class` when `mode=2`
real(kind=wp),	private		::	eps	=	1.0e-9_wp	tolerance parameter for diff
real(kind=wp),	private		::	acc	=	0.0_wp	tolerance parameter for diff
type(function_cache),	private		::	cache			if using the function cache
logical,	private		::	use_diff	=	.false.	if we are using the Neville's process method, rather than finite differences
procedure(func),	private,	pointer	::	problem_func	=>	null()	the user-defined function
procedure(func),	private,	pointer	::	compute_function	=>	null()	compute the user-defined function this can point to the `problem_func` or, if using the cache, it points to compute_function_with_cache.
procedure(spars_f),	private,	pointer	::	compute_sparsity	=>	null()	for computing the sparsity pattern
procedure(info_f),	private,	pointer	::	info_function	=>	null()	an optional function the user can define which is called when each column of the jacobian is computed. It can be used to perform any setup operations.
procedure(jacobian_f),	private,	pointer	::	jacobian_function	=>	null()	the low-level function called by compute_jacobian that actually computes the jacobian.

Type-Bound Procedures

procedure, public :: initialize => initialize_numdiff

initialize the class

private subroutine initialize_numdiff(me, n, m, xlow, xhigh, perturb_mode, dpert, problem_func, sparsity_mode, jacobian_method, jacobian_methods, class, classes, info, chunk_size, partition_sparsity_pattern, cache_size, xlow_for_sparsity, xhigh_for_sparsity, dpert_for_sparsity, sparsity_perturb_mode, linear_sparsity_tol, function_precision_tol, num_sparsity_points)

Initialize a numdiff_type class. This must be called first.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(in)			::	n	number of `x` variables
integer,	intent(in)			::	m	number of `f` functions
real(kind=wp),	intent(in),		dimension(n)	::	xlow	lower bounds on `x`
real(kind=wp),	intent(in),		dimension(n)	::	xhigh	upper bounds on `x`
integer,	intent(in)			::	perturb_mode	perturbation mode: 1 - perturbation is `dx=dpert`, 2 - perturbation is `dx=dpertx`, 3 - perturbation is `dx=dpert(1+x)`
real(kind=wp),	intent(in),		dimension(n)	::	dpert	perturbation vector for `x`
procedure(func)				::	problem_func	the user function that defines the problem (returns `m` functions)
integer,	intent(in)			::	sparsity_mode	the sparsity computation method: 1 - assume dense, 2 - three-point simple method, 3 - will be specified by the user in a subsequent call to set_sparsity_pattern. 4 - computes a two-point jacobian at `num_sparsity_points` points.
integer,	intent(in),	optional		::	jacobian_method	`id` code for the finite difference method to use for all `n` variables. see get_finite_difference_method
integer,	intent(in),	optional,	dimension(n)	::	jacobian_methods	`id` codes for the finite difference method to use for each variable. see get_finite_difference_method
integer,	intent(in),	optional		::	class	method class for the finite difference method to use for all `n` variables. see get_finite_difference_method
integer,	intent(in),	optional,	dimension(n)	::	classes	method class for the finite difference methods to use for each variable. see get_finite_difference_method
procedure(info_f),		optional		::	info	a function the user can define which is called when each column of the jacobian is computed. It can be used to perform any setup operations.
integer,	intent(in),	optional		::	chunk_size	chunk size for allocating the arrays (must be >0) [default is 100]
logical,	intent(in),	optional		::	partition_sparsity_pattern	if the sparisty pattern is to be partitioned using DSM [default is False]
integer,	intent(in),	optional		::	cache_size	if present, this is the cache size for the function cache (default is not to enable cache)
real(kind=wp),	intent(in),	optional,	dimension(n)	::	xlow_for_sparsity	lower bounds on `x` used for sparsity computation (when `sparsity_mode` is 2). If not present, then `xlow` is used.
real(kind=wp),	intent(in),	optional,	dimension(n)	::	xhigh_for_sparsity	upper bounds on `x` used for sparsity computation (when `sparsity_mode` is 2). If not present, then `xhigh` is used.
real(kind=wp),	intent(in),	optional,	dimension(n)	::	dpert_for_sparsity	for `sparsity_mode=4`, the perturbation
integer,	intent(in),	optional		::	sparsity_perturb_mode	perturbation mode (required if `sparsity_mode=4`): 1 - perturbation is `dx=dpert`, 2 - perturbation is `dx=dpertx`, 3 - perturbation is `dx=dpert(1+x)`
real(kind=wp),	intent(in),	optional		::	linear_sparsity_tol	the equality tolerance for derivatives to indicate a constant jacobian element (linear sparsity)
real(kind=wp),	intent(in),	optional		::	function_precision_tol	the function precision. two functions values that are the within this tolerance are considered the same value. This is used when estimating the sparsity pattern when `sparsity_mode=2` in compute_sparsity_random
integer,	intent(in),	optional		::	num_sparsity_points	for `sparsity_mode=4`, the number of jacobian evaluations used to estimate the sparsity pattern.

procedure, public :: diff_initialize => initialize_numdiff_for_diff

initialize the class

private subroutine initialize_numdiff_for_diff(me, n, m, xlow, xhigh, problem_func, sparsity_mode, info, chunk_size, eps, acc, cache_size, xlow_for_sparsity, xhigh_for_sparsity, dpert_for_sparsity, sparsity_perturb_mode, linear_sparsity_tol, function_precision_tol, num_sparsity_points, print_messages)

Alternate version of initialize_numdiff routine when using diff to compute the Jacobian.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(in)			::	n	number of `x` variables
integer,	intent(in)			::	m	number of `f` functions
real(kind=wp),	intent(in),		dimension(n)	::	xlow	lower bounds on `x`
real(kind=wp),	intent(in),		dimension(n)	::	xhigh	upper bounds on `x`
procedure(func)				::	problem_func	the user function that defines the problem (returns `m` functions)
integer,	intent(in)			::	sparsity_mode	the sparsity computation method: Read more…
procedure(info_f),		optional		::	info	a function the user can define which is called when each column of the jacobian is computed. It can be used to perform any setup operations.
integer,	intent(in),	optional		::	chunk_size	chunk size for allocating the arrays (must be >0) [default is 100]
real(kind=wp),	intent(in),	optional		::	eps	tolerance parameter for diff if not present, default is `1.0e-9_wp`
real(kind=wp),	intent(in),	optional		::	acc	tolerance parameter for diff if not present, default is `0.0_wp`
integer,	intent(in),	optional		::	cache_size	if present, this is the cache size for the function cache (default is not to enable cache)
real(kind=wp),	intent(in),	optional,	dimension(n)	::	xlow_for_sparsity	lower bounds on `x` used for sparsity computation (when `sparsity_mode` is 2). If not present, then `xlow` is used.
real(kind=wp),	intent(in),	optional,	dimension(n)	::	xhigh_for_sparsity	upper bounds on `x` used for sparsity computation (when `sparsity_mode` is 2). If not present, then `xhigh` is used.
real(kind=wp),	intent(in),	optional,	dimension(n)	::	dpert_for_sparsity	required if `sparsity_mode=4`
integer,	intent(in),	optional		::	sparsity_perturb_mode	perturbation mode (required if `sparsity_mode=4`): Read more…
real(kind=wp),	intent(in),	optional		::	linear_sparsity_tol	the equality tolerance for derivatives to indicate a constant jacobian element (linear sparsity)
real(kind=wp),	intent(in),	optional		::	function_precision_tol	the function precision. two functions values that are the within this tolerance are considered the same value. This is used when estimating the sparsity pattern when `sparsity_mode=2` in compute_sparsity_random
integer,	intent(in),	optional		::	num_sparsity_points	for `sparsity_mode=4`, the number of jacobian evaluations used to estimate the sparsity pattern.
logical,	intent(in),	optional		::	print_messages	if true, print error messages to `error_unit`. default is True.

procedure, public :: compute_jacobian

main routine to compute the Jacobian using the selected options. It returns the sparse (vector) form.

private subroutine compute_jacobian(me, x, jac)

Compute the Jacobian.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(out),	dimension(:), allocatable	::	jac	sparse jacobian vector

procedure, public :: compute_jacobian_dense

return the dense size(m,n) matrix form of the Jacobian.

private subroutine compute_jacobian_dense(me, x, jac)

just a wrapper for compute_jacobian, that returns a dense (m x n) matrix.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(out),	dimension(:,:), allocatable	::	jac	the jacobian matrix

procedure, public :: compute_jacobian_times_vector

returns the product of the Jacobian matrix and an input vector

private subroutine compute_jacobian_times_vector(me, x, v, z)

Returns the product J*v, where J is the m x n Jacobian matrix and v is an n x 1 vector.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(in),	dimension(:)	::	v	vector (size `n`)
real(kind=wp),	intent(out),	dimension(:)	::	z	The product `J*v` (size `m`)

procedure, public :: destroy => destroy_numdiff_type

destroy the class

private subroutine destroy_numdiff_type(me)

destroy the numdiff_type class.

Arguments

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

procedure, public :: print_sparsity_pattern

print the sparsity pattern in vector form to a file

private subroutine print_sparsity_pattern(me, iunit)

Print the sparsity pattern in vector form (irow, icol).

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(in)			::	iunit	file unit to write to. (assumed to be already opened)

procedure, public :: print_sparsity_matrix

print the sparsity pattern in matrix form to a file

private subroutine print_sparsity_matrix(me, iunit)

Print the sparsity pattern in matrix form.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(in)			::	iunit	file unit to write to. (assumed to be already opened)

procedure, public :: set_sparsity_pattern

manually set the sparsity pattern

private subroutine set_sparsity_pattern(me, irow, icol, linear_irow, linear_icol, linear_vals, maxgrp, ngrp)

To specify the sparsity pattern directly if it is already known.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(in),		dimension(:)	::	irow	sparsity pattern nonzero elements row indices
integer,	intent(in),		dimension(:)	::	icol	sparsity pattern nonzero elements column indices
integer,	intent(in),	optional,	dimension(:)	::	linear_irow	linear sparsity pattern nonzero elements row indices
integer,	intent(in),	optional,	dimension(:)	::	linear_icol	linear sparsity pattern nonzero elements column indices
real(kind=wp),	intent(in),	optional,	dimension(:)	::	linear_vals	linear sparsity values (constant elements of the Jacobian)
integer,	intent(in),	optional		::	maxgrp	DSM sparsity partition [only used if `me%partition_sparsity_pattern=True`]
integer,	intent(in),	optional,	dimension(:)	::	ngrp	DSM sparsity partition (size `n`) [only used if `me%partition_sparsity_pattern=True`]

procedure, public :: select_finite_diff_method

select a method in a specified class so that the variable bounds are not violated when by the perturbations.

private subroutine select_finite_diff_method(me, x, xlow, xhigh, dx, list_of_methods, fd, status_ok)

Select a finite diff method of a given class so that the perturbations of x will not violate the variable bounds.

Arguments

Type	Intent		Name
class(numdiff_type),	intent(inout)	::	me
real(kind=wp),	intent(in)	::	x	the variable value
real(kind=wp),	intent(in)	::	xlow	the variable lower bound
real(kind=wp),	intent(in)	::	xhigh	the variable upper bound
real(kind=wp),	intent(in)	::	dx	the perturbation value (>0)
type(meth_array),	intent(in)	::	list_of_methods	list of available methods to choose from
type(finite_diff_method),	intent(out)	::	fd	this method can be used
logical,	intent(out)	::	status_ok	true if it really doesn't violate the bounds (say, the bounds are very close to each other) if `status_ok=False`, then the first method in the given class is returned in `fd`.

procedure, public :: select_finite_diff_method_for_partition_group

version of select_finite_diff_method for partitioned sparsity pattern.

private subroutine select_finite_diff_method_for_partition_group(me, x, xlow, xhigh, dx, list_of_methods, fd, status_ok)

Select a finite diff method of a given class so that the perturbations of x will not violate the variable bounds for any variable in the group.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	x	the variable values
real(kind=wp),	intent(in),	dimension(:)	::	xlow	the variable lower bounds
real(kind=wp),	intent(in),	dimension(:)	::	xhigh	the variable upper bounds
real(kind=wp),	intent(in),	dimension(:)	::	dx	the perturbation values (>0)
type(meth_array),	intent(in)		::	list_of_methods	list of available methods to choose from
type(finite_diff_method),	intent(out)		::	fd	this method can be used
logical,	intent(out)		::	status_ok	true if it really doesn't violate the bounds (say, the bounds are very close to each other) if `status_ok=False`, then the first method in the given class is returned in `fd`.

procedure, public :: set_numdiff_bounds

can be called to change the variable bounds.

private subroutine set_numdiff_bounds(me, xlow, xhigh)

Change the variable bounds in a numdiff_type.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	xlow	lower bounds on `x`
real(kind=wp),	intent(in),	dimension(:)	::	xhigh	upper bounds on `x`

procedure, public :: compute_sparsity_pattern

if the user needs to compute the sparsity pattern manually. (otherwise it will be done the first time the Jacobian is computed)

private subroutine compute_sparsity_pattern(me, x, irow, icol, linear_irow, linear_icol, linear_vals)

Computes the sparsity pattern and return it. Uses the settings currently in the class.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
real(kind=wp),	intent(in),		dimension(:)	::	x	vector of variables (size `n`)
integer,	intent(out),		dimension(:), allocatable	::	irow	sparsity pattern nonzero elements row indices
integer,	intent(out),		dimension(:), allocatable	::	icol	sparsity pattern nonzero elements column indices
integer,	intent(out),	optional,	dimension(:), allocatable	::	linear_irow	linear sparsity pattern nonzero elements row indices
integer,	intent(out),	optional,	dimension(:), allocatable	::	linear_icol	linear sparsity pattern nonzero elements column indices
real(kind=wp),	intent(out),	optional,	dimension(:), allocatable	::	linear_vals	linear sparsity values (constant elements of the Jacobian)

procedure, public :: get_sparsity_pattern

returns the sparsity pattern (if it is allocated)

private subroutine get_sparsity_pattern(me, irow, icol, linear_irow, linear_icol, linear_vals, maxgrp, ngrp)

Returns the sparsity pattern from the class. If it hasn't been computed, the output arrays will not be allocated.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(out),		dimension(:), allocatable	::	irow	sparsity pattern nonzero elements row indices
integer,	intent(out),		dimension(:), allocatable	::	icol	sparsity pattern nonzero elements column indices
integer,	intent(out),	optional,	dimension(:), allocatable	::	linear_irow	linear sparsity pattern nonzero elements row indices
integer,	intent(out),	optional,	dimension(:), allocatable	::	linear_icol	linear sparsity pattern nonzero elements column indices
real(kind=wp),	intent(out),	optional,	dimension(:), allocatable	::	linear_vals	linear sparsity values (constant elements of the Jacobian)
integer,	intent(out),	optional		::	maxgrp	DSM sparsity partition
integer,	intent(out),	optional,	dimension(:), allocatable	::	ngrp	DSM sparsity partition

procedure, public :: terminate

can be called by user to stop the computation

private subroutine terminate(me)

A user-callable routine. When called, it will terminate all computations and return. The istat return code will be set to -1. This can be called in the function or the info function.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me

procedure, public :: failed

to check if an exception was raised.

private pure function failed(me)

Returns True if an exception has been raised.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(in)			::	me

Return Value logical

procedure, public :: get_error_status

the status of error condition

private subroutine get_error_status(me, istat, error_msg)

Returns the current error code and message.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(in)			::	me
integer,	intent(out),	optional		::	istat	error code (`istat=0` means no errors).
character(len=:),	intent(out),	optional,	allocatable	::	error_msg	error message string.

procedure, public :: set_dpert

change the dpert values

private subroutine set_dpert(me, dpert)

Change the dpert vector. Can be used after the class has been initialized to change the perturbation step sizes (e.g., after an iteration).

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
real(kind=wp),	intent(in),		dimension(:)	::	dpert	perturbation vector for `x`

procedure, private :: destroy_sparsity_pattern

destroy the sparsity pattern

private subroutine destroy_sparsity_pattern(me)

Destroy the sparsity pattern in the class.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me

procedure, private :: compute_perturb_vector

private subroutine compute_perturb_vector(me, x, dpert, perturb_mode, dx)

Compute dx, the perturbation vector for x

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(me%n)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(in),	dimension(me%n)	::	dpert
integer,	intent(in)		::	perturb_mode
real(kind=wp),	intent(out),	dimension(me%n)	::	dx	absolute perturbation (>0) for each variable

procedure, private :: compute_perturbation_vector

computes the variable perturbation factor

private subroutine compute_perturbation_vector(me, x, dx)

Compute dx, the perturbation vector for x used when computing the gradients.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(me%n)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(out),	dimension(me%n)	::	dx	absolute perturbation (>0) for each variable

procedure, private :: compute_sparsity_perturbation_vector

private subroutine compute_sparsity_perturbation_vector(me, x, dx)

Compute dx, the perturbation vector for x used when computing the sparsity pattern.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(me%n)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(out),	dimension(me%n)	::	dx	absolute perturbation (>0) for each variable

procedure, private :: perturb_x_and_compute_f

private subroutine perturb_x_and_compute_f(me, x, dx_factor, dx, df_factor, column, idx, df)

Perturb the specified optimization variable, and compute the function. This routine is designed so that df is accumulated as each function evaluation is done, to avoid having to allocate more temporary storage.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	x	nominal variable vector
real(kind=wp),	intent(in)		::	dx_factor	factor to multiply `dx`
real(kind=wp),	intent(in),	dimension(:)	::	dx	the perturbation value for this column
real(kind=wp),	intent(in)		::	df_factor	factor to multiply function value
integer,	intent(in)		::	column	the variable to perturb
integer,	intent(in),	dimension(:)	::	idx	the elements in this column of the Jacobian to compute (passed to function)
real(kind=wp),	intent(inout),	dimension(me%m)	::	df	the accumulated function value note: for the first call, this should be set to zero

procedure, private :: perturb_x_and_compute_f_partitioned

private subroutine perturb_x_and_compute_f_partitioned(me, x, dx_factor, dx, df_factor, columns, idx, df)

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	x	nominal variable vector
real(kind=wp),	intent(in)		::	dx_factor	factor to multiply `dx`
real(kind=wp),	intent(in),	dimension(:)	::	dx	the perturbation value for this column
real(kind=wp),	intent(in)		::	df_factor	factor to multiply function value
integer,	intent(in),	dimension(:)	::	columns	the variables to perturb
integer,	intent(in),	dimension(:)	::	idx	the elements in this column of the Jacobian to compute (passed to function)
real(kind=wp),	intent(inout),	dimension(me%m)	::	df	the accumulated function value note: for the first call, this should be set to zero

procedure, private :: set_numdiff_sparsity_bounds

private subroutine set_numdiff_sparsity_bounds(me, xlow, xhigh)

Sets the variable bounds for sparsity in a numdiff_type. These are only used for sparsity_mode=2.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
real(kind=wp),	intent(in),	optional,	dimension(:)	::	xlow	lower bounds on `x` to be used for sparsity computation. If not present, then then `xlow` values in the class are used.
real(kind=wp),	intent(in),	optional,	dimension(:)	::	xhigh	upper bounds on `x` to be used for sparsity computation. If not present, then then `xhigh` values in the class are used.

procedure, private :: set_sparsity_mode

private subroutine set_sparsity_mode(me, sparsity_mode, xlow_for_sparsity, xhigh_for_sparsity)

Set sparsity mode.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me
integer,	intent(in)			::	sparsity_mode	the sparsity computation method: 1 - assume dense, 2 - three-point simple method, 3 - will be specified by the user in a subsequent call to set_sparsity_pattern. 4 - computes a two-point jacobian at `num_sparsity_points` points.
real(kind=wp),	intent(in),	optional,	dimension(:)	::	xlow_for_sparsity	lower bounds on `x` used for sparsity computation (when `sparsity_mode` is 2). If not present, then `xlow` is used.
real(kind=wp),	intent(in),	optional,	dimension(:)	::	xhigh_for_sparsity	upper bounds on `x` used for sparsity computation (when `sparsity_mode` is 2). If not present, then `xhigh` is used.

procedure, private :: generate_dense_sparsity_partition

private subroutine generate_dense_sparsity_partition(me)

Generate a "dense" sparsity partition.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me

procedure, private :: compute_jacobian_for_sparsity

private subroutine compute_jacobian_for_sparsity(me, i, class_meths, x, jac)

A separate version of compute_jacobian to be used only when computing the sparsity pattern in compute_sparsity_random_2. It uses class_meths and the sparsity dperts and bounds.

Arguments

Type	Intent	Attributes		Name
class(numdiff_type),	intent(inout)		::	me
integer,	intent(in)		::	i	the column being computed
type(meth_array),	intent(in)		::	class_meths	set of finite diff methods to use
real(kind=wp),	intent(in),	dimension(:)	::	x	vector of variables (size `n`)
real(kind=wp),	intent(out),	dimension(:), allocatable	::	jac	sparse jacobian vector

procedure, private :: resize_sparsity_vectors

private subroutine resize_sparsity_vectors(me, n_icol, n_irow, n_linear_icol, n_linear_irow, n_linear_vals)

Resize the sparsity arrays after accumulating them.

Arguments

Type	Intent		Name
class(numdiff_type),	intent(inout)	::	me
integer,	intent(inout)	::	n_icol
integer,	intent(inout)	::	n_irow
integer,	intent(inout)	::	n_linear_icol
integer,	intent(inout)	::	n_linear_irow
integer,	intent(inout)	::	n_linear_vals

procedure, private :: raise_exception

private subroutine raise_exception(me, istat, routine, error_msg)

Raise an exception.

Arguments

Type	Intent		Name
class(numdiff_type),	intent(inout)	::	me
integer,	intent(in)	::	istat	error code.
character(len=*),	intent(in)	::	routine	the routine where the error was raised.
character(len=*),	intent(in)	::	error_msg	error message string.

procedure, private :: clear_exceptions

private subroutine clear_exceptions(me)

Clear all exceptions.

Arguments

Type	Intent	Optional	Attributes		Name
class(numdiff_type),	intent(inout)			::	me

Source Code

    type,public :: numdiff_type

        !! base type for sparsity and Jacobian computations.

        private

        logical :: exception_raised = .false. !! if true, an exception has been raised
        integer :: istat = 0   !! a non-zero value will cause all routine to exit.
                               !! this can be set to `-1` by calling [[terminate]].
        character(len=:),allocatable :: error_msg !! error message (if `istat/=0`)

        integer :: n = 0 !! number of `x` variables
        integer :: m = 0 !! number of `f` functions

        real(wp),dimension(:),allocatable :: xlow  !! lower bounds on `x`
        real(wp),dimension(:),allocatable :: xhigh !! upper bounds on `x`

        logical :: print_messages = .true. !! if true, warning messages are printed
                                           !! to the `error_unit` for any errors.

        integer :: chunk_size = 100  !! chuck size for allocating the arrays (>0)

        integer :: perturb_mode = 1  !! perturbation mode:
                                     !!
                                     !! **1** - perturbation is `dx=dpert`,
                                     !! **2** - perturbation is `dx=dpert*x`,
                                     !! **3** - perturbation is `dx=dpert*(1+x)`
        real(wp),dimension(:),allocatable :: dpert !! perturbation vector for `x`

        logical :: partition_sparsity_pattern = .false.  !! to partition the sparsity pattern using [[dsm]]
        type(sparsity_pattern) :: sparsity  !! the sparsity pattern
        real(wp),dimension(:),allocatable :: xlow_for_sparsity  !! lower bounds on `x` for computing sparsity (optional)
        real(wp),dimension(:),allocatable :: xhigh_for_sparsity !! upper bounds on `x` for computing sparsity (optional)
        real(wp),dimension(:),allocatable :: dpert_for_sparsity !! perturbation vector for `x` when computing
                                                                !! sparsity for `sparsity_mode=4`
        integer :: num_sparsity_points = 3    !! for `sparsity_mode=4`, the number of jacobian
                                              !! evaluations used to estimate the sparsity pattern.
                                              !! See [[compute_sparsity_random_2]]
        integer :: sparsity_perturb_mode = 1  !! perturbation mode (if `sparsity_mode=4`):
                                              !!
                                              !! **1** - perturbation is `dx=dpert`,
                                              !! **2** - perturbation is `dx=dpert*x`,
                                              !! **3** - perturbation is `dx=dpert*(1+x)`

        ! if computing the sparsity pattern, we also have an option to
        ! compute the linear pattern, which indicates constant elements
        ! of the jacobian. these elements don't need to be computed again.
        logical :: compute_linear_sparsity_pattern = .false.  !! to also compute the linear sparsity pattern
        real(wp) :: linear_sparsity_tol = epsilon(1.0_wp)    !! the equality tolerance for derivatives to
                                                             !! indicate a constant jacobian element (linear sparsity)
        real(wp) :: function_precision_tol = epsilon(1.0_wp) !! the function precision. two functions values
                                                             !! that are the within this tolerance are
                                                             !! considered the same value. This is used
                                                             !! when estimating the sparsity pattern when
                                                             !! `sparsity_mode=2` in [[compute_sparsity_random]]

        integer :: mode = 1 !! **1** = use `meth` (specified methods),
                            !! **2** = use `class` (specified class, method is selected on-the-fly).
        type(finite_diff_method),dimension(:),allocatable :: meth   !! the finite difference method to use
                                                                    !! compute the `n`th column of the Jacobian
                                                                    !! `size(n)`.  Either this or `class` is used
        integer,dimension(:),allocatable :: class  !! the class of method to use to
                                                   !! compute the `n`th column of the Jacobian
                                                   !! `size(n)`. Either this or `meth` is used
        type(meth_array),dimension(:),allocatable :: class_meths !! array of methods for the specified classes.
                                                                 !! used with `class` when `mode=2`

        ! parameters when using diff:
        real(wp) :: eps = 1.0e-9_wp !! tolerance parameter for [[diff]]
        real(wp) :: acc = 0.0_wp    !! tolerance parameter for [[diff]]

        type(function_cache) :: cache    !! if using the function cache

        logical :: use_diff = .false. !! if we are using the Neville's process method,
                                      !! rather than finite differences

        procedure(func),pointer :: problem_func => null()
            !! the user-defined function

        ! these are required to be defined by the user:
        procedure(func),pointer :: compute_function => null()
            !! compute the user-defined function
            !! this can point to the `problem_func` or, if using
            !! the cache, it points to [[compute_function_with_cache]].

        procedure(spars_f),pointer :: compute_sparsity => null()
            !! for computing the sparsity pattern

        procedure(info_f),pointer :: info_function => null()
            !! an optional function the user can define
            !! which is called when each column of the jacobian is computed.
            !! It can be used to perform any setup operations.

        procedure(jacobian_f),pointer :: jacobian_function => null()
            !! the low-level function called by [[compute_jacobian]]
            !! that actually computes the jacobian.

    contains

        private

        procedure,public :: initialize => initialize_numdiff !! initialize the class
        procedure,public :: diff_initialize => initialize_numdiff_for_diff !! initialize the class

        procedure,public :: compute_jacobian              !! main routine to compute the Jacobian
                                                          !! using the selected options. It
                                                          !! returns the sparse (vector) form.
        procedure,public :: compute_jacobian_dense        !! return the dense `size(m,n)`
                                                          !! matrix form of the Jacobian.
        procedure,public :: compute_jacobian_times_vector !! returns the product of the Jacobian
                                                          !! matrix and an input vector
        procedure,public :: destroy => destroy_numdiff_type  !! destroy the class
        procedure,public :: print_sparsity_pattern  !! print the sparsity pattern in vector form to a file
        procedure,public :: print_sparsity_matrix   !! print the sparsity pattern in matrix form to a file
        procedure,public :: set_sparsity_pattern    !! manually set the sparsity pattern
        procedure,public :: select_finite_diff_method  !! select a method in a specified class so
                                                       !! that the variable bounds are not violated
                                                       !! when by the perturbations.
        procedure,public :: select_finite_diff_method_for_partition_group  !! version of [[select_finite_diff_method]]
                                                                           !! for partitioned sparsity pattern.
        procedure,public :: set_numdiff_bounds  !! can be called to change the variable bounds.

        procedure,public :: compute_sparsity_pattern !! if the user needs to compute the sparsity pattern manually.
                                                     !! (otherwise it will be done the first time the Jacobian is computed)
        procedure,public :: get_sparsity_pattern     !! returns the sparsity pattern (if it is allocated)

        procedure,public :: terminate !! can be called by user to stop the computation

        procedure,public :: failed !! to check if an exception was raised.
        procedure,public :: get_error_status !! the status of error condition
        procedure,public :: set_dpert !! change the `dpert` values

        ! internal routines:
        procedure :: destroy_sparsity_pattern      !! destroy the sparsity pattern
        procedure :: compute_perturb_vector
        procedure :: compute_perturbation_vector   !! computes the variable perturbation factor
        procedure :: compute_sparsity_perturbation_vector
        procedure :: perturb_x_and_compute_f
        procedure :: perturb_x_and_compute_f_partitioned
        procedure :: set_numdiff_sparsity_bounds
        procedure :: set_sparsity_mode
        procedure :: generate_dense_sparsity_partition
        procedure :: compute_jacobian_for_sparsity
        procedure :: resize_sparsity_vectors
        procedure :: raise_exception
        procedure :: clear_exceptions

    end type numdiff_type

numdiff_type Derived Type

Contents

Variables

Type-Bound Procedures

Source Code

type, public :: numdiff_type

Inherits

Components

Type-Bound Procedures

procedure, public :: initialize => initialize_numdiff

Arguments

procedure, public :: diff_initialize => initialize_numdiff_for_diff

Arguments

procedure, public :: compute_jacobian

private subroutine compute_jacobian(me, x, jac)

Arguments

procedure, public :: compute_jacobian_dense

private subroutine compute_jacobian_dense(me, x, jac)

Arguments

procedure, public :: compute_jacobian_times_vector

private subroutine compute_jacobian_times_vector(me, x, v, z)

Arguments

procedure, public :: destroy => destroy_numdiff_type

private subroutine destroy_numdiff_type(me)

Arguments

procedure, public :: print_sparsity_pattern

private subroutine print_sparsity_pattern(me, iunit)

Arguments

procedure, public :: print_sparsity_matrix

private subroutine print_sparsity_matrix(me, iunit)

Arguments

procedure, public :: set_sparsity_pattern

private subroutine set_sparsity_pattern(me, irow, icol, linear_irow, linear_icol, linear_vals, maxgrp, ngrp)

Arguments

procedure, public :: select_finite_diff_method

private subroutine select_finite_diff_method(me, x, xlow, xhigh, dx, list_of_methods, fd, status_ok)

Arguments

procedure, public :: select_finite_diff_method_for_partition_group

private subroutine select_finite_diff_method_for_partition_group(me, x, xlow, xhigh, dx, list_of_methods, fd, status_ok)

Arguments

procedure, public :: set_numdiff_bounds

private subroutine set_numdiff_bounds(me, xlow, xhigh)

Arguments

procedure, public :: compute_sparsity_pattern

private subroutine compute_sparsity_pattern(me, x, irow, icol, linear_irow, linear_icol, linear_vals)

Arguments

procedure, public :: get_sparsity_pattern

private subroutine get_sparsity_pattern(me, irow, icol, linear_irow, linear_icol, linear_vals, maxgrp, ngrp)

Arguments

procedure, public :: terminate

private subroutine terminate(me)

Arguments

procedure, public :: failed

private pure function failed(me)

Arguments

Return Value logical

procedure, public :: get_error_status

private subroutine get_error_status(me, istat, error_msg)

Arguments

procedure, public :: set_dpert

private subroutine set_dpert(me, dpert)

Arguments

procedure, private :: destroy_sparsity_pattern

private subroutine destroy_sparsity_pattern(me)

Arguments

procedure, private :: compute_perturb_vector

private subroutine compute_perturb_vector(me, x, dpert, perturb_mode, dx)

Arguments

procedure, private :: compute_perturbation_vector

private subroutine compute_perturbation_vector(me, x, dx)

Arguments

procedure, private :: compute_sparsity_perturbation_vector

private subroutine compute_sparsity_perturbation_vector(me, x, dx)

Arguments

procedure, private :: perturb_x_and_compute_f

private subroutine perturb_x_and_compute_f(me, x, dx_factor, dx, df_factor, column, idx, df)

Arguments

procedure, private :: perturb_x_and_compute_f_partitioned

private subroutine perturb_x_and_compute_f_partitioned(me, x, dx_factor, dx, df_factor, columns, idx, df)

Arguments