TRIQS/nda 1.3.0
Multi-dimensional array library for C++
Loading...
Searching...
No Matches
Linear algebra tools

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::matmul (A &&a, B &&b)
 Perform a matrix-matrix multiplication.
template<Matrix A, Vector X>
auto nda::matvecmul (A &&a, X &&x)
 Perform a matrix-vector multiplication.
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
VVector type.
Parameters
xLeft hand side vector.
yRight 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
MType of the matrix/view.
Parameters
mMatrix/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
MType of the matrix/view.
Parameters
mMatrix/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
Xnda::Vector or nda::Scalar type.
Ynda::Vector or nda::Scalar type.
Parameters
xInput vector/scalar.
yInput 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
starIf true, conjugate the first operand (for complex types only).
XVector type.
YVector type.
Parameters
xFirst input vector.
ySecond 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
Xnda::Vector or nda::Scalar type.
Ynda::Vector or nda::Scalar type.
Parameters
xInput vector/scalar.
yInput 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
MType of the matrix/view.
Parameters
mMatrix/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
MType of the matrix/view.
Parameters
mMatrix/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
MType of the matrix/view.
Parameters
mMatrix/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
Mnda::MemoryMatrix type.
Parameters
mnda::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
Mnda::MemoryMatrix type.
Parameters
mnda::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
Mnda::MemoryMatrix type.
Parameters
mnda::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
Mnda::MemoryMatrix type.
Parameters
mnda::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
Mnda::MemoryMatrix type.
Parameters
mnda::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
AArray/View type.
Parameters
aArray/View object.
print_errorIf 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
AArray/View type.
Parameters
aArray/View object.
print_errorIf 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::matmul ( A && a,
B && b )

#include <nda/linalg/matmul.hpp>

Perform a matrix-matrix multiplication.

It is generic in the sense that it allows the input matrices to belong to a different nda::mem::AddressSpace (as long as they are compatible).

If possible, it uses nda::blas::gemm, otherwise it calls nda::blas::gemm_generic.

Template Parameters
Anda::Matrix type of lhs operand.
Bnda::Matrix type of rhs operand.
Parameters
aLeft hand side matrix operand.
bRight hand side matrix operand.
Returns
Result of the matrix-matrix multiplication.

Definition at line 76 of file matmul.hpp.

◆ matvecmul()

template<Matrix A, Vector X>
auto nda::matvecmul ( A && a,
X && x )

#include <nda/linalg/matmul.hpp>

Perform a matrix-vector multiplication.

It is generic in the sense that it allows the input matrix and vector to belong to a different nda::mem::AddressSpace (as long as they are compatible).

If possible, it uses nda::blas::gemv, otherwise it calls nda::blas::gemv_generic.

Template Parameters
Anda::Matrix type of lhs operand.
Xnda::Vector type of rhs operand.
Parameters
aLeft hand side matrix operand.
xRight hand side vector operand.
Returns
Result of the matrix-vector multiplication.

Definition at line 141 of file matmul.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
Anda::MemoryArray type.
Bnda::MemoryArray type.
Parameters
aInput array/view \( \mathbf{A} \).
bInput array/view \( \mathbf{B} \).
Returns
Outer product \( \mathbf{A} \otimes \mathbf{B} \).

Definition at line 52 of file outer_product.hpp.