Namespaces
	transforms

Classes
struct	floating_point_limits
	Traits class for floating point types to query limit information. The floating_point_limits class is traits class providing a standardized and portable way to query the following information on either host/device: More...

struct	floating_point_limits< float >

struct	floating_point_limits< double >

struct	floating_point_limits< long double >

class	Matrix
	The Matrix class is used to represent \( M \times N \) matrices. It provides common matrix operations and allows accessing matrix elements in a more natural way, using row and column indices, regardless of the underlying flat array storage layout. More...

struct	PolynomialRootResult
	Result metadata for an all-roots polynomial solve. More...

class	QuadratureRule
	Stores fixed views to arrays of 1D quadrature points and weights. More...

Enumerations
enum	ReturnCodes { LU_SUCCESS , LU_SINGULAR_MATRIX , LU_NONSQUARE_MATRIX }

Enumerators & Types
enum	MatrixNorm { P1_NORM , INF_NORM , FROBENIUS_NORM }
	Enumerates the supported matrix norms. More...

Functions
template<typename T >
int	eigen_solve (Matrix< T > &A, int k, T u, T lambdas, int numIterations=160)
	Approximates k eigenvectors and eigenvalues of the passed in square matrix using the power method. More...

template<typename T >
bool	eigen_sort (T *lamdas, Matrix< T > &eigen_vectors)
	Sorts the supplied eigenvalues and eigenvectors in ascending order. More...

template<typename T >
int	jacobi_eigensolve (Matrix< T > A, Matrix< T > &V, T lambdas, int maxIterations=JACOBI_DEFAULT_MAX_ITERATIONS, int numIterations=nullptr, T TOL=JACOBI_DEFAULT_TOLERANCE)
	Computes the eigenvalues and eigenvectors of a real symmetric matrix using the Jacobi iteration method. More...

template<typename T >
int	linear_solve (Matrix< T > &A, const T b, T x)
	Solves a linear system of the form \( Ax=b \). More...

template<typename ScalarType >
ScalarType	evaluate_polynomial (ArrayView< const double > coeffs_descending, const ScalarType &x)
	Evaluate a polynomial with coefficients stored in descending power order. More...

axom::Array< double >	bernstein_to_monomial (ArrayView< const double > bernstein_coeffs)
	Convert Bernstein coefficients to monomial coefficients on `[0, 1]`. More...

int	effective_polynomial_degree (ArrayView< const double > coeffs_ascending, double tol=1e-12)
	Return the effective degree of a polynomial after trimming near-zero leading coefficients. More...

axom::Array< std::complex< double > >	solve_polynomial_durand_kerner (ArrayView< const double > coeffs_ascending, double tol=1e-12, std::complex< double > seed=std::complex< double > {0.4, 0.9})
	Approximate all roots of a polynomial using the Durand-Kerner method. More...

PolynomialRootResult	solve_polynomial_durand_kerner_checked (ArrayView< const double > coeffs_ascending, double tol=1e-12, int max_iters=200, std::complex< double > seed=std::complex< double > {0.4, 0.9})
	Approximate all roots of a polynomial using the Durand-Kerner method, including convergence diagnostics. More...

int	solve_linear (ArrayView< const double > coeff, ArrayView< double > roots, int &numRoots)
	Find the real root for a linear equation of form \( ax + b = 0 \). More...

int	solve_linear (const double coeff, double roots, int &numRoots)
	Deprecated pointer-based overload for solve_linear(ArrayView<const double>, ArrayView<double>, int&). More...

int	solve_quadratic (ArrayView< const double > coeff, ArrayView< double > roots, int &numRoots)
	For a quadratic equation of the form \( ax^2 + bx + c = 0 \), find real roots using the quadratic formula. More...

int	solve_quadratic (const double coeff, double roots, int &numRoots)
	Deprecated pointer-based overload for solve_quadratic(ArrayView<const double>, ArrayView<double>, int&). More...

int	solve_cubic (ArrayView< const double > coeff, ArrayView< double > roots, int &numRoots)
	For a cubic equation of the form \( ax^3 + bx^2 + cx + d = 0 \), find real roots. More...

int	solve_cubic (const double coeff, double roots, int &numRoots)
	Deprecated pointer-based overload for solve_cubic(ArrayView<const double>, ArrayView<double>, int&). More...

void	compute_gauss_legendre_data (int npts, axom::Array< double > &nodes, axom::Array< double > &weights, int allocatorID=axom::getDefaultAllocatorID())
	Computes a 1D quadrature rule of Gauss-Legendre points. More...

QuadratureRule	get_gauss_legendre (int npts, int allocatorID=axom::getDefaultAllocatorID())
	Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points. More...

Matrix Operators
template<typename real >
AXOM_HOST_DEVICE real	determinant (const real &a00, const real &a01, const real &a10, const real &a11)
	Computes a 2X2 determinant of the given matrix. More...

template<typename real >
AXOM_HOST_DEVICE real	determinant (const real &a00, const real &a01, const real &a02, const real &a10, const real &a11, const real &a12, const real &a20, const real &a21, const real &a22)
	Computes the 3x3 determinant of the given matrix. More...

template<typename real >
AXOM_HOST_DEVICE real	determinant (const real &a00, const real &a01, const real &a02, const real &a03, const real &a10, const real &a11, const real &a12, const real &a13, const real &a20, const real &a21, const real &a22, const real &a23, const real &a30, const real &a31, const real &a32, const real &a33)
	Computes the 4x4 determinant of the given matrix. More...

template<typename real >
real	determinant (const Matrix< real > &A)
	Computes the determinant of the given square matrix. More...

template<typename T >
int	lu_decompose (Matrix< T > &A, int *pivots)
	Perform LU-decomposition on a given square matrix, \( A \) using partial pivoting. More...

template<typename T >
int	lu_solve (const Matrix< T > &A, const int pivots, const T b, T *x)
	Solve the system \( Ax=b \) by back-substitution, where, A is an LU decomposed matrix produced via a call to lu_decompose(). More...

template<typename T >
Matrix< T >	lower_triangular (const Matrix< T > &A, bool unit_diagonal=false)
	Extracts the lower triangular part of a square matrix. More...

template<typename T >
Matrix< T >	upper_triangular (const Matrix< T > &A, bool unit_diagonal=true)
	Extract the upper triangular part of a square matrix. More...

template<typename T >
bool	matrix_add (const Matrix< T > &A, const Matrix< T > &B, Matrix< T > &C)
	Computes \( \mathcal{C} = \mathcal{A} + \mathcal{B} \). More...

template<typename T >
bool	matrix_subtract (const Matrix< T > &A, const Matrix< T > &B, Matrix< T > &C)
	Computes \( \mathcal{C} = \mathcal{A} - \mathcal{B} \). More...

template<typename T >
bool	matrix_multiply (const Matrix< T > &A, const Matrix< T > &B, Matrix< T > &C)
	Computes the matrix-matrix product of \( \mathcal{A} \) and \( \mathcal{B} \) and stores the result in \( \mathcal{C} \). More...

template<typename T >
void	matrix_scalar_multiply (Matrix< T > &A, const T &c)
	Computes a scalar-matrix produect of a matrix \( \mathcal{A} \) and a scalar \( c \) and stores the result in \( \mathcal{A} \). More...

template<typename T >
void	matrix_vector_multiply (const Matrix< T > &A, const T vec, T output)
	Computes the matrix-vector product of a matrix \(\mathcal{A}\) and a vector \(\mathbf{x}\) and store the result in the user-supplied output vector. More...

template<typename T >
bool	matrix_transpose (const Matrix< T > &A, Matrix< T > &M)
	Computes the matrix transpose of a given matrix \( \mathcal{A} \). More...

template<typename T >
T	matrix_norm (const Matrix< T > &A, MatrixNorm normType)
	Computes the specified norm of a given \( M \times N \) matrix. More...

Overloaded Matrix Operators
template<typename T >
std::ostream &	operator<< (std::ostream &os, const Matrix< T > &A)
	Overloaded output stream operator. Outputs the matrix coefficients in to the given output stream. More...

template<typename T >
bool	operator== (const Matrix< T > &lhs, const Matrix< T > &rhs)
	Checks if two matrices are component-wise equal. More...

template<typename T >
bool	operator!= (const Matrix< T > &lhs, const Matrix< T > &rhs)
	Checks if two matrices are not component-wise equal. More...

Vector Operators
template<typename T >
bool	linspace (const T &x0, const T &x1, T *v, int N)
	Generates a vector consisting of a sequence of N uniformly spaced values over the interval [x0,x1]. More...

template<typename T >
AXOM_HOST_DEVICE void	cross_product (const T u, const T v, T *w)
	Computes the vector cross-product of two vectors, u and v. More...

template<typename T >
AXOM_HOST_DEVICE T	dot_product (const T u, const T v, int dim)
	Computes the dot product of the arrays u and v. More...

template<typename T >
void	make_orthogonal (T u, T v, int dim, double tol=1E-16)
	Makes u orthogonal to v. More...

template<typename T >
bool	orthonormalize (T *basis, int size, int dim, double eps=1E-16)
	Performs Gram-Schmidt orthonormalization in-place on a 2D array of shape size,dim where it treats the rows as the individual vectors. More...

template<typename T >
AXOM_HOST_DEVICE bool	normalize (T *v, int dim, double eps=1e-16)
	Normalizes the passed in array. More...

Variables
constexpr double	JACOBI_DEFAULT_TOLERANCE = 1.e-18

constexpr int	JACOBI_DEFAULT_MAX_ITERATIONS = 20

constexpr int	JACOBI_EIGENSOLVE_SUCCESS = 0

constexpr int	JACOBI_EIGENSOLVE_FAILURE = -1

Enumeration Type Documentation

◆ ReturnCodes

enum axom::numerics::ReturnCodes

Enumerator
LU_SUCCESS
LU_SINGULAR_MATRIX
LU_NONSQUARE_MATRIX

◆ MatrixNorm

enum axom::numerics::MatrixNorm

Enumerates the supported matrix norms.

Enumerator
P1_NORM	P1_NORM.
INF_NORM	INF_NORM.
FROBENIUS_NORM	FROBENIUS_NORM.

Function Documentation

◆ determinant() [1/4]

template<typename real >

AXOM_HOST_DEVICE real axom::numerics::determinant	(	const real &	a00,
		const real &	a01,
		const real &	a10,
		const real &	a11
	)

inline

Computes a 2X2 determinant of the given matrix.

Parameters

[in]	a00	matrix element
[in]	a01	matrix element
[in]	a10	matrix element
[in]	a11	matrix element

Returns: det the determinant of the 2X2 matrix

References axom::utilities::fma().

◆ determinant() [2/4]

template<typename real >

AXOM_HOST_DEVICE real axom::numerics::determinant	(	const real &	a00,
		const real &	a01,
		const real &	a02,
		const real &	a10,
		const real &	a11,
		const real &	a12,
		const real &	a20,
		const real &	a21,
		const real &	a22
	)

inline

Computes the 3x3 determinant of the given matrix.

Parameters

[in]	a00	matrix element
[in]	a01	matrix element
[in]	a02	matrix element
[in]	a10	matrix element
[in]	a11	matrix element
[in]	a12	matrix element
[in]	a20	matrix element
[in]	a21	matrix element
[in]	a22	matrix element

Returns: det the determinant of the 3x3 matrix.

References axom::utilities::fma().

◆ determinant() [3/4]

template<typename real >

AXOM_HOST_DEVICE real axom::numerics::determinant	(	const real &	a00,
		const real &	a01,
		const real &	a02,
		const real &	a03,
		const real &	a10,
		const real &	a11,
		const real &	a12,
		const real &	a13,
		const real &	a20,
		const real &	a21,
		const real &	a22,
		const real &	a23,
		const real &	a30,
		const real &	a31,
		const real &	a32,
		const real &	a33
	)

inline

Computes the 4x4 determinant of the given matrix.

Parameters

[in]	a00	matrix element
[in]	a01	matrix element
[in]	a02	matrix element
[in]	a03	matrix element
[in]	a10	matrix element
[in]	a11	matrix element
[in]	a12	matrix element
[in]	a13	matrix element
[in]	a20	matrix element
[in]	a21	matrix element
[in]	a22	matrix element
[in]	a23	matrix element
[in]	a30	matrix element
[in]	a31	matrix element
[in]	a32	matrix element
[in]	a33	matrix element

Returns: det the determinant of the 4x4 matrix

2x2 minors

3x3 minors

det = m123*a03 - m023*a13 + m013*a23 - m012*a33

References axom::utilities::fma().

◆ determinant() [4/4]

template<typename real >

real axom::numerics::determinant ( const Matrix< real > & A )

Computes the determinant of the given square matrix.

Parameters

[in] A an \( N \times N \) input matrix

Returns: det the computed determinant.

Note: if \( A \) is not square or empty, this function will return 0.0

References determinant(), axom::numerics::Matrix< T >::empty(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::isSquare(), lu_decompose(), and LU_SUCCESS.

◆ eigen_solve()

template<typename T >

int axom::numerics::eigen_solve	(	Matrix< T > &	A,
		int	k,
		T *	u,
		T *	lambdas,
		int	numIterations = `160`
	)

Approximates k eigenvectors and eigenvalues of the passed in square matrix using the power method.

We compute eigenvectors and eigenvalues from the power method as follows:

Generate a random vector
Make orthonormal to previous eigenvectors if any
Run depth iterations of power method
Store the result

Because step 1 and step 3 require random numbers, this function requires the caller to properly seed the randomizer (for example, by calling srand()).

Parameters

[in]	A	a square input matrix
[in]	k	number of eigenvalue-eigenvectors to find
[out]	u	pointer to k eigenvectors in order by magnitude of eigenvalue
[out]	lambdas	pointer to k eigenvalues in order by size
[in]	numIterations	optional number of iterations for the power method

Note: if k <= 0, the solve is declared successful

Returns: rc return value, nonzero if the solve is successful.

Precondition: A.isSquare() == true; u != nullptr; lambdas != nullptr; k <= A.getNumRows(); T is a floating point type; System randomizer has been initialized (for example, with srand())

References AXOM_STATIC_ASSERT_MSG, dot_product(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), and matrix_vector_multiply().

◆ eigen_sort()

template<typename T >

bool axom::numerics::eigen_sort	(	T *	lamdas,
		Matrix< T > &	eigen_vectors
	)

Sorts the supplied eigenvalues and eigenvectors in ascending order.

Parameters

[in,out]	lambdas	pointer to buffer consisting of the eigenvalues
[in,out]	eigen_vectors	matrix of the corresponding eigenvectors

Note: The eigen_vectors matrix is an orthogonal matrix consisting of the eigenvectors. Each column in the matrix stores an eigenvector that is associated with the an eigenvalue in the corresponding entry in the supplied lambdas array.; lambdas should point to a buffer that holds eigen_vectors.getNumCols(); The sorting algorithm in the current implementation is O( \( N^2 \)) rather than O( \( NlogN \) )

Precondition: lambdas != nullptr; eigen_vectors.getNumRows() >= 1; eigen_vectors.getNumCols() >= 1

References axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), axom::utilities::swap(), and axom::numerics::Matrix< T >::swapColumns().

◆ jacobi_eigensolve()

template<typename T >

int axom::numerics::jacobi_eigensolve	(	Matrix< T >	A,
		Matrix< T > &	V,
		T *	lambdas,
		int	maxIterations = `JACOBI_DEFAULT_MAX_ITERATIONS`,
		int *	numIterations = `nullptr`,
		T	TOL = `JACOBI_DEFAULT_TOLERANCE`
	)

Computes the eigenvalues and eigenvectors of a real symmetric matrix using the Jacobi iteration method.

Parameters

[in]	A	the input matrix whose eigenpairs will be computed
[out]	V	matrix to store the eigenvectors of the input matrix
[out]	lambdas	buffer where the computed eigenvalues will be stored.
[in]	maxIterations	the maximum number of iterations (optional)
[out]	numIterations	the number of actual iterations (optional)
[in]	TOL	convergence tolerance. Default is set to 1.e-18 (optional)

Returns: status returns JACOBI_EIGENSOLVE_SUCCESS on success, otherwise, JACOBI_EIGENSOLVE_FAILURE is returned.

Note: A return status of JACOBI_EIGENSOLVE_FAILURE, can indicate: ( a ) a problem with the supplied input, e.g., nullptr etc. ( b ) the method did not converge in the specified number of max iterations.; The supplied matrix is expected to be a real symmetric matrix.; The Jacobi method is absolutely foolproof for all real symmetric matrices. It is efficient for moderate size matrices and generally can evaluate the smaller eigenvalues with better relative accuracy than other methods that convert the matrix to tridiagonal form, e.g., QR. However, for matrices of order greater than, e.g., 20, the Jacobi algorithm is generally slower by a significant constant factor.; The Jacobi iteration will generally converge within 6-15 iterations. The default max number of iterations is set to 20, but, the caller may specify a different value if necessary.; The implementation of this method is based on the algorithm described in Numerical Recipes, Chapter 11.1, "Jacobi Transformations of a Symmetric Matrix", p. 570.

Precondition: A.isSquare() == true; V.isSquare() == true; lambdas != nullptr

References axom::utilities::abs(), AXOM_STATIC_ASSERT_MSG, axom::deallocate(), eigen_sort(), axom::numerics::Matrix< T >::getNumRows(), axom::utilities::isNearlyEqual(), axom::numerics::Matrix< T >::isSquare(), JACOBI_EIGENSOLVE_FAILURE, and JACOBI_EIGENSOLVE_SUCCESS.

◆ linear_solve()

template<typename T >

int axom::numerics::linear_solve	(	Matrix< T > &	A,
		const T *	b,
		T *	x
	)

Solves a linear system of the form \( Ax=b \).

Parameters

[in]	A	a square input matrix
[in]	b	the right-hand side
[out]	x	the solution vector (computed)

Returns: rc return value, 0 if the solve is successful.

Precondition: A.isSquare() == true; b != nullptr; x != nullptr

Note: The input matrix is destroyed (modified) in the process.

References determinant(), axom::numerics::Matrix< T >::getNumColumns(), axom::utilities::isNearlyEqual(), axom::numerics::Matrix< T >::isSquare(), lu_decompose(), LU_NONSQUARE_MATRIX, lu_solve(), and LU_SUCCESS.

◆ lu_decompose()

template<typename T >

int axom::numerics::lu_decompose	(	Matrix< T > &	A,
		int *	pivots
	)

Perform LU-decomposition on a given square matrix, \( A \) using partial pivoting.

This method decomposes the square matrix \( A \) into a lower-triangular matrix \( L \) and an upper-triangular matrix \( U \), such that \( P A = L U \), where \( P \) is a permutation matrix that encodes any row interchanges. The resulting \( L \) and \( U \) matrices are stored in-place in the input matrix.

Parameters

[in,out]	A	the matrix to decompose. Matrix is decomposed in-place.
[in,out]	pivots	buffer to store pivoting, e.g., row-interchanges.

Returns: rc return code, LU_SUCCESS if operation is successful

Note

The method returns the following error codes if an error occurs:

LU_NONSQUARE_MATRIX when input matrix is not square
LU_SINGULAR_MATRIX when the input matrix is singular

The matrix is decomposed in-place.

Precondition: pivots != nullptr; pivots must be able to hold A.getNumRows() entries; A.isSquare()==true

References axom::utilities::abs(), axom::numerics::Matrix< T >::getNumRows(), axom::utilities::isNearlyEqual(), axom::numerics::Matrix< T >::isSquare(), LU_NONSQUARE_MATRIX, LU_SINGULAR_MATRIX, LU_SUCCESS, and axom::numerics::Matrix< T >::swapRows().

◆ lu_solve()

template<typename T >

int axom::numerics::lu_solve	(	const Matrix< T > &	A,
		const int *	pivots,
		const T *	b,
		T *	x
	)

Solve the system \( Ax=b \) by back-substitution, where, A is an LU decomposed matrix produced via a call to lu_decompose().

Parameters

[in]	A	an LU decomposed square matrix ( output from lu_decompose() )
[in]	pivots	stores row-interchanges ( output from lu_decompose() )
[in]	b	the vector on the right-hand side.
[out]	x	the solution vector ( computed ).

Returns: rc return code, LU_SUCCESS if operation is successful

Precondition: A.isSquare()==true; pivots != nullptr; b != nullptr; x != nullptr

References axom::deallocate(), axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), LU_NONSQUARE_MATRIX, LU_SUCCESS, and axom::utilities::swap().

◆ operator<<()

template<typename T >

std::ostream & axom::numerics::operator<<	(	std::ostream &	os,
		const Matrix< T > &	A
	)

Overloaded output stream operator. Outputs the matrix coefficients in to the given output stream.

Parameters

[in,out]	os	output stream object.
[in]	A	user-supplied matrix instance.

Returns: os the updated output stream object.

References axom::numerics::Matrix< T >::getNumColumns(), and axom::numerics::Matrix< T >::getNumRows().

◆ operator==()

template<typename T >

bool axom::numerics::operator==	(	const Matrix< T > &	lhs,
		const Matrix< T > &	rhs
	)

Checks if two matrices are component-wise equal.

Parameters

[in]	lhs	first matrix
[in]	rhs	second matrix

Returns: status true if lhs==rhs, otherwise, false.

References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::utilities::isNearlyEqual().

◆ operator!=()

template<typename T >

bool axom::numerics::operator!=	(	const Matrix< T > &	lhs,
		const Matrix< T > &	rhs
	)

Checks if two matrices are not component-wise equal.

Parameters

[in]	lhs	first matrix
[in]	rhs	second matrix

Returns: status true if lhs!=rhs, otherwise, false.

◆ lower_triangular()

template<typename T >

Matrix< T > axom::numerics::lower_triangular	(	const Matrix< T > &	A,
		bool	unit_diagonal = `false`
	)

Extracts the lower triangular part of a square matrix.

Parameters

[in]	A	a square matrix
[in]	unit_diagonal	indicates if the diagonal entries are implicitly set to ones (1s) or copied from the original matrix, A. Default is false.

Returns: L a new matrix instance of the same dimensions as A consisting of the lower triangular part of A and all the rest of its entries set to zero.

Precondition: A.isSquare() == true

Postcondition: upper-triangular of matrix L is set to zeros (0s)

References axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), and axom::numerics::Matrix< T >::zeros().

◆ upper_triangular()

template<typename T >

Matrix< T > axom::numerics::upper_triangular	(	const Matrix< T > &	A,
		bool	unit_diagonal = `true`
	)

Extract the upper triangular part of a square matrix.

Parameters

[in]	A	a square matrix
[in]	unit_diagonal	indicates if the diagonal entries are implicitly set to ones (1s) or copied from the original matrix, A. Default is true.

Returns: U a new matrix instance of the same dimensions as A consisting of the upper triangulart part of A and all the rest of its entries set to zero.

Precondition: A.isSquare() == true

Postcondition: lower-triangular of matrix U is set to zeros (0s)

References axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), and axom::numerics::Matrix< T >::zeros().

◆ linspace()

template<typename T >

bool axom::numerics::linspace	(	const T &	x0,
		const T &	x1,
		T *	v,
		int	N
	)

inline

Generates a vector consisting of a sequence of N uniformly spaced values over the interval [x0,x1].

Parameters

[in]	x0	scalar, the start value of the sequence.
[in]	x1	scalar, the end value of the sequence.
[out]	v	user-supplied buffer where to store the sequence of numbers
[in]	N	the size of the computed vector sequence.

Returns: status true if successful, otherwise, false.

Note: The output vector, v, must be able to hold at least N values.; if x0 < x1, the sequence of values will be in ascending order.; if x0 > x1, the sequence of values will be in descending order.

Precondition: v != nullptr; N > 1

References AXOM_STATIC_ASSERT_MSG.

◆ cross_product()

template<typename T >

AXOM_HOST_DEVICE void axom::numerics::cross_product	(	const T *	u,
		const T *	v,
		T *	w
	)

inline

Computes the vector cross-product of two vectors, u and v.

Parameters

[in]	u	array pointer to the vector u
[in]	v	array pointer to the vector v
[out]	w	array pointer where to store the cross-product

Note: The u, v, and w are expected to be 3D vectors and are expected to be pointing to array of at least size 3.

Precondition: u != nullptr; v != nullptr; w != nullptr

References determinant().

◆ dot_product()

template<typename T >

AXOM_HOST_DEVICE T axom::numerics::dot_product	(	const T *	u,
		const T *	v,
		int	dim
	)

inline

Computes the dot product of the arrays u and v.

Template Parameters

T data type

Parameters

[in]	u	pointer to an array of size dim
[in]	v	pointer to an array of size dim
[in]	dim	the dimension of the arrays at u and v

Returns: the dot product of the arrays u and v

Precondition: dim >= 1; u != nullptr; v != nullptr; u has at least dim entries; v has at least dim entries

◆ make_orthogonal()

template<typename T >

void axom::numerics::make_orthogonal	(	T *	u,
		T *	v,
		int	dim,
		double	tol = `1E-16`
	)

Makes u orthogonal to v.

Template Parameters

T data type

Parameters

[in,out]	u	vector to be made orthogonal to other; saves in-place
[in]	v	vector that u is made orthogonal to
[in]	dim	dimension of vectors
[in]	tol	tolerance; if the norm of v is less than tol we do nothing

Precondition: dim >= 1; u != nullptr; v != nullptr; T is a floating point type

References AXOM_STATIC_ASSERT_MSG, and dot_product().

◆ orthonormalize()

template<typename T >

bool axom::numerics::orthonormalize	(	T *	basis,
		int	size,
		int	dim,
		double	eps = `1E-16`
	)

Performs Gram-Schmidt orthonormalization in-place on a 2D array of shape size,dim where it treats the rows as the individual vectors.

Template Parameters

T data type

Parameters

[in,out]	basis	vectors to be made orthonormal; saves them in-place
[in]	size	number of vectors
[in]	dim	dimension of vectors
[in]	eps	If one of the vectors, after being made orthogonal to the others, has norm less than eps, then the orthonormalization is declared a failure. Note that this may well destory the data pointed to by basis.

Returns: true if the orthonormalization is successful, false otherwise

Precondition: dim >= 1; 1 <= size <= dim; basis != nullptr; T is a floating point type

References AXOM_STATIC_ASSERT_MSG, make_orthogonal(), and normalize().

◆ normalize()

template<typename T >

AXOM_HOST_DEVICE bool axom::numerics::normalize	(	T *	v,
		int	dim,
		double	eps = `1e-16`
	)

inline

Normalizes the passed in array.

Template Parameters

T data type

Parameters

[in]	v	pointer the array
[in]	dim	the dimension of v
[in]	eps	fuzz parameter

Note: If squared norm of v is less than eps, the normalization fails and we do nothing to v

Returns: true if normalization is successful, false otherwise.

Precondition: dim >= 1; v != nullptr; T is a floating point type

References AXOM_STATIC_ASSERT_MSG, and dot_product().

◆ matrix_add()

template<typename T >

bool axom::numerics::matrix_add	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		Matrix< T > &	C
	)

inline

Computes \( \mathcal{C} = \mathcal{A} + \mathcal{B} \).

Parameters

[in]	A	\( M \times N \) matrix on the left-hand side.
[in]	B	\( M \times N \) matrix on the right-hand side.
[out]	C	\( M \times N \) output matrix.

Postcondition: \( c_{ij} = a_{ij} + b{ij} \), \(\forall c_{ij} \in \mathcal{C}\)

Note: Matrix addition is undefined for matrices that have different dimensions. If the input matrices, \( \mathcal{A} \), \( \mathcal{B} \) or \( \mathcal{C} \) have different dimensions, a \( 1 \times 1 \) null matrix is returned in \( \mathcal{C} \)

Returns: status true if addition is successful, otherwise, false.

References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().

◆ matrix_subtract()

template<typename T >

bool axom::numerics::matrix_subtract	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		Matrix< T > &	C
	)

inline

Computes \( \mathcal{C} = \mathcal{A} - \mathcal{B} \).

Parameters

[in]	A	\( M \times N \) matrix on the left-hand side.
[in]	B	\( M \times N \) matrix on the right-hand side.
[out]	C	\( M \times N \) output matrix.

Postcondition: \( c_{ij} = a_{ij} + b{ij} \), \(\forall c_{ij} \in \mathcal{C}\)

Note: Matrix subtraction is undefined for matrices that have different dimensions. If the input matrices, \( \mathcal{A} \), \( \mathcal{B} \) or \( \mathcal{C} \) have different dimensions, a \( 1 \times 1 \) null matrix is returned in \( \mathcal{C} \)

Returns: status true if addition is successful, otherwise, false.

References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().

◆ matrix_multiply()

template<typename T >

bool axom::numerics::matrix_multiply	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		Matrix< T > &	C
	)

inline

Computes the matrix-matrix product of \( \mathcal{A} \) and \( \mathcal{B} \) and stores the result in \( \mathcal{C} \).

Parameters

[in]	A	\( M \times N \) matrix on the left hand-side.
[in]	B	\( N \times K \) matrix on the right hand-side.
[out]	C	\( M \times K \) output matrix

Precondition: The inner dimensions of matrices A, B must match; Output matrix should be an \( M \times K \) matrix

Note: Matrix multiplication is undefined for matrices with different inner dimension. If the inner dimensions are not matching, the code returns a \( 1 \times 1 \) null matrix in \( \mathcal{C} \)

Note
Matrix \( \mathcal{C} \) should be a different Matrix instance than \( \mathcal{A} \) , \( \mathcal{B} \)

Returns
status true if the multiplication is successful, otherwise, false.

References axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().

◆ matrix_scalar_multiply()

template<typename T >

void axom::numerics::matrix_scalar_multiply	(	Matrix< T > &	A,
		const T &	c
	)

inline

Computes a scalar-matrix produect of a matrix \( \mathcal{A} \) and a scalar \( c \) and stores the result in \( \mathcal{A} \).

Parameters

[in,out]	A	an \( M \times N \) matrix
[in]	c	scalar value to multiply the matrix coefficients

Postcondition: \( a_{ij} = c \cdot a_{ij} \), \( \forall a_{ij} \in mathcal{A} \)

References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), and axom::numerics::Matrix< T >::getNumRows().

◆ matrix_vector_multiply()

template<typename T >

void axom::numerics::matrix_vector_multiply	(	const Matrix< T > &	A,
		const T *	vec,
		T *	output
	)

inline

Computes the matrix-vector product of a matrix \(\mathcal{A}\) and a vector \(\mathbf{x}\) and store the result in the user-supplied output vector.

Parameters

[in]	A	an \( M \times N \) matrix
[in]	vec	pointer to user-supplied vector storing the vector
[out]	output	pointer to the user-supplied output vector

Precondition: vec != nullptr; vec must be of dimension \( N \); output != nullptr; output must be of dimension \( M \); vec != output

Postcondition: \( b_i = \sum\limits_{j=0}^N a_{ij} \cdot x_j \), \(\forall i \in [0,M-1] \)

References axom::numerics::Matrix< T >::getNumColumns(), and axom::numerics::Matrix< T >::getNumRows().

◆ matrix_transpose()

template<typename T >

bool axom::numerics::matrix_transpose	(	const Matrix< T > &	A,
		Matrix< T > &	M
	)

inline

Computes the matrix transpose of a given matrix \( \mathcal{A} \).

Parameters

[in]	A	an \( m \times n \) matrix
[out]	M	an \( n \times m \) matrix where to store the transpose.

Note: If the supplied matrix does not have the correct dimensions, the code returns a \( 1 \times 1 \) null matrix in \( \mathcal{M} \); Matrix \( \mathcal{M} \) should be a different Matrix instance than \( \mathcal{A} \)

Returns: status true if the matrix transpose is successful, otherwise, false.

References axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().

◆ matrix_norm()

template<typename T >

T axom::numerics::matrix_norm	(	const Matrix< T > &	A,
		MatrixNorm	normType
	)

inline

Computes the specified norm of a given \( M \times N \) matrix.

Parameters

[in]	A	the matrix whose norm is computed
[in]	normType	the type of norm to compute

Returns: norm the computed norm or -1.0 on error.

Note: The computed norm is a non-negative scalar value. This method will return a negative return value to indicate an error, e.g., the supplied matrix has incorrect dimensions.; The second argument specifies the type of norm to compute. Three matrix norm types are supported:

\( P_1\) Norm

The \( P_1\) norm is computed by taking the maximum absolute column sum, given by:
\[ ||A||_1 = \max\limits_{1 \le j \le N}( \sum_{i=1}^M | a_{ij} | ) \]
\(\infty\) Norm

The \(\infty\) norm is computed by taking the maximum absolute row sum, given by:
\[ ||A||_\infty = \max\limits_{1 \le j \le M}( \sum_{i=1}^N | a_{ij} | ) \]
Frobenius Norm

The frobenius norm of an \( M \times N \) matrix, is given by:
\[ ||A||_F = \sqrt{ \sum_{i=1}^M \sum_{j=1}^N ({a_{ij}})^2 } \]

Precondition: A.getNumRows() >= 2; A.getNumCols() >= 2

See also: MatrixNorm

References FROBENIUS_NORM, axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), INF_NORM, and P1_NORM.

◆ evaluate_polynomial()

template<typename ScalarType >

ScalarType axom::numerics::evaluate_polynomial	(	ArrayView< const double >	coeffs_descending,
		const ScalarType &	x
	)

Evaluate a polynomial with coefficients stored in descending power order.

Template Parameters

ScalarType The scalar type used for evaluation, e.g. double or std::complex<double>.

Parameters

[in]	coeffs_descending	Polynomial coefficients in descending power order.
[in]	x	The evaluation point.

Returns: The polynomial value at x.

◆ bernstein_to_monomial()

axom::Array<double> axom::numerics::bernstein_to_monomial ( ArrayView< const double > bernstein_coeffs )

Convert Bernstein coefficients to monomial coefficients on [0, 1].

Given Bernstein coefficients \( b_i \) of degree \( n \), this returns monomial coefficients \( a_i \) such that \( \sum_i b_i B_i^n(t) = \sum_i a_i t^i \).

Parameters

[in] bernstein_coeffs Coefficients in the Bernstein basis.

Returns: Coefficients in ascending monomial order.

◆ effective_polynomial_degree()

int axom::numerics::effective_polynomial_degree	(	ArrayView< const double >	coeffs_ascending,
		double	tol = `1e-12`
	)

Return the effective degree of a polynomial after trimming near-zero leading coefficients.

Parameters

[in]	coeffs_ascending	Polynomial coefficients in ascending power order.
[in]	tol	Relative trimming tolerance. Coefficients with magnitude at or below `tol * max(abs(coeffs_ascending))` are treated as zero when trimming the highest-degree terms.

Returns: The largest exponent whose coefficient exceeds the scaled trimming threshold, or 0 if the polynomial is constant to within the requested tolerance.

◆ solve_polynomial_durand_kerner()

axom::Array<std::complex<double> > axom::numerics::solve_polynomial_durand_kerner	(	ArrayView< const double >	coeffs_ascending,
		double	tol = `1e-12`,
		std::complex< double >	seed = `std::complex< double > {0.4, 0.9}`
	)

Approximate all roots of a polynomial using the Durand-Kerner method.

Parameters

[in]	coeffs_ascending	Polynomial coefficients in ascending power order.
[in]	tol	Iteration tolerance and relative effective-zero threshold used when trimming leading coefficients.
[in]	seed	Seed used to generate deterministic, distinct initial guesses. This should be a complex number on or near the unit circle (and not purely real), so that `seed^(i+1)` spreads the starting points around a reasonable radius.

Returns: The polynomial roots, sorted by real part and then imaginary part, or an empty array if the iteration does not converge.

◆ solve_polynomial_durand_kerner_checked()

PolynomialRootResult axom::numerics::solve_polynomial_durand_kerner_checked	(	ArrayView< const double >	coeffs_ascending,
		double	tol = `1e-12`,
		int	max_iters = `200`,
		std::complex< double >	seed = `std::complex< double > {0.4, 0.9}`
	)

Approximate all roots of a polynomial using the Durand-Kerner method, including convergence diagnostics.

Parameters

[in]	coeffs_ascending	Polynomial coefficients in ascending power order.
[in]	tol	Iteration tolerance and relative effective-zero threshold used when trimming leading coefficients.
[in]	max_iters	Maximum number of Durand-Kerner iterations. Default is 200, which provides reliable convergence for well-conditioned polynomials up to degree ~15-20. Higher-degree or ill-conditioned polynomials (repeated roots, clustered roots) may require more iterations or specialized methods.
[in]	seed	Seed used to generate deterministic, distinct initial guesses. This should be a complex number on or near the unit circle (and not purely real), so that `seed^(i+1)` spreads the starting points around a reasonable radius.

Returns: The roots and convergence metadata for the solve.

◆ solve_linear() [1/2]

int axom::numerics::solve_linear	(	ArrayView< const double >	coeff,
		ArrayView< double >	roots,
		int &	numRoots
	)

Find the real root for a linear equation of form \( ax + b = 0 \).

Parameters

[in]	coeff	Equation coefficients: coeff[i] multiplies \( x^i \). So, coeff[0] = b and coeff[1] = a.
[out]	roots	The real roots of the equation.
[out]	numRoots	The number of distinct, real roots found (max: 1). If the line lies on the X-axis, numRoots is assigned -1 to indicate infinitely many solutions. If the line does not intersect the X-axis, numRoots is assigned 0.

Returns: 0 for success, or -1 to indicate an inconsistent equation (all coefficients 0 except for coeff[0]).

Precondition: coeff.size() >= 2; roots.size() >= 1

◆ solve_linear() [2/2]

int axom::numerics::solve_linear	(	const double *	coeff,
		double *	roots,
		int &	numRoots
	)

Deprecated pointer-based overload for solve_linear(ArrayView<const double>, ArrayView<double>, int&).

◆ solve_quadratic() [1/2]

int axom::numerics::solve_quadratic	(	ArrayView< const double >	coeff,
		ArrayView< double >	roots,
		int &	numRoots
	)

For a quadratic equation of the form \( ax^2 + bx + c = 0 \), find real roots using the quadratic formula.

Parameters

[in]	coeff	Equation coefficients: coeff[i] multiplies \( x^i \). So, coeff[0] = c, coeff[1] = b, coeff[2] = a.
[out]	roots	The real roots of the equation.
[out]	numRoots	The number of distinct, real roots found (max: 2). If the equation degenerates to a line lying on the X-axis, numRoots is assigned -1 to indicate infinitely many solutions. If the equation does not intersect the X-axis, numRoots is assigned 0.

Returns: 0 for success, or -1 to indicate an inconsistent equation (all coefficients 0 except for coeff[0]).

Precondition: coeff.size() >= 3; roots.size() >= 2

◆ solve_quadratic() [2/2]

int axom::numerics::solve_quadratic	(	const double *	coeff,
		double *	roots,
		int &	numRoots
	)

Deprecated pointer-based overload for solve_quadratic(ArrayView<const double>, ArrayView<double>, int&).

◆ solve_cubic() [1/2]

int axom::numerics::solve_cubic	(	ArrayView< const double >	coeff,
		ArrayView< double >	roots,
		int &	numRoots
	)

For a cubic equation of the form \( ax^3 + bx^2 + cx + d = 0 \), find real roots.

A closed-form solution for cubic equations was published in Cardano's Ars Magna of 1545. This can be summarized as follows:

Start with the cubic equation

\[ ax^3 + bx^2 + cx + d = 0. \]

Divide through by a to normalize, then define the following:

\[ p = b^2 - 3c \]

\[ q = -\frac{27}{2}d - b^3 + \frac{9}{2}cb \]

\[ t = 2p^{-3/2}q \]

\[ y = \frac{3}{\sqrt{p}}(x + \frac{1}{3}b) \]

Then the original cubic equation can be written as \( y^3 - 3y = t \) with a solution \( y = \frac{1}{u} + u \), with

\[ u = \sqrt[3]{\frac{t}{2} \pm \sqrt{\frac{t^2}{4} - 1}}. \]

Because the cubic is an odd function, there can be either one, two, or three distinct real roots. If there is one distinct real root, the root can either have multiplicity 3, or have multiplicity 1 with two complex roots. In the case of two real roots, one root will have multiplicity 2. If the discriminant \( d = -27t^2 - 4(-3^3) \) is zero, the equation has at least one root with multiplicity greater than 1. If \( d > 0, \) there are three distinct real roots, and if \( d < 0, \) there is one real root.

See J. Kopp, Efficient numerical diagonalization of hermitian 3x3 matrices, Int.J.Mod.Phys. C19:523-548, 2008 (https://arxiv.org/abs/physics/0610206) and G. A. Korn and T. M. Korn, "Mathematical Handbook for Scientists and Engineers," QA37 K84 1968 in the library; section 1.8-3 "Cubic Equations", p. 23.

Parameters

[in]	coeff	Equation coefficients: coeff[i] multiplies \( x^i \). So, coeff[0] = d, coeff[1] = c, coeff[2] = b, coeff[3] = a.
[out]	roots	The real roots of the equation.
[out]	numRoots	The number of distinct, real roots found (max: 3). If the equation degenerates to a line lying on the X-axis, numRoots is assigned -1 to indicate infinitely many solutions. If the equation does not intersect the X-axis, numRoots is assigned 0.

Returns: 0 for success, or -1 to indicate an inconsistent equation (all coefficients 0 except for coeff[0]).

Precondition: coeff.size() >= 4; roots.size() >= 3

◆ solve_cubic() [2/2]

int axom::numerics::solve_cubic	(	const double *	coeff,
		double *	roots,
		int &	numRoots
	)

Deprecated pointer-based overload for solve_cubic(ArrayView<const double>, ArrayView<double>, int&).

◆ compute_gauss_legendre_data()

void axom::numerics::compute_gauss_legendre_data	(	int	npts,
		axom::Array< double > &	nodes,
		axom::Array< double > &	weights,
		int	allocatorID = `axom::getDefaultAllocatorID()`
	)

Computes a 1D quadrature rule of Gauss-Legendre points.

Parameters

[in]	npts	The number of points in the rule
[out]	nodes	The array of 1D nodes
[out]	weights	The array of weights

A Gauss-Legendre rule with npts points can exactly integrate polynomials of order 2 * npts - 1

Algorithm adapted from the MFEM implementation in mfem/fem/intrules.cpp

Note: This method constructs the points by scratch each time, without caching

See also: get_gauss_legendre(int)

◆ get_gauss_legendre()

QuadratureRule axom::numerics::get_gauss_legendre	(	int	npts,
		int	allocatorID = `axom::getDefaultAllocatorID()`
	)

Computes or accesses a precomputed 1D quadrature rule of Gauss-Legendre points.

Parameters

[in] npts The number of points in the rule

A Gauss-Legendre rule with npts points can exactly integrate polynomials of order 2 * npts - 1

Note: If this method has already been called for a given order, it will reuse the same quadrature points without needing to recompute them

Returns: The QuadratureRule object which contains axom::ArrayView<double>'s of stored nodes and weights

Variable Documentation

◆ JACOBI_DEFAULT_TOLERANCE

constexpr double axom::numerics::JACOBI_DEFAULT_TOLERANCE = 1.e-18

constexpr

◆ JACOBI_DEFAULT_MAX_ITERATIONS

constexpr int axom::numerics::JACOBI_DEFAULT_MAX_ITERATIONS = 20

constexpr

◆ JACOBI_EIGENSOLVE_SUCCESS

constexpr int axom::numerics::JACOBI_EIGENSOLVE_SUCCESS = 0

constexpr

◆ JACOBI_EIGENSOLVE_FAILURE

constexpr int axom::numerics::JACOBI_EIGENSOLVE_FAILURE = -1

constexpr

Namespaces

Classes

Enumerations

Functions

Variables

Enumeration Type Documentation

◆ ReturnCodes

◆ MatrixNorm

Function Documentation

◆ determinant() [1/4]

◆ determinant() [2/4]

◆ determinant() [3/4]

◆ determinant() [4/4]

◆ eigen_solve()

◆ eigen_sort()

◆ jacobi_eigensolve()

◆ linear_solve()

◆ lu_decompose()

◆ lu_solve()

◆ operator<<()

◆ operator==()

◆ operator!=()

◆ lower_triangular()

◆ upper_triangular()

◆ linspace()

◆ cross_product()

◆ dot_product()

◆ make_orthogonal()

◆ orthonormalize()

◆ normalize()

◆ matrix_add()

◆ matrix_subtract()

◆ matrix_multiply()

◆ matrix_scalar_multiply()

◆ matrix_vector_multiply()

◆ matrix_transpose()

◆ matrix_norm()

◆ evaluate_polynomial()

◆ bernstein_to_monomial()

◆ effective_polynomial_degree()

◆ solve_polynomial_durand_kerner()

◆ solve_polynomial_durand_kerner_checked()

◆ solve_linear() [1/2]

◆ solve_linear() [2/2]

◆ solve_quadratic() [1/2]

◆ solve_quadratic() [2/2]

◆ solve_cubic() [1/2]

◆ solve_cubic() [2/2]

◆ compute_gauss_legendre_data()

◆ get_gauss_legendre()

Variable Documentation

◆ JACOBI_DEFAULT_TOLERANCE

◆ JACOBI_DEFAULT_MAX_ITERATIONS

◆ JACOBI_EIGENSOLVE_SUCCESS

◆ JACOBI_EIGENSOLVE_FAILURE