Detailed Description

Generic functions that wrap the BLAS and LAPACK interfaces to provide a more user-friendly approach for the most common linear algebra related tasks.

Functions
template<typename V>
auto	nda::linalg::cross_product (V const &x, V const &y)
	Compute the cross product of two 3-dimensional vectors.
template<typename... A> requires (nda::clef::is_any_lazy<A...>)
auto	nda::clef::determinant (A &&...__a)
	Lazy version of nda::determinant.
template<typename M>
auto	nda::determinant (M const &m)
	Compute the determinant of a square matrix/view.
template<typename M> requires (is_matrix_or_view_v<M>)
auto	nda::determinant_in_place (M &m)
	Compute the determinant of a square matrix/view.
template<typename X, typename Y> requires ((Scalar<X> and Scalar<Y>) or (Vector<X> and Vector<Y>))
auto	nda::linalg::dot (X const &x, Y const &y)
	Compute the dot product of two nda::vector objects or the product of two scalars.
template<bool star = false, Vector X, Vector Y> requires (Scalar<get_value_t<X>> and Scalar<get_value_t<Y>> and mem::have_host_compatible_addr_space<X, Y>)
auto	nda::linalg::dot_generic (X const &x, Y const &y)
	Generic loop-based dot product implementation for vectors.
template<typename X, typename Y> requires ((Scalar<X> and Scalar<Y>) or (Vector<X> and Vector<Y>))
auto	nda::linalg::dotc (X const &x, Y const &y)
	Compute the dotc (LHS operand is conjugated) product of two nda::vector objects or the product of two scalars.
template<typename M>
auto	nda::linalg::eigenelements (M const &m)
	Find the eigenvalues and eigenvectors of a symmetric (real) or hermitian (complex) matrix/view.
template<typename M>
auto	nda::linalg::eigenvalues (M const &m)
	Find the eigenvalues of a symmetric (real) or hermitian (complex) matrix/view.
template<typename M>
auto	nda::linalg::eigenvalues_in_place (M &m)
	Find the eigenvalues of a symmetric (real) or hermitian (complex) matrix/view.
template<typename... A> requires (nda::clef::is_any_lazy<A...>)
auto	nda::clef::inverse (A &&...__a)
	Lazy version of nda::inverse.
template<Matrix M> requires (get_algebra<M> == 'M')
auto	nda::inverse (M const &m)
	Compute the inverse of an n-by-n matrix.
template<MemoryMatrix M> requires (get_algebra<M> == 'M' and mem::on_host<M>)
void	nda::inverse1_in_place (M &&m)
	Compute the inverse of a 1-by-1 matrix.
template<MemoryMatrix M> requires (get_algebra<M> == 'M' and mem::on_host<M>)
void	nda::inverse2_in_place (M &&m)
	Compute the inverse of a 2-by-2 matrix.
template<MemoryMatrix M> requires (get_algebra<M> == 'M' and mem::on_host<M>)
void	nda::inverse3_in_place (M &&m)
	Compute the inverse of a 3-by-3 matrix.
template<MemoryMatrix M> requires (get_algebra<M> == 'M')
void	nda::inverse_in_place (M &&m)
	Compute the inverse of an n-by-n matrix.
template<typename A>
bool	nda::is_matrix_diagonal (A const &a, bool print_error=false)
	Check if a given array/view is diagonal, i.e. if it is square (see nda::is_matrix_square) and all the the off-diagonal elements are zero.
template<typename A>
bool	nda::is_matrix_square (A const &a, bool print_error=false)
	Check if a given array/view is square, i.e. if the first dimension has the same extent as the second dimension.
template<Matrix A, Matrix B>
auto	nda::linalg::matmul (A &&a, B &&b)
	Compute the matrix-matrix product of two nda::matrix objects.
template<Matrix A, Vector X> requires (mem::have_compatible_addr_space<A, X>)
auto	nda::linalg::matvecmul (A const &a, X const &x)
	Compute the matrix-vector product of an nda::matrix and an nda::vector object.
template<MemoryArray A, MemoryArray B> requires ((nda::blas::has_C_layout<A> or nda::blas::has_F_layout<A>) and nda::blas::has_C_layout<A> == nda::blas::has_C_layout<B>)
auto	nda::linalg::outer_product (A const &a, B const &b)
	Outer product of two arrays/views.

Function Documentation

◆ cross_product()

template<typename V>

auto nda::linalg::cross_product	(	V const &	x,
		V const &	y )

#include <nda/linalg/cross_product.hpp>

Compute the cross product of two 3-dimensional vectors.

Template Parameters

V	Vector type.

Parameters

x	Left hand side vector.
y	Right hand side vector.

Returns: nda::array of rank 1 containing the cross product of the two vectors.

Definition at line 29 of file cross_product.hpp.

◆ determinant()

template<typename M>

auto nda::determinant ( M const & m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the determinant of a square matrix/view.

The given matrix/view is not modified. It first makes a copy of the given matrix/view and then calls nda::determinant_in_place with the copy.

Template Parameters

M	Type of the matrix/view.

Parameters

m	Matrix/view object.

Returns: Determinant of the matrix/view.

Definition at line 132 of file det_and_inverse.hpp.

◆ determinant_in_place()

template<typename M>
requires (is_matrix_or_view_v<M>)

auto nda::determinant_in_place ( M & m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the determinant of a square matrix/view.

It uses nda::lapack::getrf to compute the LU decomposition of the matrix and then calculates the determinant by multiplying the diagonal elements of the \( \mathbf{U} \) matrix and taking into account that getrf may change the ordering of the rows/columns of the matrix.

The given matrix/view is modified in place.

Template Parameters

M	Type of the matrix/view.

Parameters

m	Matrix/view object.

Returns: Determinant of the matrix/view.

Definition at line 89 of file det_and_inverse.hpp.

◆ dot()

template<typename X, typename Y>
requires ((Scalar<X> and Scalar<Y>) or (Vector<X> and Vector<Y>))

auto nda::linalg::dot	(	X const &	x,
		Y const &	y )

#include <nda/linalg/dot.hpp>

Compute the dot product of two nda::vector objects or the product of two scalars.

The behaviour of this function is identical to nda::blas::dot, except that it allows

the two input objects to be scalars,
lazy expressions as input vectors,
the value types of the input vectors to be different from each other and
the value types of the input vectors to be different from nda::is_double_or_complex_v.

For vectors, it calls nda::blas::dot if possible, otherwise it simply loops over the input arrays/views and sums up the element-wise products.

Note: The first argument is never conjugated. Even for complex types. Use nda::linalg::dotc for that.

Template Parameters

X	nda::Vector or nda::Scalar type.
Y	nda::Vector or nda::Scalar type.

Parameters

x	Input vector/scalar.
y	Input vector/scalar.

Returns: Result of the dot product.

Definition at line 104 of file dot.hpp.

◆ dot_generic()

template<bool star = false, Vector X, Vector Y>
requires (Scalar<get_value_t<X>> and Scalar<get_value_t<Y>> and mem::have_host_compatible_addr_space<X, Y>)

auto nda::linalg::dot_generic	(	X const &	x,
		Y const &	y )

#include <nda/linalg/dot.hpp>

Generic loop-based dot product implementation for vectors.

Computes the dot product of two vector objects with optional conjugation:

For star = false: result = sum(x[i] * y[i])
For star = true: result = sum(conj(x[i]) * y[i]) for complex types

Template Parameters

star	If true, conjugate the first operand (for complex types only).
X	Vector type.
Y	Vector type.

Parameters

x	First input vector.
y	Second input vector.

Returns: The computed dot product.

Definition at line 45 of file dot.hpp.

◆ dotc()

template<typename X, typename Y>
requires ((Scalar<X> and Scalar<Y>) or (Vector<X> and Vector<Y>))

auto nda::linalg::dotc	(	X const &	x,
		Y const &	y )

#include <nda/linalg/dot.hpp>

Compute the dotc (LHS operand is conjugated) product of two nda::vector objects or the product of two scalars.

The behaviour of this function is identical to nda::blas::dotc, except that it allows

the two input objects to be scalars,
lazy expressions as input vectors,
the value types of the input vectors to be different from each other and
the value types of the input vectors to be different from nda::is_double_or_complex_v.

For vectors, it calls nda::blas::dotc if possible, otherwise it simply loops over the input arrays/views and sums up the element-wise products with the LHS operand conjugated.

Template Parameters

X	nda::Vector or nda::Scalar type.
Y	nda::Vector or nda::Scalar type.

Parameters

x	Input vector/scalar.
y	Input vector/scalar.

Returns: Result of the dotc product.

Definition at line 135 of file dot.hpp.

◆ eigenelements()

template<typename M>

auto nda::linalg::eigenelements ( M const & m )

#include <nda/linalg/eigenelements.hpp>

Find the eigenvalues and eigenvectors of a symmetric (real) or hermitian (complex) matrix/view.

For a real symmetric matrix/view, it calls the LAPACK routine syev. For a complex hermitian matrix/view, it calls the LAPACK routine heev.

The given matrix/view is copied and the original is not modified.

Template Parameters

M	Type of the matrix/view.

Parameters

m	Matrix/View to diagonalize.

Returns: std::pair consisting of the array of eigenvalues in ascending order and the matrix containing the eigenvectors in its columns.

Definition at line 87 of file eigenelements.hpp.

◆ eigenvalues()

template<typename M>

auto nda::linalg::eigenvalues ( M const & m )

#include <nda/linalg/eigenelements.hpp>

Find the eigenvalues of a symmetric (real) or hermitian (complex) matrix/view.

For a real symmetric matrix/view, it calls the LAPACK routine syev. For a complex hermitian matrix/view, it calls the LAPACK routine heev.

The given matrix/view is copied and the original is not modified.

Template Parameters

M	Type of the matrix/view.

Parameters

m	Matrix/View to diagonalize.

Returns: An nda::array of rank 1 containing the eigenvalues in ascending order.

Definition at line 106 of file eigenelements.hpp.

◆ eigenvalues_in_place()

template<typename M>

auto nda::linalg::eigenvalues_in_place ( M & m )

#include <nda/linalg/eigenelements.hpp>

Find the eigenvalues of a symmetric (real) or hermitian (complex) matrix/view.

For a real symmetric matrix/view, it calls the LAPACK routine syev. For a complex hermitian matrix/view, it calls the LAPACK routine heev.

The given matrix/view will be modified by the diagonalization process.

Template Parameters

M	Type of the matrix/view.

Parameters

m	Matrix/View to diagonalize.

Returns: An nda::array of rank 1 containing the eigenvalues in ascending order.

Definition at line 124 of file eigenelements.hpp.

◆ inverse()

template<Matrix M>
requires (get_algebra<M> == 'M')

auto nda::inverse ( M const & m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the inverse of an n-by-n matrix.

The given matrix/view is not modified. It first makes copy of the given matrix/view and then calls nda::inverse_in_place with the copy.

Template Parameters

M	nda::MemoryMatrix type.

Parameters

m	nda::MemoryMatrix object to be inverted.

Returns: Inverse of the matrix.

Definition at line 286 of file det_and_inverse.hpp.

◆ inverse1_in_place()

template<MemoryMatrix M>
requires (get_algebra<M> == 'M' and mem::on_host<M>)

void nda::inverse1_in_place ( M && m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the inverse of a 1-by-1 matrix.

The inversion is performed in place.

Template Parameters

M	nda::MemoryMatrix type.

Parameters

m	nda::MemoryMatrix object to be inverted.

Definition at line 158 of file det_and_inverse.hpp.

◆ inverse2_in_place()

template<MemoryMatrix M>
requires (get_algebra<M> == 'M' and mem::on_host<M>)

void nda::inverse2_in_place ( M && m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the inverse of a 2-by-2 matrix.

The inversion is performed in place.

Template Parameters

M	nda::MemoryMatrix type.

Parameters

m	nda::MemoryMatrix object to be inverted.

Definition at line 173 of file det_and_inverse.hpp.

◆ inverse3_in_place()

template<MemoryMatrix M>
requires (get_algebra<M> == 'M' and mem::on_host<M>)

void nda::inverse3_in_place ( M && m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the inverse of a 3-by-3 matrix.

The inversion is performed in place.

Template Parameters

M	nda::MemoryMatrix type.

Parameters

m	nda::MemoryMatrix object to be inverted.

Definition at line 199 of file det_and_inverse.hpp.

◆ inverse_in_place()

template<MemoryMatrix M>
requires (get_algebra<M> == 'M')

void nda::inverse_in_place ( M && m )

#include <nda/linalg/det_and_inverse.hpp>

Compute the inverse of an n-by-n matrix.

The inversion is performed in place.

For small matrices (1-by-1, 2-by-2, 3-by-3), we directly compute the matrix inversion using the optimized routines: nda::inverse1_in_place, nda::inverse2_in_place, nda::inverse3_in_place.

For larger matrices, it uses nda::lapack::getrf and nda::lapack::getri.

Template Parameters

M	nda::MemoryMatrix type.

Parameters

m	nda::MemoryMatrix object to be inverted.

Definition at line 243 of file det_and_inverse.hpp.

◆ is_matrix_diagonal()

template<typename A>

bool nda::is_matrix_diagonal	(	A const &	a,
		bool	print_error = false )

#include <nda/linalg/det_and_inverse.hpp>

Check if a given array/view is diagonal, i.e. if it is square (see nda::is_matrix_square) and all the the off-diagonal elements are zero.

Note: It does not check if the array/view has rank 2.

Template Parameters

A	Array/View type.

Parameters

a	Array/View object.
print_error	If true, print an error message if the matrix is not diagonal.

Returns: True if the array/view is diagonal, false otherwise.

Definition at line 69 of file det_and_inverse.hpp.

◆ is_matrix_square()

template<typename A>

bool nda::is_matrix_square	(	A const &	a,
		bool	print_error = false )

#include <nda/linalg/det_and_inverse.hpp>

Check if a given array/view is square, i.e. if the first dimension has the same extent as the second dimension.

Note: It does not check if the array/view has rank 2.

Template Parameters

A	Array/View type.

Parameters

a	Array/View object.
print_error	If true, print an error message if the matrix is not square.

Returns: True if the array/view is square, false otherwise.

Definition at line 50 of file det_and_inverse.hpp.

◆ matmul()

template<Matrix A, Matrix B>

auto nda::linalg::matmul	(	A &&	a,
		B &&	b )

#include <nda/linalg/matmul.hpp>

Compute the matrix-matrix product of two nda::matrix objects.

This function computes the matrix-matrix product \( \mathrm{op}_1(\mathbf{A}) \mathrm{op}_2(\mathbf{B}) \), where \( \mathbf{A} \) and \( \mathbf{B} \) are \( m \times k \) and \( k \times n \) matrices respectively, \( \mathrm{op}_i \) can be some lazy operation, e.g. nda::conj, nda::sin, etc.

We try to call nda::blas::gemm whenever possible, i.e. when the value type of the result is compatible with nda::is_blas_lapack_v, even if this requires to make copies of the input arrays/views. Otherwise, we perform a very naive and inefficient matrix-matrix multiplication manually.

Therefore, if performance is important, users should make sure to pass input arrays/views which are compatible with nda::blas::gemm.

Note: The layout of the returned matrix depends on the layout of the input matrices. If both input matrices are in Fortran-layout, the returned matrix is also in Fortran layout. Otherwise, it is in C-layout.

Warning: This function might make copies of the input arrays/views. When working on the device memory space, this may lead to runtime errors if the copying fails.

Template Parameters

A	nda::Matrix type.
B	nda::Matrix type.

Parameters

a	Input matrix \( \mathbf{A} \) of size \( m \times k \).
b	Input matrix \( \mathbf{B} \) of size \( k \times n \).

Returns: Resulting matrix of the matrix-matrix multiplication of size \( m \times n \).

Definition at line 124 of file matmul.hpp.

◆ matvecmul()

template<Matrix A, Vector X>
requires (mem::have_compatible_addr_space<A, X>)

auto nda::linalg::matvecmul	(	A const &	a,
		X const &	x )

#include <nda/linalg/matvecmul.hpp>

Compute the matrix-vector product of an nda::matrix and an nda::vector object.

This function computes the matrix-vector product \( \mathrm{op}_1(\mathbf{A}) \mathrm{op}_2(\mathbf{x}) \), where \( \mathbf{A} \) is an \( m \times n \) matrix, \( \mathbf{x} \) is a vector of size \( n \), and \( \mathrm{op}_i \) can be some lazy operation, e.g. nda::conj, nda::sin, etc.

We try to call nda::blas::gemv whenever possible, i.e. when the value type of the result is compatible with nda::is_blas_lapack_v, even if this requires to make copies of the input arrays/views. Otherwise, we perform a very naive and inefficient matrix-vector multiplication manually.

Therefore, if performance is important, users should make sure to pass input arrays/views which are compatible with nda::blas::gemv.

Warning: This function might make copies of the input arrays/views. When working on the device memory space, this may lead to runtime errors if the copying fails.

Template Parameters

A	nda::Matrix type.
X	nda::Vector type.

Parameters

a	Input matrix \( \mathbf{A} \) of size \( m \times n \).
x	Input vector \( \mathbf{x} \) of size \( n \).

Returns: Resulting vector of the matrix-vector multiplication of size \( m \).

Definition at line 116 of file matvecmul.hpp.

◆ outer_product()

template<MemoryArray A, MemoryArray B>
requires ((nda::blas::has_C_layout<A> or nda::blas::has_F_layout<A>) and nda::blas::has_C_layout<A> == nda::blas::has_C_layout<B>)

auto nda::linalg::outer_product	(	A const &	a,
		B const &	b )

#include <nda/linalg/outer_product.hpp>

Outer product of two arrays/views.

It calculates the outer product \( \mathbf{C} = \mathbf{A} \otimes \mathbf{B} \), such that

\[ \mathbf{C}_{i_1 \ldots i_k j_1 \ldots j_l} = \mathbf{A}_{i_1 \ldots i_k} \mathbf{B}_{j_1 \ldots j_l} \; . \]

Here, \( \mathbf{A} \) and \( \mathbf{B} \) are the input arrays with shape \( (m_1, \ldots, m_k) \) and \((n_1, \ldots, n_l) \), respectively. The resulting array \( \mathbf{C} \) has shape \( (m_1, \ldots, m_k, n_1, \ldots, n_l) \).

The outer product is performed by calling nda::blas::ger, which imposes various constraints on the supported input arrays/views, e.g.

their memory layouts have to be the same and either C- or Fortran,
they have to be contiguous in memory,
etc.

See nda::blas::ger for more details.

Template Parameters

A	nda::MemoryArray type.
B	nda::MemoryArray type.

Parameters

a	Input array/view \( \mathbf{A} \).
b	Input array/view \( \mathbf{B} \).

Returns: Outer product \( \mathbf{A} \otimes \mathbf{B} \).

Definition at line 52 of file outer_product.hpp.

Detailed Description

Functions

Function Documentation

◆ cross_product()

◆ determinant()

◆ determinant_in_place()

◆ dot()

◆ dot_generic()

◆ dotc()

◆ eigenelements()

◆ eigenvalues()

◆ eigenvalues_in_place()

◆ inverse()

◆ inverse1_in_place()

◆ inverse2_in_place()

◆ inverse3_in_place()

◆ inverse_in_place()

◆ is_matrix_diagonal()

◆ is_matrix_square()

◆ matmul()

◆ matvecmul()

◆ outer_product()