mathematicalsubroutines Module Reference

contains either true mathematical procedures or wrappers of to the standard libraries (lapack) to assure communication with the cartesius data structures More...

Data Types
interface	real_f_one_var
	Abstract interface required to pass any real function of one variable into the numerical integration algorithm. More...

Functions/Subroutines
real(doubleprecision) function	det (rin)
	calculates the determinant of 3x3 matrix rin More...
integer function	gcd (M, N)
	calculates the greast common divisor of (m,n) More...
integer function	lcm (M, N)
	calculates the Least common factor of (m,n) More...
integer(8) function	factorial (n)
	calculates the Factorial(n) More...
real function	factorial_real (n)
	Real version of the factorial function. More...
real function	doublefactorial_real (n)
integer function	binomialcoefficient (n, m)
	calculates the BinomialCoefficient(n,m) of integers n and m More...
real function	binomialcoefficient_real (n, m)
	Real number version of the Binomial Coefficient function. More...
real function	lorentzian (w, w0, width)
	calculates the Lorentzian shape at frequency w centered at w0 with width width; integral wrt w equals to unity. More...
real function	lorentzim (w, w0, width)
	calculates the Lorentzian shape at frequency w centered at w0 with width width; integral wrt w equals to unity. \( \frac{1}{\pi} \frac{\gamma}{(\omega - \omega_0)^2 - \gamma^2} \) More...
real function	lorentzre (w, w0, width)
	calculates the real part of Lorentzian shape at frequency w centered at w0 with width width. Magnitude conforms with the imaginary part having integral of unity. More...
real function	gaussian (w, w0, width)
	calculates the real Gaussian shape at frequency w centered at w0 with width width; integral wrt w equals to unity. Magnitude conforms with the integral of unity. More...
real function	wigner3jm (j1, j2, j3, m1, m2, m3)
	calculates the Wigner 3jm symbol for integer arguments j1 j2 j3 m1 m2 m3 More...
real function	gauntcoefficient (k, l1, m1, l2, m2)
	calculates the Gaunt coefficient for intermediate momentum k for the momenta l1,l2 and projections of momenta m1,m2 More...
real function	pleg (l, m, x)
	calculates the m-th adjoint of n-th Legendre polynomial More...
real function	legendre_polynomial (l, m, x)
	General Legendre polynomial function which accepts negative m. More...
complex function	spherical_function (l, m, teta, phi)
	Calculates complex spherical function with Condon-Shortley phase. More...
real function	upperincompletegammaf (n, x)
	Calculates the "upper" incomplete gamma function defined as. More...
real function	upperincompletegammaf_real (s, x)
	Calculates the upper incomplete gamma function defined above for non-integer values of n (renamed to s) More...
real function	lowerincompletegammaf (n, x)
	Calculates the lower incomplete gamma function, the complement to the upper incomplete gamma function, for integer n. More...
real function	lowerincompletegammaf_real (s, x)
	Calculates the lower incomplete gamma function for non-integer values of n (renamed to s) More...
complex function	ex (M, S, C)
integer function	isig (N)
integer function	ihvside (n)
	Heaviside function on integers. Returns 1 if n >= 0 and 0 otherwise. More...
subroutine	sdiagv (N, A, D, X)
	MATRIX DIAGONALIZATION PROCEDURE WE USED FOR DECADES... IS GOOD DUE TO FACT THAT IT SORTS THE EIGENVALUES. More...
subroutine	simplex (start, n, EPSILON, scale, iprint, func, finish, valmin)
	gradientless minimization of function func. Our DirtyTricks are used to make the interface modern (transmission fo the function by pointer More...
subroutine	svdcmp (a, m, n, mp, np, w, v)
	Singular value decomposition of matrix A. More...
subroutine	gradient_minimization (pFunction, pFirstDerivative, init_variables, step, prec, reset)
	Function allows to find local minimum of function by using gradient optimization methid. More...
real function	adapt_simpson_int (f, x_min, x_max, iprec, pars)
	Performs adaptive Simpson's numerical integration of a real function of one variable. This algorithm finds the optimal number of points required to reach the given accuracy of integration, which usually allows to decrease number of function evaluations compared to standard integration algorithm. More...
real function, dimension(2)	optimize_f_goldensection (f, xmin, xmax, prec, pars)
	Finds minimum of a real function of one variable in a given interval using golden-section search method. Can be used for any unimodal function, when the gradients are not known. Returns an interval where the minimum is located with required precision. More...
real function	opt_scalar_brent (f, xmin, xmax, tol, stat, pars)
	Finds minimum of scalar function f bracketed by [xmin, xmax] More...
real function	solve_scalar_brent (f, xmin, xmax, tol, stat, pars)
	Finds simple root of scalar function f bracketed by [xmin, xmax] More...
recursive subroutine	generalizedlaguerrepolynomialcoefficients (n, l, coef_array)
subroutine	getjacobirotationmatrix (space_dim, dim1, dim2, angle, JacobiMatrix)
	produces Jacobi rotation matrix for angle angle of the size space_dim x space_dim for a pair of components with numbers dim1, dim2 More...
subroutine	getfirstderivativeofjacobirotationmatrix (space_dim, dim1, dim2, angle, DerivativeMatrix)
	produces first derivative wrt angle of the Jacobi rotation matrix for angle angle of the size space_dim x space_dim for a pair of components with numbers dim1, dim2 More...
subroutine	getsecondderivativeofjacobirotationmatrix (space_dim, dim1, dim2, angle, DerivativeMatrix)
	produces second derivative wrt angle of the Jacobi rotation matrix for angle angle of the size space_dim x space_dim for a pair of components with numbers dim1, dim2 More...
subroutine	rotate_matrix_left (space_dim, dim1, dim2, angle, mat)
	Multiplies the given matrix by a Jacobi rotation matrix from the left. More...
subroutine	rotate_matrix_right (space_dim, dim1, dim2, angle, mat)
	Multiplies the given matrix by a Jacobi rotation matrix from the left. More...
subroutine	jacobi (n, a, d, v)
	Diagonalization of matrix a(n,n) by Jacobi method n is the matrix dimension a(n,n) - matrix d(n) - eigenvalues array v(n,n) – eigenvectors array. More...
real function, dimension(n, n)	adj_singular (n, M)
	Calculates adjugate matrix of the degenerate NxN matrix M (detM = 0). More...

Variables
logical, private	ifreach
	Private variable required for the adaptive numerical integration: More...

Detailed Description

contains either true mathematical procedures or wrappers of to the standard libraries (lapack) to assure communication with the cartesius data structures

Function/Subroutine Documentation

adapt_simpson_int()

real function mathematicalsubroutines::adapt_simpson_int	(	procedure(real_f_one_var), intent(in), pointer	f,
		real, intent(in)	x_min,
		real, intent(in)	x_max,
		real, intent(in)	iprec,
		class(*), intent(inout), optional, pointer	pars
	)

Performs adaptive Simpson's numerical integration of a real function of one variable. This algorithm finds the optimal number of points required to reach the given accuracy of integration, which usually allows to decrease number of function evaluations compared to standard integration algorithm.

Parameters

f	- abstract integrand;
pars	- optional polymorphic variable required to pass additional parameters into "f";
x_min	- left integration limit;
x_max	- right integration limit;
iprec	- required precision of numerical integration.

Here is the call graph for this function:

adj_singular()

real function, dimension(n,n) mathematicalsubroutines::adj_singular	(	integer, intent(in)	n,
		real, dimension(n,n), intent(in)	M
	)

Calculates adjugate matrix of the degenerate NxN matrix M (detM = 0).

Note

Here we can not use the usual algorithm where inverse matrix is calculated and then multiplied by determinant. Therefore we use less efficient method.

Warning

This algorithm is not the most efficient one, because it scales as O((1/3)*Nˆ5). There are more involved approaches for computing adjugates of degenerate matrices. For example: https://core.ac.uk/download/pdf/82551348.pdf

@TODO: to code more efficient algorithm.

binomialcoefficient()

integer function mathematicalsubroutines::binomialcoefficient	(	integer, intent(in)	n,
		integer, intent(in)	m
	)

calculates the BinomialCoefficient(n,m) of integers n and m

Here is the call graph for this function:

binomialcoefficient_real()

real function mathematicalsubroutines::binomialcoefficient_real	(	integer, intent(in)	n,
		integer, intent(in)	m
	)

Real number version of the Binomial Coefficient function.

Here is the call graph for this function:

det()

real(doubleprecision) function mathematicalsubroutines::det

(

real(doubleprecision), dimension(1:3,1:3), intent(in)

rin

)

calculates the determinant of 3x3 matrix rin

doublefactorial_real()

real function mathematicalsubroutines::doublefactorial_real

(

integer, intent(in)

n

)

ex()

complex function mathematicalsubroutines::ex	(	integer, intent(in)	M,
		real, intent(in)	S,
		real, intent(in)	C
	)

Here is the call graph for this function:

factorial()

integer(8) function mathematicalsubroutines::factorial

(

integer, intent(in)

n

)

calculates the Factorial(n)

Attention

Changed output kind from int4 to int8 to allow correct evaluation of factorials for "n" higher than 12. - Ilya Popov

factorial_real()

real function mathematicalsubroutines::factorial_real

(

integer, intent(in)

n

)

Real version of the factorial function.

gauntcoefficient()

real function mathematicalsubroutines::gauntcoefficient	(	integer, intent(in)	k,
		integer, intent(in)	l1,
		integer, intent(in)	m1,
		integer, intent(in)	l2,
		integer, intent(in)	m2
	)

calculates the Gaunt coefficient for intermediate momentum k for the momenta l1,l2 and projections of momenta m1,m2

Here is the call graph for this function:

gaussian()

real function mathematicalsubroutines::gaussian	(	real, intent(in)	w,
		real, intent(in)	w0,
		real, intent(in)	width
	)

calculates the real Gaussian shape at frequency w centered at w0 with width width; integral wrt w equals to unity. Magnitude conforms with the integral of unity.

Here is the call graph for this function:

gcd()

integer function mathematicalsubroutines::gcd	(	integer, intent(in)	M,
		integer, intent(in)	N
	)

calculates the greast common divisor of (m,n)

Here is the call graph for this function:

generalizedlaguerrepolynomialcoefficients()

recursive subroutine mathematicalsubroutines::generalizedlaguerrepolynomialcoefficients	(	integer, intent(in)	n,
		integer, intent(in)	l,
		real, dimension(:,:), intent(inout), allocatable	coef_array
	)

Here is the call graph for this function:

getfirstderivativeofjacobirotationmatrix()

subroutine mathematicalsubroutines::getfirstderivativeofjacobirotationmatrix	(	integer, intent(in)	space_dim,
		integer, intent(in)	dim1,
		integer, intent(in)	dim2,
		real, intent(in)	angle,
		real, dimension(:,:), intent(out), allocatable	DerivativeMatrix
	)

produces first derivative wrt angle of the Jacobi rotation matrix for angle angle of the size space_dim x space_dim for a pair of components with numbers dim1, dim2

getjacobirotationmatrix()

subroutine mathematicalsubroutines::getjacobirotationmatrix	(	integer, intent(in)	space_dim,
		integer, intent(in)	dim1,
		integer, intent(in)	dim2,
		real, intent(in)	angle,
		real, dimension(:,:), intent(out), allocatable	JacobiMatrix
	)

produces Jacobi rotation matrix for angle angle of the size space_dim x space_dim for a pair of components with numbers dim1, dim2

getsecondderivativeofjacobirotationmatrix()

subroutine mathematicalsubroutines::getsecondderivativeofjacobirotationmatrix	(	integer, intent(in)	space_dim,
		integer, intent(in)	dim1,
		integer, intent(in)	dim2,
		real, intent(in)	angle,
		real, dimension(:,:), intent(out), allocatable	DerivativeMatrix
	)

produces second derivative wrt angle of the Jacobi rotation matrix for angle angle of the size space_dim x space_dim for a pair of components with numbers dim1, dim2

gradient_minimization()

subroutine mathematicalsubroutines::gradient_minimization	(	procedure(realfunctionofrealptraray), intent(inout), pointer	pFunction,
		procedure(derivativeoffunctionofrealptraray), intent(inout), pointer	pFirstDerivative,
		type(real_ptr), dimension(:), intent(inout)	init_variables,
		real, intent(in)	step,
		real, intent(in)	prec,
		procedure(realroutineofrealptraray), intent(inout), optional, pointer	reset
	)

Function allows to find local minimum of function by using gradient optimization methid.

ihvside()

integer function mathematicalsubroutines::ihvside

(

integer, intent(in)

n

)

Heaviside function on integers. Returns 1 if n >= 0 and 0 otherwise.

Here is the call graph for this function:

isig()

integer function mathematicalsubroutines::isig

(

integer, intent(in)

N

)

Here is the call graph for this function:

jacobi()

subroutine mathematicalsubroutines::jacobi	(	integer	n,
		real, dimension(n,n)	a,
		real, dimension(n)	d,
		real, dimension(n,n)	v
	)

Diagonalization of matrix a(n,n) by Jacobi method n is the matrix dimension a(n,n) - matrix d(n) - eigenvalues array v(n,n) – eigenvectors array.

Remarks

- repeats functionality of sdiagv

lcm()

integer function mathematicalsubroutines::lcm	(	integer, intent(in)	M,
		integer, intent(in)	N
	)

calculates the Least common factor of (m,n)

Here is the call graph for this function:

legendre_polynomial()

real function mathematicalsubroutines::legendre_polynomial	(	integer, intent(in)	l,
		integer, intent(in)	m,
		real, intent(in)	x
	)

General Legendre polynomial function which accepts negative m.

Here is the call graph for this function:

lorentzian()

real function mathematicalsubroutines::lorentzian	(	real, intent(in)	w,
		real, intent(in)	w0,
		real, intent(in)	width
	)

calculates the Lorentzian shape at frequency w centered at w0 with width width; integral wrt w equals to unity.

Todo:

to be replaced by LorentzIm throughout all applications

Here is the call graph for this function:

lorentzim()

real function mathematicalsubroutines::lorentzim	(	real, intent(in)	w,
		real, intent(in)	w0,
		real, intent(in)	width
	)

calculates the Lorentzian shape at frequency w centered at w0 with width width; integral wrt w equals to unity. \( \frac{1}{\pi} \frac{\gamma}{(\omega - \omega_0)^2 - \gamma^2} \)

Here is the call graph for this function:

lorentzre()

real function mathematicalsubroutines::lorentzre	(	real, intent(in)	w,
		real, intent(in)	w0,
		real, intent(in)	width
	)

calculates the real part of Lorentzian shape at frequency w centered at w0 with width width. Magnitude conforms with the imaginary part having integral of unity.

Here is the call graph for this function:

lowerincompletegammaf()

real function mathematicalsubroutines::lowerincompletegammaf	(	integer, intent(in)	n,
		real, intent(in)	x
	)

Calculates the lower incomplete gamma function, the complement to the upper incomplete gamma function, for integer n.

\[ \gamma\left(n,x\right)=\Gamma\left(n\right) - \Gamma\left(n,x\right) \]

Here is the call graph for this function:

lowerincompletegammaf_real()

real function mathematicalsubroutines::lowerincompletegammaf_real	(	real, intent(in)	s,
		real, intent(in)	x
	)

Calculates the lower incomplete gamma function for non-integer values of n (renamed to s)

Here is the call graph for this function:

opt_scalar_brent()

real function mathematicalsubroutines::opt_scalar_brent	(	procedure(real_f_one_var), intent(in), pointer	f,
		real, intent(in)	xmin,
		real, intent(in)	xmax,
		real, intent(in), optional	tol,
		integer, intent(out), optional	stat,
		class(*), intent(inout), optional, pointer	pars
	)

Finds minimum of scalar function f bracketed by [xmin, xmax]

Parameters

f	Function to be minimized
xmin	Left bracket
xmax	Right bracket
[in,optional]	default = 32 * epsilon] tol The accuracy of the minimum
[out,optional]	stat Status of solver. 0 - success, 1 - did not converge
[inout,optional]	pars Additional parameters for the function f

Note

Adapted from Numerical Recipes.

optimize_f_goldensection()

real function, dimension(2) mathematicalsubroutines::optimize_f_goldensection	(	procedure(real_f_one_var), intent(in), pointer	f,
		real, intent(in)	xmin,
		real, intent(in)	xmax,
		real, intent(in)	prec,
		class(*), intent(inout), optional, pointer	pars
	)

Finds minimum of a real function of one variable in a given interval using golden-section search method. Can be used for any unimodal function, when the gradients are not known. Returns an interval where the minimum is located with required precision.

Parameters

f	- function of interest;
pars	- optional polymorphic variable required to pass additional parameters into "f";
x_min	- left boundary of the interval;
x_max	- right boundary of the interval;
prec	- required precision.

pleg()

real function mathematicalsubroutines::pleg	(	integer, intent(in)	l,
		integer, intent(in)	m,
		real, intent(in)	x
	)

calculates the m-th adjoint of n-th Legendre polynomial

Note

There exist two definitions of these polynomials in the literature - with and without the Condon-Shortley factor (-1)^m. See corresponding discussions via the following links: https://mathworld.wolfram.com/AssociatedLegendrePolynomial.html https://mathworld.wolfram.com/Condon-ShortleyPhase.html In Wolfram Mathematica Condon-Shortley phase is included into "LegendreP". Moreover an equation of spherical harmonics from the Varshalovich's book (eq.1 paragraph 5.2) implies that this phase is included into the Legendre polynomial. Therefore, here we include the Condon-Shortley phase into this function.

Tests against Wolfram Mathematica benchmarks are available in "cartesius_tests" repository.

Here is the call graph for this function:

rotate_matrix_left()

subroutine mathematicalsubroutines::rotate_matrix_left	(	integer, intent(in)	space_dim,
		integer, intent(in)	dim1,
		integer, intent(in)	dim2,
		real, intent(in)	angle,
		real, dimension(space_dim, space_dim), intent(inout)	mat
	)

Multiplies the given matrix by a Jacobi rotation matrix from the left.

Parameters

[in]	space_dim	Dimensions of the given square matrix mat
[in]	dim1,dim2	Dimensions of the rotation plane. Are sorted to be in ascending order.
[in]	angle	The rotation angle
[in,out]	mat	The matrix to be rotated. The matrix chenges in-place.

rotate_matrix_right()

subroutine mathematicalsubroutines::rotate_matrix_right	(	integer, intent(in)	space_dim,
		integer, intent(in)	dim1,
		integer, intent(in)	dim2,
		real, intent(in)	angle,
		real, dimension(space_dim, space_dim), intent(inout)	mat
	)

Multiplies the given matrix by a Jacobi rotation matrix from the left.

Parameters

[in]	space_dim	Dimensions of the given square matrix mat
[in]	dim1,dim2	Dimensions of the rotation plane. Are sorted to be in ascending order.
[in]	angle	The rotation angle
[in,out]	mat	The matrix to be rotated. The matrix chenges in-place.

sdiagv()

subroutine mathematicalsubroutines::sdiagv	(	integer, intent(in)	N,
		real, dimension(:,:), intent(in)	A,
		real, dimension(:), intent(out)	D,
		real, dimension(:,:), intent(out)	X
	)

MATRIX DIAGONALIZATION PROCEDURE WE USED FOR DECADES... IS GOOD DUE TO FACT THAT IT SORTS THE EIGENVALUES.

simplex()

subroutine mathematicalsubroutines::simplex	(	real, dimension(0:n-1), intent(in)	start,
		integer, intent(in)	n,
		real, intent(in)	EPSILON,
		real, intent(in)	scale,
		integer, intent(in)	iprint,
		procedure (realfunctionofrealarray), pointer	func,
		real, dimension(0:n-1), intent(out)	finish,
		real, intent(out)	valmin
	)

gradientless minimization of function func. Our DirtyTricks are used to make the interface modern (transmission fo the function by pointer

Here is the call graph for this function:

solve_scalar_brent()

real function mathematicalsubroutines::solve_scalar_brent	(	procedure(real_f_one_var), intent(in), pointer	f,
		real, intent(in)	xmin,
		real, intent(in)	xmax,
		real, intent(in), optional	tol,
		integer, intent(out), optional	stat,
		class(*), intent(inout), optional, pointer	pars
	)

Finds simple root of scalar function f bracketed by [xmin, xmax]

Parameters

f	Function to be solved
xmin	Left bracket
xmax	Right bracket
[in,optional]	default = 32 * epsilon] tol The accuracy of the root
[out,optional]	stat Status of solver. 0 - success, 1 - did not converge, 2 - incorrect input
[inout,optional]	pars Addition parameters for the function f

Note

Adapted from Numerical Recipes.

spherical_function()

complex function mathematicalsubroutines::spherical_function	(	integer, intent(in)	l,
		integer, intent(in)	m,
		real, intent(in)	teta,
		real, intent(in)	phi
	)

Calculates complex spherical function with Condon-Shortley phase.

Note

Equation from the Varshalovich's book (eq.1 paragraph 5.2).

For negative m the associated Legendre polynomial is defined by

\[ P_{l}^{-\left|m\right|}\left(x\right)=\left(-1\right)^{\left|m\right|} \frac{\left(l-\left|m\right|\right)!}{\left(l+\left|m\right|\right)!} P_{l}^{\left|m\right|}\left(x\right) \]

Note

Tests against Wolfram Mathematica benchmarks are available in "cartesius_tests" repository.

Here is the call graph for this function:

svdcmp()

subroutine mathematicalsubroutines::svdcmp	(	real(doubleprecision), dimension(mp,np)	a,
		integer	m,
		integer	n,
		integer	mp,
		integer	np,
		real(doubleprecision), dimension(np)	w,
		real(doubleprecision), dimension(np,np)	v
	)

Singular value decomposition of matrix A.

upperincompletegammaf()

real function mathematicalsubroutines::upperincompletegammaf	(	integer, intent(in)	n,
		real, intent(in)	x
	)

Calculates the "upper" incomplete gamma function defined as.

\[ \Gamma\left(n,x\right)=\int_{x}^{\infty}t^{n-1}\exp\left(-t\right)dt \]

Here is the call graph for this function:

upperincompletegammaf_real()

real function mathematicalsubroutines::upperincompletegammaf_real	(	real, intent(in)	s,
		real, intent(in)	x
	)

Calculates the upper incomplete gamma function defined above for non-integer values of n (renamed to s)

Note

Uses continued fraction approximation:

\[ \Gamma\left(n,x\right)=\frac{e^{-x}x^{s}}{x+1-s-}\frac{1\cdot\left(1-s\right)}{x+3-s-}\frac{2\cdot\left(2-s\right)}{x+5-s-}.. \]

Here is the call graph for this function:

wigner3jm()

real function mathematicalsubroutines::wigner3jm	(	integer, intent(in)	j1,
		integer, intent(in)	j2,
		integer, intent(in)	j3,
		integer, intent(in)	m1,
		integer, intent(in)	m2,
		integer, intent(in)	m3
	)

calculates the Wigner 3jm symbol for integer arguments j1 j2 j3 m1 m2 m3

Author

A. Tchougreeff ca. 1981

D. Raenko (updated summer 2022)

Note

Function works properly now for a large number of values. Tested against Wolfram Mathematica and GNU Scientific Library function. It is also 3 times faster now.

Here is the call graph for this function:

Variable Documentation

ifreach

logical, private mathematicalsubroutines::ifreach

private

Private variable required for the adaptive numerical integration: