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

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

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

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

## ◆ 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

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

## ◆ 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

## ◆ 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

## ◆ 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] 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 eigenvales 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())

## ◆ 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

## ◆ 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

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

## ◆ 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

## ◆ 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

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

## ◆ 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)

## ◆ 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)

## ◆ 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().

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.

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

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

## ◆ 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}$$

## ◆ 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]$$

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

## ◆ 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
MatrixNorm

## ◆ 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] 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 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.

 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] 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 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:

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

## ◆ 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