Backend-agnostic linear algebra API and LAPACK wrappers. More...

Collaboration diagram for Linear Algebra:

Files
file	src/linalg/fn_append.hpp
	Appends one vector onto another, in place.
file	src/linalg/fn_combine.hpp
	Concatenates two, three or four vectors into one.
file	src/linalg/fn_cross.hpp
	Cross product of two vectors.
file	src/linalg/fn_crossmat.hpp
	Cross product of a normal vector with every column of a matrix.
file	src/linalg/fn_ddpolyval.hpp
	Evaluates the second derivative of a polynomial.
file	src/linalg/fn_det.hpp
	Determinant of a square matrix.
file	src/linalg/fn_dot.hpp
	Inner products of vectors and matrices.
file	src/linalg/fn_dpolyval.hpp
	Evaluates the first derivative of a polynomial.
file	src/linalg/fn_eigen.hpp
	Eigenvalues of a general square matrix.
file	src/linalg/fn_inv.hpp
	Inverse of a square matrix.
file	src/linalg/fn_inv2.hpp
	Closed-form inverse of a 2x2 matrix.
file	src/linalg/fn_inv3.hpp
	Closed-form inverse of a 3x3 matrix.
file	src/linalg/fn_linspace.hpp
	Evenly spaced values over an interval.
file	src/linalg/fn_max.hpp
	Largest entry of a vector, matrix or column view.
file	src/linalg/fn_min.hpp
	Smallest entry of a vector, matrix or column view.
file	src/linalg/fn_norm.hpp
	Euclidean length of a vector.
file	src/linalg/fn_polyfit.hpp
	Least-squares polynomial fit.
file	src/linalg/fn_polyval.hpp
	Evaluates a polynomial.
file	src/linalg/fn_r2.hpp
	Coefficient of determination, R^2.
file	src/linalg/fn_reverse.hpp
	Reverses the order of a vector's entries.
file	src/linalg/fn_sort.hpp
	Sorts a vector in place, ascending.
file	src/linalg/fn_sum.hpp
	Sum of the entries of a vector.
file	src/linalg/fn_to_cell.hpp
	Converts a Vector into a Cell.
file	src/linalg/fn_to_vector.hpp
	Converts a Cell into a Vector.
file	src/linalg/fn_trans.hpp
	Matrix transpose.
file	src/linalg/fn_unique.hpp
	Sorts a vector in place and removes duplicate entries.
file	src/linalg/lapack/fn_gees.hpp
	Schur decomposition of a general square matrix (LAPACK ?gees).
file	src/linalg/lapack/fn_geev.hpp
	Eigenvalues and eigenvectors of a general square matrix (LAPACK ?geev).
file	src/linalg/lapack/fn_gels.hpp
	Least-squares or minimum-norm solution of a full-rank system (LAPACK ?gels).
file	src/linalg/lapack/fn_gemm.hpp
	Matrix-matrix product C := alphaop(A)op(B) + beta*C (BLAS ?gemm).
file	src/linalg/lapack/fn_gesv.hpp
	Solves a square linear system by LU factorization (LAPACK ?gesv).
file	src/linalg/lapack/fn_gesvd.hpp
	Singular value decomposition A = U * diag(S) * VT (LAPACK ?gesvd).
file	src/linalg/lapack/fn_getrf.hpp
	LU factorization with partial pivoting, A = PLU (LAPACK ?getrf).
file	src/linalg/lapack/fn_getri.hpp
	Inverts a square matrix in place from its LU factors (LAPACK ?getri).
file	src/linalg/lapack/fn_posv.hpp
	Solves a real symmetric or complex Hermitian positive-definite system by Cholesky factorization (LAPACK ?posv).
file	src/linalg/lapack/lapacktools.hpp
	Shared helpers for the LAPACK wrappers, such as leading_dimension().

Classes
class	belfem::Matrix< T >
	Dense column-major matrix. More...
class	belfem::Vector< T >
	Column vector. More...

Functions
template<typename T>
void	belfem::append (Vector< T > &aA, Vector< T > &aB)
	Appends one vector onto the end of another, in place.
void	belfem::crossmat (const Vector< real > &aN, const Matrix< real > &aA, Vector< real > &aNxA)
	Crosses a normal vector with every column of a matrix, 2D.
void	belfem::crossmat (const Vector< real > &aN, const Matrix< real > &aA, const real aScale, Vector< real > &aNxA)
	Crosses a normal vector with every column of a matrix, 2D, scaled.
void	belfem::crossmat (const Vector< real > &aN, const Matrix< real > &aA, Matrix< real > &aNxA)
	Crosses a normal vector with every column of a matrix, 3D.
void	belfem::crossmat (const Vector< real > &aN, const Matrix< real > &aA, const real aScale, Matrix< real > &aNxA)
	Crosses a normal vector with every column of a matrix, 3D, scaled.
template<typename T>
T	belfem::ddpolyval (const Vector< T > &aCoeffs, const T aX)
	Evaluates the second derivative of a polynomial at one point.
template<typename T>
T	belfem::dpolyval (const Vector< T > &aCoeffs, const T aX)
	Evaluates the first derivative of a polynomial at one point.
real	belfem::inv2 (const Matrix< real > &aA, Matrix< real > &aB)
	Inverse of a 2x2 matrix from the closed-form adjugate.
real	belfem::inv3 (const Matrix< real > &aA, Matrix< real > &aB)
	Inverse of a 3x3 matrix from the closed-form adjugate.
template<typename T>
T	belfem::polyval (const Vector< T > &aCoeffs, const T aX)
	Evaluates a polynomial at one point, by Horner's scheme.
template<typename T>
void	belfem::sort (Vector< T > &aVector)
	Sorts a vector in place, in ascending order.
template<typename T>
void	belfem::unique (Vector< T > &aVector)
	Sorts a vector in place and drops duplicate entries.
template<typename T>
auto	belfem::cross (const Vector< T > &aA, const Vector< T > &aB) -> decltype(cross(aA.vector_data(), aB.vector_data()))
	Cross product of two length-3 vectors.
template<typename T>
auto	belfem::det (const Matrix< T > &aA) -> decltype(det(aA.matrix_data())) const
	Determinant of a square matrix.
template<typename T>
auto	belfem::dot (const Vector< T > &aA, const Vector< T > &aB) -> decltype(arma::dot(aA.vector_data(), aB.vector_data()))
	Scalar product of two vectors.
template<typename T>
auto	belfem::dot (const Matrix< T > &aA, const Vector< T > &aB) -> decltype(arma::dot(aA.matrix_data(), aB.vector_data()))
	Scalar product of a matrix and a vector, both read as flat sequences.
template<typename T>
auto	belfem::dot (const Vector< T > &aA, const Matrix< T > &aB) -> decltype(arma::dot(aA.vector_data(), aB.matrix_data()))
	Scalar product of a vector and a matrix, both read as flat sequences.
int_t	belfem::eigen (const Matrix< real > &aMatrix, Vector< real > &aValues, const bool aAbortOnComplex=true)
	Eigenvalues of a general (not necessarily symmetric) square matrix.
void	belfem::eigen_sym (const Matrix< real > &aMatrix, Vector< real > &aValues)
	Eigenvalues of a symmetric matrix.
template<typename T>
auto	belfem::inv (const Matrix< T > &aA) -> decltype(inv(aA.matrix_data())) const
	Inverse of a square matrix.
template<typename T>
auto	belfem::norm (const Vector< T > &aA) -> decltype(norm(aA.vector_data()))
	Euclidean (L2) norm of a vector.

Detailed Description

Backend-agnostic linear algebra API and LAPACK wrappers.

Source: src/linalg. Hand-written documentation: Linear Algebra Module Documentation

Function Documentation

◆ append()

template<typename T>

void belfem::append	(	Vector< T > &	aA,
		Vector< T > &	aB )

Appends one vector onto the end of another, in place.

Warning: aA and aB must be different objects. In assertion-enabled builds the Blaze path checks this with BELFEM_ASSERT; the Armadillo path has no such check, and in a release build neither does. Self-appending would silently produce a doubled vector. Do not rely on either backend's aliasing behaviour.

The overload taking aB by const reference carries exactly the same hazard.

Parameters

aA	grown by the length of `aB`, with `aB's` entries copied onto its end
aB	the vector to append; not modified, despite the non-const reference

◆ cross()

template<typename T>

auto belfem::cross	(	const Vector< T > &	aA,
		const Vector< T > &	aB ) -> decltype(cross(aA.vector_data(), aB.vector_data()))

Cross product of two length-3 vectors.

Parameters

aA	left input vector; must have length 3
aB	right input vector; must have length 3

Returns: the right-hand-rule vector aA x aB, with magnitude equal to the area of the parallelogram they span

◆ crossmat() [1/4]

void belfem::crossmat	(	const Vector< real > &	aN,
		const Matrix< real > &	aA,
		const real	aScale,
		Matrix< real > &	aNxA )

inline

Crosses a normal vector with every column of a matrix, 3D, scaled.

Parameters

aN	normal vector, length 3
aA	matrix whose first three rows hold the vector components; further rows are ignored
aScale	multiplier applied to each column cross product before accumulation
aNxA	accumulated into; 3 rows, one column per column of `aA`; initialise it first

◆ crossmat() [2/4]

void belfem::crossmat	(	const Vector< real > &	aN,
		const Matrix< real > &	aA,
		const real	aScale,
		Vector< real > &	aNxA )

inline

Crosses a normal vector with every column of a matrix, 2D, scaled.

Parameters

aN	normal vector, length 2
aA	matrix whose first two rows hold the vector components; further rows are ignored
aScale	multiplier applied to each column cross product before accumulation
aNxA	accumulated into, one entry per column of `aA`; initialise it first

◆ crossmat() [3/4]

void belfem::crossmat	(	const Vector< real > &	aN,
		const Matrix< real > &	aA,
		Matrix< real > &	aNxA )

inline

Crosses a normal vector with every column of a matrix, 3D.

Parameters

aN	normal vector, length 3
aA	matrix whose first three rows hold the vector components; further rows are ignored
aNxA	assigned to; must have 3 rows and one column per column of `aA`

◆ crossmat() [4/4]

void belfem::crossmat	(	const Vector< real > &	aN,
		const Matrix< real > &	aA,
		Vector< real > &	aNxA )

inline

Crosses a normal vector with every column of a matrix, 2D.

Parameters

aN	normal vector, length 2
aA	matrix whose first two rows hold the vector components; further rows are ignored
aNxA	assigned to, one entry per column of `aA`; must already have that length

◆ ddpolyval()

template<typename T>

T belfem::ddpolyval	(	const Vector< T > &	aCoeffs,
		const T	aX )

Evaluates the second derivative of a polynomial at one point.

Differentiates the polynomial described by aCoeffs twice, analytically; nothing is approximated by finite differences.

Parameters

aCoeffs	coefficients of the polynomial itself, highest power first – not the coefficients of either derivative. Must not be empty.
aX	the point to evaluate at

Returns: the second derivative at aX

◆ det()

template<typename T>

auto belfem::det ( const Matrix< T > & aA ) -> decltype(det(aA.matrix_data())) const

Determinant of a square matrix.

Parameters

aA	square matrix; not modified

Returns: the determinant

◆ dot() [1/3]

template<typename T>

auto belfem::dot	(	const Matrix< T > &	aA,
		const Vector< T > &	aB ) -> decltype(arma::dot(aA.matrix_data(), aB.vector_data()))

Scalar product of a matrix and a vector, both read as flat sequences.

This is not a matrix-vector product: it forwards to the backend's own dot, which walks both operands as flat element sequences and returns a single number. The matrix must therefore hold exactly as many elements as the vector.

Parameters

aA	matrix, read element by element
aB	vector with as many entries as `aA` has elements

Returns: scalar sum of the elementwise products

◆ dot() [2/3]

template<typename T>

auto belfem::dot	(	const Vector< T > &	aA,
		const Matrix< T > &	aB ) -> decltype(arma::dot(aA.vector_data(), aB.matrix_data()))

Scalar product of a vector and a matrix, both read as flat sequences.

The mirror of the overload above, and equally not a vector-matrix product.

Parameters

aA	vector with as many entries as `aB` has elements
aB	matrix, read element by element

Returns: scalar sum of the elementwise products

◆ dot() [3/3]

template<typename T>

auto belfem::dot	(	const Vector< T > &	aA,
		const Vector< T > &	aB ) -> decltype(arma::dot(aA.vector_data(), aB.vector_data()))

Scalar product of two vectors.

Parameters

aA	left input vector; must have the same length as `aB`
aB	right input vector; must have the same length as `aA`

Returns: scalar sum of the elementwise products

◆ dpolyval()

template<typename T>

T belfem::dpolyval	(	const Vector< T > &	aCoeffs,
		const T	aX )

Evaluates the first derivative of a polynomial at one point.

Differentiates the polynomial described by aCoeffs analytically and evaluates the result; nothing is approximated by finite differences.

Parameters

aCoeffs	coefficients of the polynomial itself, highest power first – not the coefficients of its derivative. Must not be empty.
aX	the point to evaluate at

Returns: the first derivative at aX

◆ eigen()

int_t belfem::eigen	(	const Matrix< real > &	aMatrix,
		Vector< real > &	aValues,
		const bool	aAbortOnComplex = true )

inline

Eigenvalues of a general (not necessarily symmetric) square matrix.

A general matrix may have complex eigenvalues, which a real Vector cannot hold. By default that is treated as a caller error and aborts with a BELFEM_ERROR naming the offending index and its imaginary part – an always-active check, not an assert, so it survives release builds.

Pass aAbortOnComplex as false to handle the case instead: the complex entries are then written as BELFEM_QUIET_NAN and the return value counts them. This mirrors the AbortOnError argument of belfem::gesv and belfem::posv.

If the matrix is symmetric, prefer belfem::eigen_sym(), whose eigenvalues are real by construction, so the question cannot arise, and whose solver is faster.

Parameters

aMatrix	square matrix; not modified
aValues	resized to the matrix order and filled with the eigenvalues, in the order the backend returns them, which is not sorted
aAbortOnComplex	abort on the first complex eigenvalue (default); pass false to receive a count instead

Returns: the number of eigenvalues found to be complex; 0 when all are real

Note: An eigenvalue counts as complex when the magnitude of its imaginary part exceeds BELFEM_EPSILON. Both backends use that same threshold.

◆ eigen_sym()

void belfem::eigen_sym	(	const Matrix< real > &	aMatrix,
		Vector< real > &	aValues )

inline

Eigenvalues of a symmetric matrix.

A real symmetric matrix has real eigenvalues, so there is no complex case to handle and no NaN to guard against. Prefer this over belfem::eigen whenever the matrix is known to be symmetric – which covers most of what a finite-element code produces.

The symmetry of aMatrix is not checked. Only the upper triangle is read, on both backends, so passing a non-symmetric matrix yields the eigenvalues of the symmetric matrix implied by that triangle rather than an error.

Parameters

aMatrix	square, symmetric matrix; not modified
aValues	resized to the matrix order and filled with the eigenvalues in ascending order

◆ inv()

template<typename T>

auto belfem::inv ( const Matrix< T > & aA ) -> decltype(inv(aA.matrix_data())) const

Inverse of a square matrix.

Parameters

aA	square matrix; not modified

Returns: the inverse of aA

◆ inv2()

real belfem::inv2	(	const Matrix< real > &	aA,
		Matrix< real > &	aB )

inline

Inverse of a 2x2 matrix from the closed-form adjugate.

Computes the determinant and the adjugate directly, without a backend or LAPACK call. Returning the determinant is deliberate: element Jacobians need it anyway, so it comes for free rather than being recomputed.

The singularity check is relative, not absolute: in debug builds BELFEM_ASSERT requires det^2 > BELFEM_EPSILON^2 times the product of the squared row norms, so the tolerance scales with the magnitude of aA. Release builds compile the check out; a singular or nearly singular matrix then produces invalid or very large floating-point values rather than a diagnostic.

Parameters

aA	the matrix to invert; not modified
aB	filled with the inverse; must already be 2x2

Returns: the determinant of aA

◆ inv3()

real belfem::inv3	(	const Matrix< real > &	aA,
		Matrix< real > &	aB )

inline

Inverse of a 3x3 matrix from the closed-form adjugate.

Computes the determinant and the adjugate directly, without a backend or LAPACK call. Returning the determinant is deliberate: element Jacobians need it anyway, so it comes for free rather than being recomputed.

The singularity check is relative, not absolute: in debug builds BELFEM_ASSERT requires det^2 > BELFEM_EPSILON^2 times the product of the squared row norms, so the tolerance scales with the magnitude of aA. Release builds compile the check out; a singular or nearly singular matrix then produces invalid or very large floating-point values rather than a diagnostic.

Parameters

aA	the matrix to invert; not modified
aB	filled with the inverse; must already be 3x3

Returns: the determinant of aA

◆ norm()

template<typename T>

auto belfem::norm ( const Vector< T > & aA ) -> decltype(norm(aA.vector_data()))

Euclidean (L2) norm of a vector.

Parameters

aA	vector whose Euclidean norm is computed

Returns: the square root of the sum of the squared entries

◆ polyval()

template<typename T>

T belfem::polyval	(	const Vector< T > &	aCoeffs,
		const T	aX )

Evaluates a polynomial at one point, by Horner's scheme.

Parameters

aCoeffs	coefficients, highest power first Must not be empty.
aX	the point to evaluate at

Returns: the value of the polynomial at aX

◆ sort()

template<typename T>

void belfem::sort ( Vector< T > & aVector )

Sorts a vector in place, in ascending order.

Parameters

aVector sorted in place; its length does not change

◆ unique()

template<typename T>

void belfem::unique ( Vector< T > & aVector )

Sorts a vector in place and drops duplicate entries.

The vector is sorted as a side effect, on both backends – this is not a duplicate filter that preserves the original order. The vector shrinks to the number of distinct values.

Parameters

aVector sorted, deduplicated and resized in place

Files

Classes

Functions

Detailed Description

Function Documentation

◆ append()

◆ cross()

◆ crossmat() [1/4]

◆ crossmat() [2/4]

◆ crossmat() [3/4]

◆ crossmat() [4/4]

◆ ddpolyval()

◆ det()

◆ dot() [1/3]

◆ dot() [2/3]

◆ dot() [3/3]

◆ dpolyval()

◆ eigen()

◆ eigen_sym()

◆ inv()

◆ inv2()

◆ inv3()

◆ norm()

◆ polyval()

◆ sort()

◆ unique()