Linear Algebra (BLAS, LAPACK, etc.)#

This module contains functions for performing high-level linear algebra operations on tensors.

The functions in this module are designed to be used with the Tensor class and its subclasses.

BLAS Level 1 Functions#

BLAS Level 1 routines perform operations of both addition and reduction on vectors of data. Typical operations include scaling and dot products.

template<TensorConcept AType, TensorConcept BType>
typename AType::data_type einsums::linear_algebra::dot(const AType &A, const BType &B)#

Performs the dot product between two tensors.

This performs \(\sum_{ijk\cdots} A_{ijk\cdots}B_{ijk\cdots}\)

Parameters:

  • A – One of the tensors

  • B – The other tensor.

template<TensorConcept AType, TensorConcept BType, TensorConcept CType>
typename AType::data_type einsums::linear_algebra::dot(const AType &A, const BType &B, const CType &C)#

Performs the dot product between three tensors.

This performs \(\sum_{ijk\cdots} A_{ijk\cdots}B_{ijk\cdots}C_{ijk\cdots}\)

Parameters:

  • A – One of the tensors.

  • B – The second tensor.

  • C – The third tensor.

template<TensorConcept AType>
void einsums::linear_algebra::scale(typename AType::data_type scale, AType *A)#

Scales a tensor by a scalar. 

auto A = einsums::create_ones_tensor("A", 3, 3);

// A is filled with 1.0
einsums::linear_algebra::scale(2.0, &A);
// A is now filled with 2.0
Template Parameters:

  • AType – The type of the tensor

  • ARank – The rank of the tensor

  • T – The underlying data type

Parameters:

  • scale – The scalar to scale the tensor by

  • A – The tensor to scale

template<MatrixConcept AType>
void einsums::linear_algebra::scale_column(size_t col, typename AType::data_type scale, AType *A)#
template<MatrixConcept AType>
void einsums::linear_algebra::scale_row(size_t row, typename AType::data_type scale, AType *A)#

BLAS Level 2 Functions#

BLAS Level 2 routines perform matrix-vector operations, such as matrix-vector multiplication, rank-1 and rank-2 matrix updates, and solution of triangular systems.

template<TensorConcept XType, TensorConcept YType>
void einsums::linear_algebra::axpy(typename XType::data_type alpha, const XType &X, YType *Y)#
template<TensorConcept XType, TensorConcept YType>
void einsums::linear_algebra::axpby(typename XType::data_type alpha, const XType &X, typename XType::data_type beta, YType *Y)#
template<bool TransA, MatrixConcept AType, VectorConcept XType, VectorConcept YType, typename U>
void einsums::linear_algebra::gemv(const U alpha, const AType &A, const XType &z, const U beta, YType *y)#

General matrix-vector multiplication.

This function performs one of the matrix-vector operations

\[ y := alpha*A*z + beta*y\mathrm{,\ or\ }y := alpha*A^{T}*z + beta*y, \]
where alpha and beta are scalars, z and y are vectors and A is an \(m\) by \(n\) matrix.

NEED TO ADD AN EXAMPLE
Template Parameters:

  • TransA – Transpose matrix A? true or false

  • AType – The type of the matrix A

  • XType – The type of the vector z

  • YType – The type of the vector y

  • ARank – The rank of the matrix A

  • XYRank – The rank of the vectors z and y

  • T – The underlying data type

Parameters:

  • alpha – Scaling factor for the product of A and z

  • A – Matrix A

  • z – Vector z

  • beta – Scaling factor for the output vector y

  • y – Output vector y

template<MatrixConcept AType, VectorConcept XYType>
void einsums::linear_algebra::ger(typename AType::data_type alpha, const XYType &X, const XYType &Y, AType *A)#

BLAS Level 3 Functions#

BLAS Level 3 routines perform matrix-matrix operations, such as matrix-matrix multiplication, rank-k update, and solutions of triangular systems.

template<bool TransA, bool TransB, MatrixConcept AType, MatrixConcept BType, MatrixConcept CType, typename U>
void einsums::linear_algebra::gemm(const U alpha, const AType &A, const BType &B, const U beta, CType *C)#

General matrix multiplication.

Takes two rank-2 tensors ( A and B ) performs the multiplication and stores the result in to another rank-2 tensor that is passed in ( C ).

In this equation, TransA is op(A) and TransB is op(B).

\[ C = \alpha \;op(A) \;op(B) + \beta C \]

auto A = einsums::create_random_tensor("A", 3, 3);
auto B = einsums::create_random_tensor("B", 3, 3);
auto C = einsums::create_tensor("C", 3, 3);

einsums::linear_algebra::gemm<false, false>(1.0, A, B, 0.0, &C);
Template Parameters:

  • TransA – Tranpose A? true or false

  • TransB – Tranpose B? true or false

  • T – the underlying data type

Parameters:

  • alpha – Scaling factor for the product of A and B

  • A – First input tensor

  • B – Second input tensor

  • beta – Scaling factor for the output tensor C

  • C – Output tensor

template<bool TransA, bool TransB, MatrixConcept AType, MatrixConcept BType, typename U>
remove_view_t<AType> einsums::linear_algebra::gemm(const U alpha, const AType &A, const BType &B)#

General matrix multiplication.

Returns new tensor.

Takes two rank-2 tensors performs the multiplication and returns the result

auto A = einsums::create_random_tensor("A", 3, 3);
auto B = einsums::create_random_tensor("B", 3, 3);

auto C = einsums::linear_algebra::gemm<false, false>(1.0, A, B);
Template Parameters:

  • TransA – Tranpose A?

  • TransB – Tranpose B?

  • T – the underlying data type

Parameters:

  • alpha – Scaling factor for the product of A and B

  • A – First input tensor

  • B – Second input tensor

Returns:

resulting tensor

LAPACK Functions#

LAPACK routines can be divided into the following groups according to the operations they perform:

  • Routines for solving systems of linear equations, factoring and inverting matrices, and estimating condition numbers.

  • Routines for solving least squares problems, eigenvalue and singular value problems, and Sylvester’s equations.

  • Auxiliary and utility routines used to perform certain subtasks, common low-level computation or related tasks.

LAPACK Linear Equation Computational Functions#

Note

These functions assume Fortran, column-major ordering.

template<MatrixConcept TensorType>
int einsums::linear_algebra::getri(TensorType *A, const std::vector<blas_int> &pivot)#

Computes the inverse of a matrix using the LU factorization computed by getrf.

The routine computes the inverse \(inv(A)\) of a general matrix \(A\). Before calling this routine, call getrf to factorize \(A\).

Template Parameters:

  • TensorType – The type of the tensor

  • T – The underlying data type

  • TensorRank – The rank of the tensor

Parameters:

  • A – The matrix to invert

  • pivot – The pivot vector from getrf

Returns:

int If 0, the execution is successful.

template<MatrixConcept TensorType>
int einsums::linear_algebra::getrf(TensorType *A, std::vector<blas_int> *pivot)#

Computes the LU factorization of a general m-by-n matrix.

The routine computes the LU factorization of a general m-by-n matrix A as

\[ A = P*L*U \]
where P is a permutation matrix, L is lower triangular with unit diagonal elements and U is upper triangular. The routine uses partial pivoting, with row interchanges.

Template Parameters:

TensorType –

Parameters:

  • A –

  • pivot –

Returns:

To be classified#

template<bool ComputeLeftRightEigenvectors = true, MatrixConcept AType, VectorConcept WType>
void einsums::linear_algebra::geev(AType *A, WType *W, AType *lvecs, AType *rvecs)#

Compute the general eigendecomposition of a matrix.

Can only be used to compute both left and right eigen vectors or neither.

template<MatrixConcept AType, MatrixConcept BType>
int einsums::linear_algebra::gesv(AType *A, BType *B)#
template<bool ComputeEigenvectors = true, MatrixConcept AType, VectorConcept WType>
void einsums::linear_algebra::heev(AType *A, WType *W)#
template<MatrixConcept TensorType>
void einsums::linear_algebra::invert(TensorType *A)#

Inverts a matrix.

Utilizes the LAPACK routines getrf and getri to invert a matrix.

Template Parameters:

  • TensorType – The type of the tensor

  • T – The underlying data type

  • TensorRank – The rank of the tensor

Parameters:

A – Matrix to invert. On exit, the inverse of A.

enum class einsums::linear_algebra::Norm : char#

Indicates the type of norm to compute.

Values:

enumerator MaxAbs#

\(val = max(abs(Aij))\), largest absolute value of the matrix A.

enumerator One#

\(val = norm1(A)\), 1-norm of the matrix A (maximum column sum)

enumerator Infinity#

\(val = normI(A)\), infinity norm of the matrix A (maximum row sum)

enumerator Frobenius#

\(val = normF(A)\), Frobenius norm of the matrix A (square root of sum of squares).

template<MatrixConcept AType>
RemoveComplexT<AType> einsums::linear_algebra::norm(Norm norm_type, const AType &a)#

Computes the norm of a matrix.

Returns the value of the one norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real matrix A.

using namespace einsums;

auto A = einsums::create_random_tensor("A", 3, 3);
auto norm = einsums::linear_algebra::norm(einsums::linear_algebra::Norm::One, A);

Note

This function assumes that the matrix is stored in column-major order.

Template Parameters:

  • AType – The type of the matrix

  • ADataType – The underlying data type of the matrix

  • ARank – The rank of the matrix

Parameters:

  • norm_type – where Norm::One denotes the one norm of a matrix (maximum column sum), Norm::Infinity denotes the infinity norm of a matrix (maximum row sum) and Norm::Frobenius denotes the Frobenius norm of a matrix (square root of sum of squares). Note that \( max(abs(A(i,j))) \) is not a consistent matrix norm.

  • a – The matrix to compute the norm of

Returns:

The norm of the matrix

template<MatrixConcept AType>
remove_view_t<AType> einsums::linear_algebra::pow(const AType &a, typename AType::data_type alpha, typename AType::data_type cutoff = std::numeric_limits<typename AType::data_type>::epsilon())#

Computes the matrix power of a to alpha.

Return a new tensor, does not destroy a.

Template Parameters:

AType –

Parameters:

  • a – Matrix to take power of

  • alpha – The power to take

  • cutoff – Values below cutoff are considered zero.

Returns:

std::enable_if_t<std::is_base_of_v<Detail::TRTensorBase<double, 2>, AType>, AType>

template<TensorConcept AType>
void einsums::linear_algebra::sum_square(const AType &a, RemoveComplexT<typename AType::data_type> *scale, RemoveComplexT<typename AType::data_type> *sumsq)#

Computes the square sum of a tensor.

returns the values scale_out and sumsq_out such that

\[ (scale_{out}^{2})*sumsq_{out} = a( 1 )^{2} +...+ a( n )^{2} + (scale_{in}^{2})*sumsq_{in}, \]

Under the hood the LAPACK routine lassq is used.

NEED TO ADD AN EXAMPLE
Template Parameters:

  • AType – The type of the tensor

  • ADataType – The underlying data type of the tensor

  • ARank – The rank of the tensor

Parameters:

  • a – The tensor to compute the sum of squares for

  • scale – scale_in and scale_out for the equation provided

  • sumsq – sumsq_in and sumsq_out for the equation provided

template<MatrixConcept AType>
std::tuple<Tensor<typename AType::data_type, 2>, Tensor<RemoveComplexT<AType>, 1>, Tensor<typename AType::data_type, 2>> einsums::linear_algebra::svd(const AType &_A)#
template<MatrixConcept AType>
Tensor<typename AType::data_type, 2> einsums::linear_algebra::svd_nullspace(const AType &_A)#
template<bool ComputeEigenvectors = true, MatrixConcept AType, VectorConcept WType>
void einsums::linear_algebra::syev(AType *A, WType *W)#

Computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix.

This routines assumes the upper triangle of A is stored. The lower triangle is not referenced.

// Create tensors A and b.
auto A = einsums::create_tensor("A", 3, 3);
auto b = einsums::create_tensor("b", 3);

// Fill A with the symmetric data.
A.vector_data() = einsums::VectorData{1.0, 2.0, 3.0, 2.0, 4.0, 5.0, 3.0, 5.0, 6.0};

// On exit, A is destroyed and replaced with the eigenvectors.
// b is replaced with the eigenvalues in ascending order.
einsums::linear_algebra::syev(&A, &b);
Template Parameters:

  • AType – The type of the tensor A

  • ARank – The rank of the tensor A (required to be 2)

  • WType – The type of the tensor W

  • WRank – The rank of the tensor W (required to be 1)

  • T – The underlying data type (required to be real)

  • ComputeEigenvectors – If true, eigenvalues and eigenvectors are computed. If false, only eigenvalues are computed. Defaults to true.

Parameters:

  • A – On entry, the symmetric matrix A in the leading N-by-N upper triangular part of A. On exit, if eigenvectors are requested, the orthonormal eigenvectors of A. Any data previously stored in A is destroyed.

  • W – On exit, the eigenvalues in ascending order.

template<bool ComputeEigenvectors = true, MatrixConcept AType>
std::tuple<remove_view_t<AType>, BasicTensorLike<AType, typename AType::data_type, 1>> einsums::linear_algebra::syev(const AType &A)#

Computes all eigenvalues and, optionally, eigenvectors of a real symmetric matrix.

This routines assumes the upper triangle of A is stored. The lower triangle is not referenced.

// Create tensors A and b.
auto A = einsums::create_tensor("A", 3, 3);

// Fill A with the symmetric data.
A.vector_data() = einsums::VectorData{1.0, 2.0, 3.0, 2.0, 4.0, 5.0, 3.0, 5.0, 6.0};

// On exit, A is not destroyed. The eigenvectors and eigenvalues are returned in a std::tuple.
auto [evecs, evals ] = einsums::linear_algebra::syev(A);
Template Parameters:

  • AType – The type of the tensor A

  • T – The underlying data type (required to be real)

  • ComputeEigenvectors – If true, eigenvalues and eigenvectors are computed. If false, only eigenvalues are computed. Defaults to true.

Parameters:

A – The symmetric matrix A in the leading N-by-N upper triangular part of A.

Returns:

std::tuple<Tensor<T, 2>, Tensor<T, 1>> The eigenvectors and eigenvalues.