Detailed Description

Linear algebra: matrices, vectors, views, and solvers.

Provides matrix and vector types, non-owning views, and a set of linear algebra solvers wrapping LAPACK and GSL.

Core types:

LinAlg::Matrix<T> : owning row-major matrix; supports real and complex element types. Provides arithmetic, GSL interop, and element-wise access.
LinAlg::Vector<T> : owning 1D array; used for eigenvectors, right-hand sides, and eigenvalue arrays.

Views:

LinAlg::View<T> : non-owning strided view over a 1D segment of an array. Used to provide row and column access into a Matrix without copying. Obtained via Matrix::row_view() and Matrix::column_view()
LinAlg::Matrix_view<T> : non-owning view of a 2D subblock of a Matrix.

Solvers:

LinAlg::solve_Axeqb() : solves Ax = b via LU decomposition (GSL).
LinAlg::symmhEigensystem() : generalised problem Av = eBv or Av = ev Several overloads, and templayte for real/complex types (LAPACK dsygv/zhegv/dsyevr/dsyev/zheev).
LinAlg::genEigensystem() : general non-symmetric real matrix (GSL).

Classes
class	Matrix
	Row-major dense matrix with arithmetic and linear algebra support. More...

class	Matrix_view
	Non-owning 2D view onto a Matrix. More...

class	Vector
	Owning 1D array; inherits from Matrix<T> with a single column. More...

class	View
	Non-owning strided view onto a 1D segment of an array. More...

Functions
CBLAS_TRANSPOSE	to_cblas_trans (bool trans)
	Converts bool to CBLAS_TRANSPOSE enum (CblasTrans if true, CblasNoTrans if false)

template<typename T >
void	GEMM (const Matrix< T > &a, const Matrix< T > &b, Matrix< T > *c, bool trans_A=false, bool trans_B=false)
	Matrix multiplication C = op(A) * op(B) via CBLAS GEMM (row-major).

template<typename T >
void	PENTA_GEMM (const Matrix< T > &A, const Matrix< T > &B, const Matrix< T > &C, const Matrix< T > &D, const Matrix< T > &E, Matrix< T > *M)
	5-matrix contraction for N*N matrices: M_ab = A_ai B_aj C_ij D_ib E_jb, with BLAS

template<typename T >
Matrix< T >	PENTA_GEMM (const Matrix< T > &A, const Matrix< T > &B, const Matrix< T > &C, const Matrix< T > &D, const Matrix< T > &E)
	M_ab = A_ai B_aj C_ij D_ib E_jb. Assuming all square matrices. Uses blas.

template<typename T , bool PARALLEL = false>
void	PENTA (const Matrix< T > &A, const Matrix< T > &B, const Matrix< T > &C, const Matrix< T > &D, const Matrix< T > &E, Matrix< T > *M)
	5-matrix contraction for N*N matrices: M_ab = A_ai B_aj C_ij D_ib E_jb, without BLAS, but with optional OpenMP parallelisation

template<typename T , bool PARALLEL = false>
Matrix< T >	PENTA (const Matrix< T > &A, const Matrix< T > &B, const Matrix< T > &C, const Matrix< T > &D, const Matrix< T > &E)
	5-index contraction for N*N matrices: M_ab = A_ai B_aj C_ij D_ib E_jb, without BLAS, but with optional OpenMP parallelisation

template<typename T >
constexpr auto	myEps ()
	Default relative tolerance for `equal()`: 1e-6 for float, 1e-12 for double.

template<typename T >
bool	equal (const Matrix< T > &lhs, const Matrix< T > &rhs, T eps=myEps< T >())
	Compares two matrices element-wise to within a relative tolerance.

template<typename T >
Matrix< T >	operator* (const Matrix< T > &a, const Matrix< T > &b)

template<typename T >
std::ostream &	operator<< (std::ostream &os, const Matrix< T > &a)

template<typename T >
Vector< T >	solve_Axeqb (Matrix< T > Am, const Vector< T > &b)
	Solves the linear system Ax = b for x, given square matrix A and vector b.

template<typename T >
std::pair< Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A)
	Solves Av = ev for all eigenvalues and eigenvectors of a real symmetric or complex Hermitian matrix A.

template<typename T >
std::pair< Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A, int number)
	Solves Av = ev for the first `number` eigenvalues and eigenvectors of a real symmetric matrix A.

template<typename T >
std::tuple< int, Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A, double all_below)
	Solves Av = ev for all eigenvalues below a given threshold of a real symmetric matrix A.

template<typename T >
std::pair< Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A, Matrix< T > B)
	Solves the generalised eigensystem Av = eBv for a real symmetric or complex Hermitian matrix pair A, B.

template<typename T >
std::pair< Vector< std::complex< double > >, Matrix< std::complex< double > > >	genEigensystem (Matrix< T > A, bool sort)
	Solves Av = ev for all eigenvalues and right eigenvectors of a general (non-symmetric) real matrix A.

Function Documentation

◆ to_cblas_trans()

CBLAS_TRANSPOSE LinAlg::to_cblas_trans ( bool trans )

inline

Converts bool to CBLAS_TRANSPOSE enum (CblasTrans if true, CblasNoTrans if false)

◆ GEMM()

template<typename T >

void LinAlg::GEMM	(	const Matrix< T > &	a,
		const Matrix< T > &	b,
		Matrix< T > *	c,
		bool	trans_A = `false`,
		bool	trans_B = `false`
	)

Matrix multiplication C = op(A) * op(B) via CBLAS GEMM (row-major).

Computes \( C = \mathrm{op}(A)\,\mathrm{op}(B) \) where \( \mathrm{op}(X) = X \) or \( X^T \) depending on trans_A / trans_B. c is overwritten with the result.

Wraps CBLAS gemm for double, float, std::complex<double>, and std::complex<float>.

Parameters

	a	Left-hand matrix.
	b	Right-hand matrix.
[out]	c	Output matrix; must be pre-allocated to the correct size.
	trans_A	If true, use \(A^T\) in the product.
	trans_B	If true, use \(B^T\) in the product.

Warning: c must not alias a or b.

◆ PENTA_GEMM() [1/2]

template<typename T >

void LinAlg::PENTA_GEMM	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		const Matrix< T > &	C,
		const Matrix< T > &	D,
		const Matrix< T > &	E,
		Matrix< T > *	M
	)

5-matrix contraction for N*N matrices: M_ab = A_ai B_aj C_ij D_ib E_jb, with BLAS

All input matrices are assumed square with the same dimension N. The contraction sums over indices i and j:

\[ M_{ab} = \sum_{i,j} A_{ai}\,B_{aj}\,C_{ij}\,D_{ib}\,E_{jb}. \]

This implementation evaluates the contraction by looping over the shared index i and using a BLAS GEMM for the inner j-summation:

\[ X^{(i)}_{jb} = C_{ij}\,E_{jb}, \qquad Y^{(i)}_{ab} = \sum_j B_{aj}\,X^{(i)}_{jb}, \]

followed by the accumulation

\[ M_{ab} \mathrel{+}= A_{ai}\,Y^{(i)}_{ab}\,D_{ib}. \]

Thus each i-iteration performs one dense matrix multiplication (B * X^{(i)}) via BLAS, while the remaining operations are elementwise scalings and accumulations.

The routine allocates temporary work matrices X and Y of size NN (reused across i) and writes the result into M, which must be preallocated and is overwritten on entry (i.e., initialised to zero).

Note: ~20% faster than Naiive implementation, BUT not easily parallelisable.

Template Parameters

T	Scalar type (float, double, std::complex<...>)

Parameters

	A,B,C,D,E	Input NN matrices
[out]	M	Output NN matrix receiving the contraction. Must be previously sized

◆ PENTA_GEMM() [2/2]

template<typename T >

Matrix< T > LinAlg::PENTA_GEMM	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		const Matrix< T > &	C,
		const Matrix< T > &	D,
		const Matrix< T > &	E
	)

M_ab = A_ai B_aj C_ij D_ib E_jb. Assuming all square matrices. Uses blas.

Calls: PENTA_GEMM(const Matrix<T>&, const Matrix<T>&, const Matrix<T>&, const Matrix<T>&, const Matrix<T>&, Matrix<T>*).

◆ PENTA() [1/2]

template<typename T , bool PARALLEL = false>

void LinAlg::PENTA	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		const Matrix< T > &	C,
		const Matrix< T > &	D,
		const Matrix< T > &	E,
		Matrix< T > *	M
	)

5-matrix contraction for N*N matrices: M_ab = A_ai B_aj C_ij D_ib E_jb, without BLAS, but with optional OpenMP parallelisation

Computes the tensor contraction

\[ M_{ab} = \sum_{i,j} A_{ai}\,B_{aj}\,C_{ij}\,D_{ib}\,E_{jb}, \]

assuming all input matrices are square with the same dimension.

This implementation uses explicit nested loops (no BLAS). If

Template Parameters

PARALLEL is true, OpenMP pragmas are enabled to parallelise the outer loops.

For a BLAS-accelerated implementation, see PENTA_GEMM(const Matrix<T>&, const Matrix<T>&, const Matrix<T>&, const Matrix<T>&, const Matrix<T>&, Matrix<T>*).

Usually, this version will be significantly faster than the BLAS one. However, if it iself is being called in a parallel region, the BLAS one will be faster.

Template Parameters

T	Scalar type: double, float, std::complex
PARALLEL	Enable OpenMP parallelisation when true (compile-time)

Parameters

	A,B,C,D,E	Input NN matrices
[out]	M	Output NN matrix receiving the contraction. Must be previously sized

◆ PENTA() [2/2]

template<typename T , bool PARALLEL = false>

Matrix< T > LinAlg::PENTA	(	const Matrix< T > &	A,
		const Matrix< T > &	B,
		const Matrix< T > &	C,
		const Matrix< T > &	D,
		const Matrix< T > &	E
	)

5-index contraction for N*N matrices: M_ab = A_ai B_aj C_ij D_ib E_jb, without BLAS, but with optional OpenMP parallelisation

calls above

◆ myEps()

template<typename T >

constexpr auto LinAlg::myEps ( )

constexpr

Default relative tolerance for equal(): 1e-6 for float, 1e-12 for double.

◆ equal()

template<typename T >

bool LinAlg::equal	(	const Matrix< T > &	lhs,
		const Matrix< T > &	rhs,
		T	eps = `myEps<T>()`
	)

Compares two matrices element-wise to within a relative tolerance.

Returns true iff all elements satisfy

\[ |A_{ij} - B_{ij}| \leq |\varepsilon\,(A_{ij} + B_{ij})|. \]

Matrices of different dimensions compare as not equal.

Parameters

lhs	Left matrix.
rhs	Right matrix.
eps	Relative tolerance (default: `myEps<T>()`).

Returns: true if all elements agree within eps; false otherwise.

◆ solve_Axeqb()

template<typename T >

Vector< T > LinAlg::solve_Axeqb	(	Matrix< T >	Am,
		const Vector< T > &	b
	)

Solves the linear system Ax = b for x, given square matrix A and vector b.

Uses LU decomposition (GSL). A is taken by value and overwritten during decomposition.

Template Parameters

T	Element type: `double` or `std::complex<double>` only.

Parameters

Am	Square matrix A (copied; overwritten internally).
b	Right-hand side vector; must satisfy `b.size() == A.rows()`.

Returns: Solution vector x.

◆ symmhEigensystem() [1/4]

template<typename T >

std::pair< Vector< double >, Matrix< T > > LinAlg::symmhEigensystem ( Matrix< T > A )

Solves Av = ev for all eigenvalues and eigenvectors of a real symmetric or complex Hermitian matrix A.

Uses LAPACK dsyev (real) or zheev (complex Hermitian). Eigenvalues are returned in ascending order. A is taken by value and overwritten during the computation.

Template Parameters

T	Element type: `double` (symmetric) or `std::complex<double>` (Hermitian).

Parameters

A	Square symmetric/Hermitian matrix (taken by value; overwritten).

Returns

Pair {e, V}, where:

e(i) is the i-th eigenvalue (always real, ascending order),
V(i,j) is the j-th element of the i-th eigenvector.

◆ symmhEigensystem() [2/4]

template<typename T >

std::pair< Vector< double >, Matrix< T > > LinAlg::symmhEigensystem	(	Matrix< T >	A,
		int	number
	)

Solves Av = ev for the first number eigenvalues and eigenvectors of a real symmetric matrix A.

Uses LAPACK dsyevx. Only available for T = double.

Template Parameters

T	Element type: `double` only.

Parameters

A	Square real symmetric matrix (taken by value; overwritten).
number	Number of lowest eigenvalues/eigenvectors to compute. Must satisfy `number < A.rows()`.

Returns: Pair {e, V}, where only the first number entries of e and rows of V are valid.

◆ symmhEigensystem() [3/4]

template<typename T >

std::tuple< int, Vector< double >, Matrix< T > > LinAlg::symmhEigensystem	(	Matrix< T >	A,
		double	all_below
	)

Solves Av = ev for all eigenvalues below a given threshold of a real symmetric matrix A.

Uses LAPACK dsyevr. Only available for T = double.

Template Parameters

T	Element type: `double` only.

Parameters

A	Square real symmetric matrix (taken by value; overwritten).
all_below	Upper bound: only eigenvalues `all_below` are returned.

Returns

Tuple {N, e, V}, where:

N is the number of eigenvalues found,
e(i) is the i-th eigenvalue (ascending),
V(i,j) is the j-th element of the i-th eigenvector. Only the first N entries of e and rows of V are valid.

◆ symmhEigensystem() [4/4]

template<typename T >

std::pair< Vector< double >, Matrix< T > > LinAlg::symmhEigensystem	(	Matrix< T >	A,
		Matrix< T >	B
	)

Solves the generalised eigensystem Av = eBv for a real symmetric or complex Hermitian matrix pair A, B.

Uses LAPACK dsygv (real) or zhegv (complex Hermitian). B must be positive definite. Eigenvalues are returned in ascending order. Both A and B are taken by value and overwritten during computation.

Template Parameters

T	Element type: `double` or `std::complex<double>`.

Parameters

A	Square symmetric/Hermitian matrix (taken by value; overwritten).
B	Square symmetric/Hermitian positive-definite matrix (same size as A; taken by value; overwritten).

Returns

Pair {e, V}, where:

e(i) is the i-th eigenvalue (always real, ascending order),
V(i,j) is the j-th element of the i-th eigenvector.

◆ genEigensystem()

template<typename T >

std::pair< Vector< std::complex< double > >, Matrix< std::complex< double > > > LinAlg::genEigensystem	(	Matrix< T >	A,
		bool	sort
	)

Solves Av = ev for all eigenvalues and right eigenvectors of a general (non-symmetric) real matrix A.

Uses GSL gsl_eigen_nonsymmv. A must be real (T = double); eigenvalues and eigenvectors are always complex.

Template Parameters

T	Element type: `double` only.

Parameters

A	Square real matrix (taken by value; overwritten).
sort	If `true`, eigenvalues (and corresponding eigenvectors) are sorted by ascending absolute value.

Returns

Pair {e, V}, where:

e(i) is the i-th complex eigenvalue,
V(i,j) is the j-th element of the i-th complex eigenvector.

Detailed Description

Core types:

Views:

Solvers:

Classes

Functions

Function Documentation

◆ to_cblas_trans()

◆ GEMM()

◆ PENTA_GEMM() [1/2]

◆ PENTA_GEMM() [2/2]

◆ PENTA() [1/2]

◆ PENTA() [2/2]

◆ myEps()

◆ equal()

◆ solve_Axeqb()

◆ symmhEigensystem() [1/4]

◆ symmhEigensystem() [2/4]

◆ symmhEigensystem() [3/4]

◆ symmhEigensystem() [4/4]

◆ genEigensystem()