num_utilities Module Reference

Numerical utilities. More...

Interfaces and Types
interface	add_arr_mult
	Add to an array (3) the product of arrays (1) and (2). More...
interface	bubble_sort
	Sorting with the bubble sort routine. More...
interface	calc_det
	Calculate determinant of a matrix. More...
interface	calc_int
	Integrates a function using the trapezoidal rule. More...
interface	calc_inv
	Calculate inverse of square matrix A. More...
interface	calc_mult
	Calculate matrix multiplication of two square matrices \(\overline{\text{AB}} = \overline{\text{A}} \ \overline{\text{B}}\). More...
interface	con
	Either takes the complex conjugate of a square matrix element A, defined on a 3-D grid, or not. More...
interface	con2dis
	Convert between points from a continuous grid to a discrete grid. More...
interface	conv_mat
	Converts a (symmetric) matrix A with the storage convention described in eq_vars.eq_2_type. More...
interface	dis2con
	Convert between points from a discrete grid to a continuous grid. More...
interface	order_per_fun
	Order a periodic function to include \(0\ldots 2\pi\) and an overlap. More...
interface	round_with_tol
	Rounds an arry of values to limits, with a tolerance \(10^{-5}\) that can optionally be modified. More...
interface	spline
	Wrapper to the pspline library, making it easier to use for 1-D applications where speed is not the main priority. If spline representations are to be reused, manually use the library. More...

Functions/Subroutines
integer function	calc_inv_0d (inv_0D, A)
	private constant version More...
subroutine, public	calc_aux_utilities (n_mod)
	Initialize utilities for fast future reference, depending on program style. More...
integer function	calc_derivs_1d_id (deriv, dims)
	Returns the 1D indices for derivatives of a certain order in certain dimensions. More...
integer function, dimension(:,:), allocatable, public	derivs (order, dims)
	Returns derivatives of certain order. More...
integer function, public	check_deriv (deriv, max_deriv, sr_name)
	checks whether the derivatives requested for a certain subroutine are valid More...
integer function, public	calc_ext_var (ext_var, var, var_points, ext_point, deriv_in)
	Extrapolates a function. More...
integer function, public	calc_coeff_fin_diff (deriv, nr, ind, coeff)
	Calculate the coefficients for finite differences. More...
integer function, public	c (ij, sym, n, lim_n)
	Convert 2-D coordinates (i,j) to the storage convention used in matrices. More...
recursive integer function, public	fac (n)
	Calculate factorial. More...
integer function, public	is_sym (n, nn, sym)
	Determines whether a matrix making use of the storage convention in eq_vars.eq_2_type is symmetric or not. More...
recursive integer function, public	lcm (u, v)
	Returns common multiple using the Euclid's algorithm. More...
recursive integer function, public	gcd (u, v)
	Returns least denominator using the GCD. More...
subroutine, public	shift_f (Al, Bl, Cl, A, B, C)
	Calculate multiplication through shifting of fourier modes A and B into C. More...
subroutine, public	solve_vand (n, a, b, x, transp)
	Solve a Vandermonde system \(\overline{\text{A}} \ \vec{X} = \vec{B}\). More...
integer function, public	calc_d2_smooth (fil_N, x, y, D2y, style)
	Calculate second derivative with smoothing formula by Holoborodko, [9]. More...

Variables
integer, dimension(:,:,:), allocatable, public	d
	1-D array indices of derivatives More...
integer, dimension(:,:), allocatable, public	m
	1-D array indices of metric indices More...
integer, dimension(:,:), allocatable, public	f
	1-D array indices of Fourier mode combination indices More...
logical, public	debug_con2dis_reg = .false.
	plot debug information for con2dis_reg() More...
logical, public	debug_calc_coeff_fin_diff = .false.
	plot debug information for calc_coeff_fin_diff() More...

Detailed Description

Numerical utilities.

Function/Subroutine Documentation

c()

integer function, public num_utilities::c	(	integer, dimension(2), intent(in)	ij,
		logical, intent(in)	sym,
		integer, intent(in), optional	n,
		integer, dimension(2,2), intent(in), optional	lim_n
	)

Convert 2-D coordinates (i,j) to the storage convention used in matrices.

Their size is by default taken to be 3:

\[ \left(\begin{array}{ccc} 1 & 4 & 7 \\ 2 & 5 & 8 \\ 3 & 6 & 9 \end{array}\right) \ \text{or} \ \left(\begin{array}{ccc} 1 & & \\ 2 & 4 & \\ 3 & 5 & 6 \end{array}\right) \ \text{for symmetic matrices} \ . \]

Optionally, this can be changed using n.

The value of \(c\) is then given by

\[c = (j-1) n + i\]

for non-symmetric matrices, and by

\[\begin{aligned} c &= (j-1)n + i - (j-1)\frac{j}{2} & \text{if} \ i > j \\ c &= (i-1)n + j - (i-1)\frac{i}{2} & \text{if} \ j > i \end{aligned}\]

since \(\sum_{k=1}^{j-1} k = \frac{(j-1)j}{2}\).

For local indices in submatrices, the limits in both dimensions have to be passed. The results for the full matrix are then subtracted by amount corresponding to the left, above and below parts with respect to the submatrix:

• left:

  \[\sum_{i=1}^{m(2)-1} n-i+1 = (m(2)-1) \left(n+1-\frac{m(2)}{2}\right) \]

• above (if positive):

  \[\sum_{i=1}^j \left(m(1)-m(2)+1-i\right) = j \left(m(1)-m(2)+\frac{1}{2} - \frac{j^*}{2}\right) \ , \ \text{with} \ j^* = \text{min}\left(0,j,m(1)-m(2)+1\right) \]

• below:

  \[\sum_{i=1}^{j-1} n-M(1) = \left(n-M(1)\right) \left(j-1\right) \]

where \(m\) = min and \(M\) = max.

Note

1. The submatrix version is not fast, so results should be saved and reused.
2. No checks are done whether the indices make sense.

Parameters

[in]	ij	2D coords. (i,j)
[in]	sym	whether symmetric
[in]	n	max of 2-D coords.
[in]	lim_n	min. and max. of 2-D coords.

Definition at line 2556 of file num_utilities.f90.

Here is the caller graph for this function:

calc_aux_utilities()

subroutine, public num_utilities::calc_aux_utilities

(

integer, intent(in), optional

n_mod

)

Initialize utilities for fast future reference, depending on program style.

Utilities initialized:

• derivatives
• metrics
• Fourier modes (optionally)

If Fourier modes are also initialized, the quantity n_mod has to be provided as well.

Parameters

[in]

n_mod

n_mod for Fourier modes

Definition at line 2063 of file num_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

calc_coeff_fin_diff()

integer function, public num_utilities::calc_coeff_fin_diff	(	integer, intent(in)	deriv,
		integer, intent(in)	nr,
		integer, intent(in)	ind,
		real(dp), dimension(:), intent(inout), allocatable	coeff
	)

Calculate the coefficients for finite differences.

These represent a derivative of degree deriv at an index ind, by the weighted sum of a number nr of points with 1<=ind<=nr.

They are found by solving a Vandermonde system, multiplied by the faculty of the derivative.

The result returned in coeff(1:nr) are used in

\[\frac{d f}{d x} = \sum_{j=1}^n c(j) f(i+j) . \]

where \(i\) = ind, \(n\) = nr and \(c\) = coeff

Note

They need to be divided by the step size before usage.

Examples:

• symmetric finite differences for derivative deriv of order ord:
  • \(m = \text{ceiling}(\frac{p+d}{2})\) to guarantee the order
  • \(n = 1+2m\)
  • \(i = 1+m\)
• left finite differences for derivative deriv of order ord:
  • \(n = 1+d+p \)
  • \(i = m\) where \(p\) = ord, \(d\) = deriv and \(m\) = nr_2.

Returns

ierr

Parameters

[in]	deriv	degree of derivative
[in]	nr	number of points
[in]	ind	position of derivative
[in,out]	coeff	output coefficients

Definition at line 2449 of file num_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

calc_d2_smooth()

integer function, public num_utilities::calc_d2_smooth	(	integer, intent(in)	fil_N,
		real(dp), dimension(:), intent(in)	x,
		real(dp), dimension(:), intent(in)	y,
		real(dp), dimension(:), intent(inout)	D2y,
		integer, intent(in)	style
	)

Calculate second derivative with smoothing formula by Holoborodko, [9].

The formula used is

\(s\left(k\right) = \left[\left(2N-10\right)s\left(k+1\right) - \left(N+2k+3\right)s\left(k+2\right)\right]/\left(N-2k-1\right)\)

with \(s\left(k>M\right) = 0\) and \(\left(k=N\right) = 1\).

for a given \(N\).

For style 1, central differences are used:

\(\frac{\partial^2 y}{\partial x^2} \approx \frac{1}{2^{N-3}} \left(\sum_{k=1}^M \alpha_k \left(y_k + y_{-k} \right) - 2 y_0 \sum_{k=1}^M \alpha_k \right)\)

and for style 2, backward differences:

\(\frac{\partial^2 y}{\partial x^2} \approx \frac{1}{2^{N-3}} \left(\sum_{k=1}^M \beta_k \left(y_{-M+k} + y_{-M-k} \right) - 2 y_{-M} \sum_{k=1}^M \beta_k \right)\)

with \(\alpha_k\) and \(\beta_k\) given by

\(\begin{aligned} \alpha_k &= \frac{4k^2 s_k}{\left(x_k-x_{-k}\right)^2} \\ \beta_k &= \frac{4k^2 s_k}{\left(x_{-M+k}-x_{-M-k}\right)^2} \end{aligned}\).

Parameters

[in]	fil_n	filter length (must be odd)
[in]	x	input abscissa
[in]	y	input ordinate
[in,out]	d2y	second derivative

Definition at line 2829 of file num_utilities.f90.

Here is the caller graph for this function:

calc_derivs_1d_id()

integer function num_utilities::calc_derivs_1d_id	(	integer, dimension(:), intent(in)	deriv,
		integer, intent(in)	dims
	)

Returns the 1D indices for derivatives of a certain order in certain dimensions.

The algorithm works by considering a structure such as the following example in three dimensions:

1. (0,0,0)
2. (1,0,0)
3. (0,1,0)
4. (0,0,1)
5. (2,0,0)
6. (1,1,0)
7. (1,0,1)
8. (0,2,0)
9. (0,1,1)
10. (0,0,2)
11. (3,0,0)
12. (2,1,0)

etc...

By then defining \(d\) as the vector of derivatives, \(n\) as the size of \(d\), \(D = \sum_i^n d(i)\) the total degree of the derivative, and \(I\) as the index of the last nonzero element in \(d\), and extending \(d\) to the left by considering \(d(0)\) = 0, the following formula can be can be deduced for the displacements in this table with respect to index 1:

displacements in this table with respect to index 1:

\[\sum_{j=0}^{I-1} \sum_{i=0}^ {D-(d(0)+...+d(j))-1} \left(\begin{array}{c}n-j+i-1\\i\end{array}\right) , \]

making use of the binomial coefficients

\[\left(\begin{array}{c}a\\b\end{array}\right) = \frac{a!}{b!(a-b)!}\]

It can be seen that each of the terms in the summation in \(j\) corresponds to the displacement in dimension \(j\) and the binomial coefficient comes from the classic stars and bars problem.

Parameters

[in]	deriv	derivatives
[in]	dims	nr. of dimensions

Definition at line 2146 of file num_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

calc_ext_var()

integer function, public num_utilities::calc_ext_var	(	real(dp), intent(inout)	ext_var,
		real(dp), dimension(:), intent(in)	var,
		real(dp), dimension(:), intent(in)	var_points,
		real(dp), intent(in)	ext_point,
		integer, intent(in), optional	deriv_in
	)

Extrapolates a function.

This is done using linear or quadratic interpolation, depending on the number of points and values given.

The data should be given sorted, in ascending order, without duplicate points in var_points.

It uses the following solution:

\[\overline{\text{A}} \ \vec{b} = \vec{c}\]

which is written out to

\[\left(\begin{array}{cccc} x_1^0 & x_1^1 & x_1^2 & \ldots \\ x_2^0 & x_2^1 & x_2^2 & \ldots \\ x_3^0 & x_3^1 & x_3^2 & \ldots \\ \vdots & \vdots & \vdots & \ddots \\ \end{array}\right) \left(\begin{array}{c} a_0 \\ a_1 \\ a_2 \\ \vdots \end{array}\right) = \left(\begin{array}{c} v_0 \\ v_1 \\ v_2 \\ \vdots \end{array}\right) \]

where \(v_{i-1}\) = var_i.

This is solved for the coefficients of the interpolating polynomial

ext_var = a_0 + a_1 x + a_2 x^2 + ...

There is an optional flag to look for the \(k^\text{th}\) derivative of the function instead of the function itself, where \(k\) should be lower than the degree of the polynomial.

Returns

ierr

Parameters

[in,out]	ext_var	output
[in]	var	abscissa
[in]	var_points	ordinate
[in]	ext_point	point to which extrapolate
[in]	deriv_in	specifies an optional derivative

Definition at line 2325 of file num_utilities.f90.

calc_inv_0d()

integer function num_utilities::calc_inv_0d	(	real(dp), dimension(:,:), intent(inout)	inv_0D,
		real(dp), dimension(:,:), intent(in)	A
	)

private constant version

Parameters

[in,out]	inv_0d	output
[in]	a	input

Definition at line 772 of file num_utilities.f90.

check_deriv()

integer function, public num_utilities::check_deriv	(	integer, dimension(3), intent(in)	deriv,
		integer, intent(in)	max_deriv,
		character(len=*), intent(in)	sr_name
	)

checks whether the derivatives requested for a certain subroutine are valid

Returns

ierr

Parameters

[in]	deriv	derivative to be checked
[in]	max_deriv	maximum derivative allowed
[in]	sr_name	name of subroutine

Definition at line 2260 of file num_utilities.f90.

Here is the caller graph for this function:

derivs()

integer function, dimension(:,:), allocatable, public num_utilities::derivs	(	integer, intent(in)	order,
		integer, intent(in), optional	dims
	)

Returns derivatives of certain order.

Parameters

[in]	order	order of derivative
[in]	dims	nr. of dimensions

Returns

array of all unique derivatives

Definition at line 2246 of file num_utilities.f90.

Here is the caller graph for this function:

fac()

recursive integer function, public num_utilities::fac

(

integer, intent(in)

n

)

Calculate factorial.

Parameters

[in]

n

input

Returns

output

Definition at line 2609 of file num_utilities.f90.

Here is the caller graph for this function:

gcd()

recursive integer function, public num_utilities::gcd	(	integer, intent(in)	u,
		integer, intent(in)	v
	)

Returns least denominator using the GCD.

See also

From https://rosettacode.org/wiki/Least_common_multiple#Fortran

Parameters

[in]	u	input
[in]	v	input

Returns

result

Definition at line 2669 of file num_utilities.f90.

Here is the caller graph for this function:

is_sym()

integer function, public num_utilities::is_sym	(	integer, intent(in)	n,
		integer, intent(in)	nn,
		logical, intent(inout)	sym
	)

Determines whether a matrix making use of the storage convention in eq_vars.eq_2_type is symmetric or not.

Returns

ierr

Parameters

[in]	n	size of matrix
[in]	nn	number of elements in matrix
[in,out]	sym	output

Definition at line 2626 of file num_utilities.f90.

Here is the caller graph for this function:

lcm()

recursive integer function, public num_utilities::lcm	(	integer, intent(in)	u,
		integer, intent(in)	v
	)

Returns common multiple using the Euclid's algorithm.

See also

From https://rosettacode.org/wiki/Greatest_common_divisor#Recursive_Euclid_algorithm_3

Parameters

[in]	u	input
[in]	v	input

Returns

result

Definition at line 2657 of file num_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

shift_f()

subroutine, public num_utilities::shift_f	(	integer, dimension(2), intent(in)	Al,
		integer, dimension(2), intent(in)	Bl,
		integer, dimension(2), intent(in)	Cl,
		real(dp), dimension(al(1):al(2),2), intent(in)	A,
		real(dp), dimension(bl(1):bl(2),2), intent(in)	B,
		real(dp), dimension(cl(1):cl(2),2), intent(inout)	C
	)

Calculate multiplication through shifting of fourier modes A and B into C.

This works by calculating \(\left[\sum_A \left(\alpha_A \cos m_A \theta + \beta_A \sin m_A \theta \right)\right] \left[\sum_B \left(\alpha_B \cos m_B \theta + \beta_B \sin m_B \theta \right)\right] = \left[\sum_C \left(\alpha_C \cos m_C \theta + \beta_C \sin m_C \theta \right)\right]\)

where the \(\alpha\) and \(\beta\) factors are provide in A, B and C.

This then boils down to finding the four combinations of cosines and sines.

Note

Modes that are larger than what C can hold are thrown away.

Parameters

[in]	al	limits on A mode numbers
[in]	bl	limits on B mode numbers
[in]	cl	limits on C mode numbers
[in]	a	input
[in]	b	input
[in,out]	c	result

Definition at line 2696 of file num_utilities.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

solve_vand()

subroutine, public num_utilities::solve_vand	(	integer, intent(in)	n,
		real(dp), dimension(n), intent(in)	a,
		real(dp), dimension(n), intent(in)	b,
		real(dp), dimension(n), intent(inout)	x,
		logical, intent(in), optional	transp
	)

Solve a Vandermonde system \(\overline{\text{A}} \ \vec{X} = \vec{B}\).

The Vandermonde matrix has the form

\[A = \left(\begin{array}{ccccc} 1 & a_1 & a_1^2 & \cdots & a_1^n \\ 1 & a_2 & a_2^2 & \cdots & a_2^n \\ 1 & a_3 & a_3^2 & \cdots & a_3^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & a_n & a_n^2 & \cdots & a_n^n \end{array}\right)\]

See also

Adapted from two routines dvand and pvand in https://people.sc.fsu.edu/~jburkardt/f_src/vandermonde/vandermonde.html based on the Björk-Pereyra Algorithm [2].

Parameters

[in]	n	matrix size
[in]	a	parameters of Vandermonde matrix
[in]	b	right-hand side
[in,out]	x	solution
[in]	transp	transposed Vandermonde matrix

Definition at line 2751 of file num_utilities.f90.

Here is the caller graph for this function:

Variable Documentation

d

integer, dimension(:,:,:), allocatable, public num_utilities::d

1-D array indices of derivatives

Definition at line 23 of file num_utilities.f90.

debug_calc_coeff_fin_diff

logical, public num_utilities::debug_calc_coeff_fin_diff = .false.

plot debug information for calc_coeff_fin_diff()

Note

Debug version only

Definition at line 28 of file num_utilities.f90.

debug_con2dis_reg

logical, public num_utilities::debug_con2dis_reg = .false.

plot debug information for con2dis_reg()

Note

Debug version only

Definition at line 27 of file num_utilities.f90.

f

integer, dimension(:,:), allocatable, public num_utilities::f

1-D array indices of Fourier mode combination indices

Definition at line 25 of file num_utilities.f90.

m

integer, dimension(:,:), allocatable, public num_utilities::m

1-D array indices of metric indices

Definition at line 24 of file num_utilities.f90.