AXOM
Axom provides a robust, flexible software infrastructure for the development of multi-physics applications and computational tools.
axom::numerics Namespace Reference

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

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...
 
int solve_linear (const double *coeff, double *roots, int &numRoots)
 Find the real root for a linear equation of form \( ax + b = 0 \). More...
 
int solve_quadratic (const double *coeff, 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_cubic (const double *coeff, double *roots, int &numRoots)
 For a cubic equation of the form \( ax^3 + bx^2 + cx + d = 0 \), find real roots. 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 >
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_DEVICEdot_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

Enumerator
LU_SUCCESS 
LU_SINGULAR_MATRIX 
LU_NONSQUARE_MATRIX 

◆ 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]a00matrix element
[in]a01matrix element
[in]a10matrix element
[in]a11matrix element
Returns
det the determinant of the 2X2 matrix

◆ 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]a00matrix element
[in]a01matrix element
[in]a02matrix element
[in]a10matrix element
[in]a11matrix element
[in]a12matrix element
[in]a20matrix element
[in]a21matrix element
[in]a22matrix element
Returns
det the determinant of the 3x3 matrix.

◆ 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]a00matrix element
[in]a01matrix element
[in]a02matrix element
[in]a03matrix element
[in]a10matrix element
[in]a11matrix element
[in]a12matrix element
[in]a13matrix element
[in]a20matrix element
[in]a21matrix element
[in]a22matrix element
[in]a23matrix element
[in]a30matrix element
[in]a31matrix element
[in]a32matrix element
[in]a33matrix element
Returns
det the determinant of the 4x4 matrix

◆ determinant() [4/4]

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

Computes the determinant of the given square matrix.

Parameters
[in]Aan \( 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:

  1. Generate a random vector
  2. Make orthonormal to previous eigenvectors if any
  3. Run depth iterations of power method
  4. 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]Aa square input matrix
[in]knumber of eigenvalue-eigenvectors to find
[out]upointer to k eigenvectors in order by magnitude of eigenvalue
[out]lambdaspointer to k eigenvales in order by size
[in]numIterationsoptional 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]lambdaspointer to buffer consisting of the eigenvalues
[in,out]eigen_vectorsmatrix 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,
TOL = JACOBI_DEFAULT_TOLERANCE 
)

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

Parameters
[in]Athe input matrix whose eigenpairs will be computed
[out]Vmatrix to store the eigenvectors of the input matrix
[out]lambdasbuffer where the computed eigenvalues will be stored.
[in]maxIterationsthe maximum number of iterations (optional)
[out]numIterationsthe number of actual iterations (optional)
[in]TOLconvergence 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]Aa square input matrix
[in]bthe right-hand side
[out]xthe 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]Athe matrix to decompose. Matrix is decomposed in-place.
[in,out]pivotsbuffer 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]Aan LU decomposed square matrix ( output from lu_decompose() )
[in]pivotsstores row-interchanges ( output from lu_decompose() )
[in]bthe vector on the right-hand side.
[out]xthe 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]osoutput stream object.
[in]Auser-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]lhsfirst matrix
[in]rhssecond 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]lhsfirst matrix
[in]rhssecond 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]Aa square matrix
[in]unit_diagonalindicates 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]Aa square matrix
[in]unit_diagonalindicates 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]x0scalar, the start value of the sequence.
[in]x1scalar, the end value of the sequence.
[out]vuser-supplied buffer where to store the sequence of numbers
[in]Nthe 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]uarray pointer to the vector u
[in]varray pointer to the vector v
[out]warray 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
Tdata type
Parameters
[in]upointer to an array of size dim
[in]vpointer to an array of size dim
[in]dimthe 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
Tdata type
Parameters
[in,out]uvector to be made orthogonal to other; saves in-place
[in]vvector that u is made orthogonal to
[in]dimdimension of vectors
[in]toltolerance; 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
Tdata type
Parameters
[in,out]basisvectors to be made orthonormal; saves them in-place
[in]sizenumber of vectors
[in]dimdimension of vectors
[in]epsIf 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
Tdata type
Parameters
[in]vpointer the array
[in]dimthe dimension of v
[in]epsfuzz 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]Aan \( M \times N \) matrix
[in]cscalar 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]Aan \( M \times N \) matrix
[in]vecpointer to user-supplied vector storing the vector
[out]outputpointer 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]Aan \( m \times n \) matrix
[out]Man \( 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]Athe matrix whose norm is computed
[in]normTypethe 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.

◆ solve_linear()

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

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

Parameters
[in]coeffEquation coefficients: coeff[i] multiplies \( x^i \). So, coeff[0] = b and coeff[1] = a.
[out]rootsThe real roots of the equation.
[out]numRootsThe 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 doesn't 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 points to an array of length at least 2.
roots points to an array of length at least 1.

◆ solve_quadratic()

int axom::numerics::solve_quadratic ( const double *  coeff,
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]coeffEquation coefficients: coeff[i] multiplies \( x^i \). So, coeff[0] = c, coeff[1] = b, coeff[2] = a.
[out]rootsThe real roots of the equation.
[out]numRootsThe 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 doesn't 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 points to an array of length at least 3.
roots points to an array of length at least 2.

◆ solve_cubic()

int axom::numerics::solve_cubic ( const double *  coeff,
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]coeffEquation coefficients: coeff[i] multiplies \( x^i \). So, coeff[0] = d, coeff[1] = c, coeff[2] = b, coeff[3] = a.
[out]rootsThe real roots of the equation.
[out]numRootsThe 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 doesn't 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 points to an array of length at least 4.
roots points to an array of length at least 3.

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