AXOM
Axom provides a robust, flexible software infrastructure for the development of multi-physics applications and computational tools.
|
Classes | |
struct | floating_point_limits |
Traits class for floating point types to query limit information. More... | |
struct | floating_point_limits< double > |
struct | floating_point_limits< float > |
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... | |
template<typename T > | |
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 > | |
T | dot_product (const T *u, const T *v, int dim) |
Computes the dot product of the arrays u and v. 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 > | |
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 > | |
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... | |
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 > | |
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 |
|
inline |
Computes a 2X2 determinant of the given matrix.
[in] | a00 | matrix element |
[in] | a01 | matrix element |
[in] | a10 | matrix element |
[in] | a11 | matrix element |
Referenced by axom::primal::Vector< SpaceCoordType, NDIMS >::cross_product(), cross_product(), determinant(), axom::primal::in_sphere(), linear_solve(), axom::primal::orientation(), axom::primal::Tetrahedron< T, NDIMS >::physToBarycentric(), and axom::primal::Triangle< T, NDIMS >::physToBarycentric().
|
inline |
Computes the 3x3 determinant of the given matrix.
[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 |
|
inline |
Computes the 4x4 determinant of the given matrix.
[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 |
real axom::numerics::determinant | ( | const Matrix< real > & | A | ) |
Computes the determinant of the given square matrix.
[in] | A | an \( N \times N \) input matrix |
References A, determinant(), axom::numerics::Matrix< T >::empty(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::isSquare(), lu_decompose(), and LU_SUCCESS.
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:
Because step 1 and step 3 require random numbers, this function requires the caller to properly seed the randomizer (for example, by calling srand()).
[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 |
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().
bool axom::numerics::eigen_sort | ( | T * | lamdas, |
Matrix< T > & | eigen_vectors | ||
) |
Sorts the supplied eigenvalues and eigenvectors in ascending order.
[in,out] | lambdas | pointer to buffer consisting of the eigenvalues |
[in,out] | eigen_vectors | matrix of the corresponding eigenvectors |
References axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), axom::utilities::swap(), and axom::numerics::Matrix< T >::swapColumns().
Referenced by jacobi_eigensolve().
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.
[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) |
References A, axom::utilities::abs(), AXOM_STATIC_ASSERT_MSG, axom::deallocate(), eigen_sort(), axom::numerics::Matrix< T >::getNumRows(), axom::utilities::isNearlyEqual(), axom::numerics::Matrix< T >::isSquare(), and JACOBI_EIGENSOLVE_FAILURE.
int axom::numerics::linear_solve | ( | Matrix< T > & | A, |
const T * | b, | ||
T * | x | ||
) |
Solves a linear system of the form \( Ax=b \).
[in] | A | a square input matrix |
[in] | b | the right-hand side |
[out] | x | the solution vector (computed) |
References A, determinant(), axom::numerics::Matrix< T >::getNumColumns(), axom::utilities::isNearlyEqual(), axom::numerics::Matrix< T >::isSquare(), lu_decompose(), LU_NONSQUARE_MATRIX, lu_solve(), and LU_SUCCESS.
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.
[in,out] | A | the matrix to decompose. Matrix is decomposed in-place. |
[in,out] | pivots | buffer to store pivoting, e.g., row-interchanges. |
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().
Referenced by determinant(), and linear_solve().
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().
[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 ). |
References A, axom::deallocate(), axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), LU_NONSQUARE_MATRIX, LU_SUCCESS, and axom::utilities::swap().
Referenced by linear_solve().
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.
[in,out] | os | output stream object. |
[in] | A | user-supplied matrix instance. |
References axom::numerics::Matrix< T >::getNumRows().
Matrix< T > axom::numerics::lower_triangular | ( | const Matrix< T > & | A, |
bool | unit_diagonal = false |
||
) |
Extracts the lower triangular part of a square matrix.
[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. |
References A, axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), and axom::numerics::Matrix< T >::zeros().
Matrix< T > axom::numerics::upper_triangular | ( | const Matrix< T > & | A, |
bool | unit_diagonal = true |
||
) |
Extract the upper triangular part of a square matrix.
[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. |
References A, axom::numerics::Matrix< T >::getNumRows(), axom::numerics::Matrix< T >::isSquare(), and axom::numerics::Matrix< T >::zeros().
|
inline |
Generates a vector consisting of a sequence of N uniformly spaced values over the interval [x0,x1].
[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. |
References AXOM_STATIC_ASSERT_MSG.
|
inline |
Computes the vector cross-product of two vectors, u and v.
[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 |
References determinant().
Referenced by axom::primal::Vector< SpaceCoordType, NDIMS >::cross_product(), and axom::primal::Plane< T, NDIMS >::Plane().
|
inline |
Computes the dot product of the arrays u and v.
T | data type |
[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 |
Referenced by axom::primal::Plane< T, NDIMS >::computeSignedDistance(), axom::primal::Vector< SpaceCoordType, NDIMS >::dot_product(), eigen_solve(), axom::primal::Sphere< T, NDIMS >::intersectsWith(), make_orthogonal(), normalize(), and axom::primal::Plane< T, NDIMS >::Plane().
void axom::numerics::make_orthogonal | ( | T * | u, |
T * | v, | ||
int | dim, | ||
double | tol = 1E-16 |
||
) |
Makes u orthogonal to v.
T | data type |
[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 |
References AXOM_STATIC_ASSERT_MSG, and dot_product().
Referenced by orthonormalize().
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.
T | data type |
[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. |
References AXOM_STATIC_ASSERT_MSG, make_orthogonal(), and normalize().
|
inline |
Normalizes the passed in array.
T | data type |
[in] | v | pointer the array |
[in] | dim | the dimension of v |
[in] | eps | fuzz parameter |
References AXOM_STATIC_ASSERT_MSG, and dot_product().
Referenced by orthonormalize(), and axom::primal::Plane< T, NDIMS >::Plane().
|
inline |
Computes \( \mathcal{C} = \mathcal{A} + \mathcal{B} \).
[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. |
References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().
|
inline |
Computes \( \mathcal{C} = \mathcal{A} - \mathcal{B} \).
[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. |
References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().
|
inline |
Computes the matrix-matrix product of \( \mathcal{A} \) and \( \mathcal{B} \) and stores the result in \( \mathcal{C} \).
[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 |
References A, B, C, axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().
|
inline |
Computes a scalar-matrix produect of a matrix \( \mathcal{A} \) and a scalar \( c \) and stores the result in \( \mathcal{A} \).
[in,out] | A | an \( M \times N \) matrix |
[in] | c | scalar value to multiply the matrix coefficients |
References axom::numerics::Matrix< T >::data(), axom::numerics::Matrix< T >::getNumColumns(), and axom::numerics::Matrix< T >::getNumRows().
Referenced by axom::primal::OrientedBoundingBox< T, NDIMS >::OrientedBoundingBox().
|
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.
[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 |
References A, axom::numerics::Matrix< T >::getNumColumns(), and axom::numerics::Matrix< T >::getNumRows().
Referenced by eigen_solve().
|
inline |
Computes the matrix transpose of a given matrix \( \mathcal{A} \).
[in] | A | an \( m \times n \) matrix |
[out] | M | an \( n \times m \) matrix where to store the transpose. |
References A, axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), and axom::numerics::Matrix< T >::zeros().
|
inline |
Computes the specified norm of a given \( M \times N \) matrix.
[in] | A | the matrix whose norm is computed |
[in] | normType | the type of norm to compute |
\( 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 } \]
References FROBENIUS_NORM, axom::numerics::Matrix< T >::getNumColumns(), axom::numerics::Matrix< T >::getNumRows(), INF_NORM, and P1_NORM.
|
inline |
Computes the vector cross-product of two vectors, u and v.
[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 |
References determinant().
Referenced by axom::primal::Vector< SpaceCoordType, NDIMS >::cross_product(), and axom::primal::Plane< T, NDIMS >::Plane().
|
inline |
Computes the dot product of the arrays u and v.
T | data type |
[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 |
Referenced by axom::primal::Plane< T, NDIMS >::computeSignedDistance(), axom::primal::Vector< SpaceCoordType, NDIMS >::dot_product(), eigen_solve(), axom::primal::Sphere< T, NDIMS >::intersectsWith(), make_orthogonal(), normalize(), and axom::primal::Plane< T, NDIMS >::Plane().
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 \).
[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. |
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.
[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. |
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.
[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. |
constexpr double axom::numerics::JACOBI_DEFAULT_TOLERANCE = 1.e-18 |
constexpr int axom::numerics::JACOBI_DEFAULT_MAX_ITERATIONS = 20 |
constexpr int axom::numerics::JACOBI_EIGENSOLVE_SUCCESS = 0 |
constexpr int axom::numerics::JACOBI_EIGENSOLVE_FAILURE = -1 |
Referenced by jacobi_eigensolve().