LinAlg Namespace Reference

Defines Matrix, Vector classes, and linear some algebra functions. More...

Classes
class	Matrix
	Matrix class; row-major. More...
class	Vector
	Vector class (inherits from Matrix) More...
class	View
	Proved a "view" onto an array. More...

Functions
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 = A*B. C is overwritten with product, and must be correct size already.
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 ()
template<typename T >
bool	equal (const Matrix< T > &lhs, const Matrix< T > &rhs, T eps=myEps< T >())
	Compares two matrices; returns true iff all elements compare relatively to better than eps.
template<typename T >
Matrix< T >	operator* (const Matrix< T > &a, const Matrix< T > &b)
CBLAS_TRANSPOSE	to_cblas_trans (bool trans)
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 matrix equation Ax=b for x, for known square matrix A and vector b.
template<typename T >
std::pair< Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A)
	Solves Av = ev for eigenvalues e and eigenvectors v of symmetric/Hermetian matrix A. Returns [e,v], where v(i,j) is the jth element of the ith eigenvector corresponding to ith eigenvalue, e(i). e is always real.
template<typename T >
std::pair< Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A, int number)
	Solves Av=ev for first N eigenvals/vecs v of 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 eigenvals of symmetric matrix A below 'all_below'. Returns {N,e,V}: N is number of solutions.
template<typename T >
std::pair< Vector< double >, Matrix< T > >	symmhEigensystem (Matrix< T > A, Matrix< T > B)
	Solves Av = eBv for eigenvalues e and eigenvectors v of symmetric/Hermetian matrix pair A,B. Returns [e,v], where v(i,j) is the jth element of the ith eigenvector corresponding to ith eigenvalue, e(i). e is always real.
template<typename T >
std::pair< Vector< std::complex< double > >, Matrix< std::complex< double > > >	genEigensystem (Matrix< T > A, bool sort)
	Solves Av = ev for eigenvalues e and eigenvectors v of non-symmetric real matrix A. Returns [e,v], where v(i,j) is the jth element of the ith eigenvector corresponding to ith eigenvalue, e(i). A must be real, while e and v are complex.

Detailed Description

Defines Matrix, Vector classes, and linear some algebra functions.

Function Documentation

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 = A*B. C is overwritten with product, and must be correct size already.

Wrapper for BLAS gemm functions:

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

equal()

template<typename T >

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

Compares two matrices; returns true iff all elements compare relatively to better than eps.

solve_Axeqb()

template<typename T >

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

Solves matrix equation Ax=b for x, for known square matrix A and vector b.

symmhEigensystem() [1/4]

template<typename T >

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

(

Matrix< T >

A

)

Solves Av = ev for eigenvalues e and eigenvectors v of symmetric/Hermetian matrix A. Returns [e,v], where v(i,j) is the jth element of the ith eigenvector corresponding to ith eigenvalue, e(i). e is always real.

symmhEigensystem() [2/4]

template<typename T >

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

Solves Av=ev for first N eigenvals/vecs v of symmetric matrix A.

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 eigenvals of symmetric matrix A below 'all_below'. Returns {N,e,V}: N is number of solutions.

symmhEigensystem() [4/4]

template<typename T >

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

Solves Av = eBv for eigenvalues e and eigenvectors v of symmetric/Hermetian matrix pair A,B. Returns [e,v], where v(i,j) is the jth element of the ith eigenvector corresponding to ith eigenvalue, e(i). e is always real.

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 eigenvalues e and eigenvectors v of non-symmetric real matrix A. Returns [e,v], where v(i,j) is the jth element of the ith eigenvector corresponding to ith eigenvalue, e(i). A must be real, while e and v are complex.