1. Blaze Avatar Blaze
  2. Blaze
  3. blaze

Wiki

Clone wiki

blaze / LAPACK Functions

View History

Introduction

The Blaze library makes extensive use of the LAPACK functionality for various compute tasks (including the decomposition, inversion and the computation of the determinant of dense matrices). For this purpose, Blaze implements several convenient C++ wrapper functions for all required LAPACK functions. The following sections give a complete overview of all available LAPACK wrapper functions. For more details on the individual LAPACK functions see the Blaze function documentation or the LAPACK online documentation browser.

Most of the wrapper functions are implemented as thin wrappers around LAPACK functions. They provide the parameters of the original LAPACK functions and thus provide maximum flexibility:

using blaze::blas_int_t;

constexpr size_t N( 100UL );

blaze::DynamicMatrix<double,blaze::columnMajor> A( N, N );
// ... Initializing the matrix

const blas_int_t m    ( numeric_cast<blas_int_t>( A.rows()    ) );  // == N
const blas_int_t n    ( numeric_cast<blas_int_t>( A.columns() ) );  // == N
const blas_int_t lda  ( numeric_cast<blas_int_t>( A.spacing() ) );  // >= N
const blas_int_t lwork( n*lda );

const std::unique_ptr<blas_int_t[]> ipiv( new blas_int_t[N] );  // No initialization required
const std::unique_ptr<double[]> work( new double[N] );          // No initialization required

blas_int_t info( 0 );

getrf( m, n, A.data(), lda, ipiv.get(), &info );                  // Reports failure via 'info'
getri( n, A.data(), lda, ipiv.get(), work.get(), lwork, &info );  // Reports failure via 'info'

In this context, blas_int_t is either a 32-bit or 64-bit signed integral type, depending on the setting of the BLAZE_BLAS_IS_64BIT compilation switch (see BLAS Mode).

Additionally, Blaze provides wrappers that provide a higher level of abstraction. These wrappers provide a maximum of convenience:

using blaze::blas_int_t;

constexpr size_t N( 100UL );

blaze::DynamicMatrix<double,blaze::columnMajor> A( N, N );
// ... Initializing the matrix

const std::unique_ptr<blas_int_t[]> ipiv( new blas_int_t[N] );  // No initialization required

getrf( A, ipiv.get() );  // Cannot fail
getri( A, ipiv.get() );  // Reports failure via exception

Note that all functions only work for general, non-adapted matrices with float, double, complex<float>, or complex<double> element type. The attempt to call the function with adaptors or matrices of any other element type results in a compile time error!

Also note that all functions can only be used if a fitting LAPACK library is available and linked to the final executable. Otherwise a call to this function will result in a linker error.

Finally, please note that for performance reasons all functions do only provide the basic exception safety guarantee, i.e. in case an exception is thrown the given matrix may already have been modified.

Matrix Decomposition

The following functions decompose/factorize the given dense matrix. Based on this decomposition the matrix can be inverted or used to solve a linear system of equations.

LU Decomposition

The following functions provide an interface for the LAPACK functions sgetrf(), dgetrf(), cgetrf(), and zgetrf(), which compute the LU decomposition for the given general matrix:

namespace blaze {

void getrf( blas_int_t m, blas_int_t n, float* A, blas_int_t lda, blas_int_t* ipiv, blas_int_t* info );

void getrf( blas_int_t m, blas_int_t n, double* A, blas_int_t lda, blas_int_t* ipiv, blas_int_t* info );

void getrf( blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, blas_int_t* ipiv, blas_int_t* info );

void getrf( blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, blas_int_t* ipiv, blas_int_t* info );

template< typename MT, bool SO >
void getrf( DenseMatrix<MT,SO>& A, blas_int_t* ipiv );

} // namespace blaze

The decomposition has the form

A = P L U,

where P is a permutation matrix, L is a lower triangular matrix, and U is an upper triangular matrix. The resulting decomposition is stored within A: In case of a column-major matrix, L is stored in the lower part of A and U is stored in the upper part. The unit diagonal elements of L are not stored. In case A is a row-major matrix the result is transposed.

Note that the LU decomposition will never fail, even for singular matrices. However, in case of a singular matrix the resulting decomposition cannot be used for a matrix inversion or solving a linear system of equations.

LDLT Decomposition

The following functions provide an interface for the LAPACK functions ssytrf(), dsytrf(), csytrf(), and zsytrf(), which compute the LDLT (Bunch-Kaufman) decomposition for the given symmetric indefinite matrix:

namespace blaze {

void sytrf( char uplo, blas_int_t n, float* A, blas_int_t lda, blas_int_t* ipiv, float* work, blas_int_t lwork, blas_int_t* info );

void sytrf( char uplo, blas_int_t n, double* A, blas_int_t lda, blas_int_t* ipiv, double* work, blas_int_t lwork, blas_int_t* info );

void sytrf( char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, blas_int_t* ipiv, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void sytrf( char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, blas_int_t* ipiv, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void sytrf( DenseMatrix<MT,SO>& A, char uplo, blas_int_t* ipiv );

} // namespace blaze

The decomposition has the form

A = U D U^{T} (if uplo = 'U'), or
A = L D L^{T} (if uplo = 'L'),

where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The resulting decomposition is stored within A: In case uplo is set to 'L' the result is stored in the lower part of the matrix and the upper part remains untouched, in case uplo is set to 'U' the result is stored in the upper part and the lower part remains untouched.

Note that the Bunch-Kaufman decomposition will never fail, even for singular matrices. However, in case of a singular matrix the resulting decomposition cannot be used for a matrix inversion or solving a linear system of equations.

LDLH Decomposition

The following functions provide an interface for the LAPACK functions chetrf() and zsytrf(), which compute the LDLH (Bunch-Kaufman) decomposition for the given Hermitian indefinite matrix:

namespace blaze {

void hetrf( char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, blas_int_t* ipiv, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void hetrf( char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, blas_int_t* ipiv, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void hetrf( DenseMatrix<MT,SO>& A, char uplo, blas_int_t* ipiv );

} // namespace blaze

The decomposition has the form

A = U D U^{H} (if uplo = 'U'), or
A = L D L^{H} (if uplo = 'L'),

where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The resulting decomposition is stored within A: In case uplo is set to 'L' the result is stored in the lower part of the matrix and the upper part remains untouched, in case uplo is set to 'U' the result is stored in the upper part and the lower part remains untouched.

Note that the Bunch-Kaufman decomposition will never fail, even for singular matrices. However, in case of a singular matrix the resulting decomposition cannot be used for a matrix inversion or solving a linear system of equations.

Cholesky Decomposition

The following functions provide an interface for the LAPACK functions spotrf(), dpotrf(), cpotrf(), and zpotrf(), which compute the Cholesky (LLH) decomposition for the given positive definite matrix:

namespace blaze {

void potrf( char uplo, blas_int_t n, float* A, blas_int_t lda, blas_int_t* info );

void potrf( char uplo, blas_int_t n, double* A, blas_int_t lda, blas_int_t* info );

void potrf( char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, blas_int_t* info );

void potrf( char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, blas_int_t* info );

template< typename MT, bool SO >
void potrf( DenseMatrix<MT,SO>& A, char uplo );

} // namespace blaze

The decomposition has the form

A = U^{T} U (if uplo = 'U'), or
A = L L^{T} (if uplo = 'L'),

where U is an upper triangular matrix and L is a lower triangular matrix. The Cholesky decomposition fails if the given matrix A is not a positive definite matrix. In this case a std::std::invalid_argument exception is thrown.

QR Decomposition

The following functions provide an interface for the LAPACK functions sgeqrf(), dgeqrf(), cgeqrf(), and zgeqrf(), which compute the QR decomposition of the given general matrix:

namespace blaze {

void geqrf( blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void geqrf( blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double* tau, double* work, blas_int_t lwork, blas_int_t* info );

void geqrf( blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void geqrf( blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void geqrf( DenseMatrix<MT,SO>& A, typename MT::ElementType* tau );

} // namespace blaze

The decomposition has the form

A = Q R,

where the Q is represented as a product of elementary reflectors

Q = H(1) H(2) . . . H(k), with k = min(m,n).

Each H(i) has the form

H(i) = I - tau * v * v^T,

where tau is a real scalar, and v is a real vector with v(0:i-1) = 0 and v(i) = 1. v(i+1:m) is stored on exit in A(i+1:m,i), and tau in tau(i). Thus on exit the elements on and above the diagonal of the matrix contain the min(m,n)-by-n upper trapezoidal matrix R (R is upper triangular if m >= n); the elements below the diagonal, with the array tau, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors.

The following functions provide an interface for the LAPACK functions sorgqr(), dorgqr(), sorg2r(), dorg2r(), cungqr(), zunqqr(), cung2r(), and zung2r(), which reconstruct the Q matrix from a QR decomposition:

namespace blaze {

void orgqr( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void orgqr( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void orgqr( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void org2r( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t* info );

void org2r( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t* info );

template< typename MT, bool SO >
void org2r( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ungqr( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void ungqr( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void ungqr( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ung2r( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t* info );

void ung2r( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t* info );

template< typename MT, bool SO >
void ung2r( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );

} // namespace blaze

The following functions provide an interface for the LAPACK functions sormqr(), dormqr(), cunmqr(), and zunmqr(), which can be used to multiply a matrix with the Q matrix from a QR decomposition:

namespace blaze {

void ormqr( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const float* A, blas_int_t lda, const float* tau, float* C, blas_int_t ldc, float* work, blas_int_t lwork, blas_int_t* info );

void ormqr( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const double* A, blas_int_t lda, const double* tau, double* C, blas_int_t ldc, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void ormqr( DenseMatrix<MT1,SO1>& C, const DenseMatrix<MT2,SO2>& A, char side, char trans, const ElementType_<MT2>* tau );


void unmqr( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* C, blas_int_t ldc, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void unmqr( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* C, blas_int_t ldc, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO, typename MT2 >
void unmqr( DenseMatrix<MT1,SO>& C, DenseMatrix<MT2,SO>& A, char side, char trans, ElementType_<MT2>* tau );

} // namespace blaze

RQ Decomposition

The following functions provide an interface for the LAPACK functions sgerqf(), dgerqf(), cgerqf(), and zgerqf(), which compute the RQ decomposition of the given general matrix:

namespace blaze {

void gerqf( blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void gerqf( blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double* tau, double* work, blas_int_t lwork, blas_int_t* info );

void gerqf( blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void gerqf( blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void gerqf( DenseMatrix<MT,SO>& A, typename MT::ElementType* tau );

} // namespace blaze

The decomposition has the form

A = R Q,

where the Q is represented as a product of elementary reflectors

Q = H(1) H(2) . . . H(k), with k = min(m,n).

Each H(i) has the form

H(i) = I - tau * v * v^T,

where tau is a real scalar, and v is a real vector with v(n-k+i+1:n) = 0 and v(n-k+i) = 1. v(1:n-k+i-1) is stored on exit in A(m-k+i,1:n-k+i-1), and tau in tau(i). Thus in case m <= n, the upper triangle of the subarray A(1:m,n-m+1:n) contains the m-by-m upper triangular matrix R and in case m >= n, the elements on and above the (m-n)-th subdiagonal contain the m-by-n upper trapezoidal matrix R; the remaining elements in combination with the array tau represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors.

The following functions provide an interface for the LAPACK functions sorgrq(), dorgrq(), sorgr2(), dorgr2(), cungrq(), zunqrq(), cungr2(), and zunqr2(), which reconstruct the Q matrix from a RQ decomposition:

namespace blaze {

void orgrq( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void orgrq( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void orgrq( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void orgr2( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t* info );

void orgr2( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t* info );

template< typename MT, bool SO >
void orgr2( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ungrq( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void ungrq( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void ungrq( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ungr2( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t* info );

void ungr2( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t* info );

template< typename MT, bool SO >
void ungr2( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );

} // namespace blaze

The following functions provide an interface for the LAPACK functions sormrq(), dormrq(), cunmrq(), and zunmrq(), which can be used to multiply a matrix with the Q matrix from a RQ decomposition:

namespace blaze {

void ormrq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const float* A, blas_int_t lda, const float* tau, float* C, blas_int_t ldc, float* work, blas_int_t lwork, blas_int_t* info );

void ormrq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const double* A, blas_int_t lda, const double* tau, double* C, blas_int_t ldc, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void ormrq( DenseMatrix<MT1,SO1>& C, const DenseMatrix<MT2,SO2>& A, char side, char trans, const ElementType_<MT2>* tau );


void unmrq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* C, blas_int_t ldc, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void unmrq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* C, blas_int_t ldc, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO, typename MT2 >
void unmrq( DenseMatrix<MT1,SO>& C, DenseMatrix<MT2,SO>& A, char side, char trans, ElementType_<MT2>* tau );

} // namespace blaze

QL Decomposition

The following functions provide an interface for the LAPACK functions sgeqlf(), dgeqlf(), cgeqlf(), and zgeqlf(), which compute the QL decomposition of the given general matrix:

namespace blaze {

void geqlf( blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void geqlf( blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double* tau, double* work, blas_int_t lwork, blas_int_t* info );

void geqlf( blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void geqlf( blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void geqlf( DenseMatrix<MT,SO>& A, typename MT::ElementType* tau );

} // namespace blaze

The decomposition has the form

A = Q L,

where the Q is represented as a product of elementary reflectors

Q = H(k) . . . H(2) H(1), with k = min(m,n).

Each H(i) has the form

H(i) = I - tau * v * v^T,

where tau is a real scalar, and v is a real vector with v(m-k+i+1:m) = 0 and v(m-k+i) = 1. v(1:m-k+i-1) is stored on exit in A(1:m-k+i-1,n-k+i), and tau in tau(i). Thus in case m >= n, the lower triangle of the subarray A(m-n+1:m,1:n) contains the n-by-n lower triangular matrix L and in case m <= n, the elements on and below the (n-m)-th subdiagonal contain the m-by-n lower trapezoidal matrix L; the remaining elements in combination with the array tau represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors.

The following functions provide an interface for the LAPACK functions sorgql(), dorgql(), sorg2l(), dorg2l(), cungql(), zungql(), cung2l(), and zung2l(), which reconstruct the Q matrix from an QL decomposition:

namespace blaze {

void orgql( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void orgql( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void orgql( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void org2l( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t* info );

void org2l( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t* info );

template< typename MT, bool SO >
void org2l( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ungql( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void ungql( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void ungql( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ung2l( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t* info );

void ung2l( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t* info );

template< typename MT, bool SO >
void ung2l( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );

} // namespace blaze

The following functions provide an interface for the LAPACK functions sormql(), dormql(), cunmql(), and zunmql(), which can be used to multiply a matrix with the Q matrix from a QL decomposition:

namespace blaze {

void ormql( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const float* A, blas_int_t lda, const float* tau, float* C, blas_int_t ldc, float* work, blas_int_t lwork, blas_int_t* info );

void ormql( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const double* A, blas_int_t lda, const double* tau, double* C, blas_int_t ldc, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void ormql( DenseMatrix<MT1,SO1>& C, const DenseMatrix<MT2,SO2>& A, char side, char trans, const ElementType_<MT2>* tau );


void unmql( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* C, blas_int_t ldc, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void unmql( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* C, blas_int_t ldc, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO, typename MT2 >
void unmql( DenseMatrix<MT1,SO>& C, DenseMatrix<MT2,SO>& A, char side, char trans, ElementType_<MT2>* tau );

} // namespace blaze

LQ Decomposition

The following functions provide an interface for the LAPACK functions sgelqf(), dgelqf(), cgelqf(), and zgelqf(), which compute the LQ decomposition of the given general matrix:

namespace blaze {

void gelqf( blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void gelqf( blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double* tau, double* work, blas_int_t lwork, blas_int_t* info );

void gelqf( blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void gelqf( blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void gelqf( DenseMatrix<MT,SO>& A, typename MT::ElementType* tau );

} // namespace blaze

The decomposition has the form

A = L Q,

where the Q is represented as a product of elementary reflectors

Q = H(k) . . . H(2) H(1), with k = min(m,n).

Each H(i) has the form

H(i) = I - tau * v * v^T,

where tau is a real scalar, and v is a real vector with v(0:i-1) = 0 and v(i) = 1. v(i+1:n) is stored on exit in A(i,i+1:n), and tau in tau(i). Thus on exit the elements on and below the diagonal of the matrix contain the m-by-min(m,n) lower trapezoidal matrix L (L is lower triangular if m <= n); the elements above the diagonal, with the array tau, represent the orthogonal matrix Q as a product of min(m,n) elementary reflectors.

The following functions provide an interface for the LAPACK functions sorglq(), dorglq(), sorgl2(), dorgl2(), cunglq(), zunqlq(), cungl2(), and zunql2(), which reconstruct the Q matrix from an LQ decomposition:

namespace blaze {

void orglq( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t lwork, blas_int_t* info );

void orglq( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void orglq( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void orgl2( blas_int_t m, blas_int_t n, blas_int_t k, float* A, blas_int_t lda, const float* tau, float* work, blas_int_t* info );

void orgl2( blas_int_t m, blas_int_t n, blas_int_t k, double* A, blas_int_t lda, const double* tau, double* work, blas_int_t* info );

template< typename MT, bool SO >
void orgl2( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void unglq( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void unglq( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void unglq( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );


void ungl2( blas_int_t m, blas_int_t n, blas_int_t k, complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* work, blas_int_t* info );

void ungl2( blas_int_t m, blas_int_t n, blas_int_t k, complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* work, blas_int_t* info );

template< typename MT, bool SO >
void ungl2( DenseMatrix<MT,SO>& A, const typename MT::ElementType* tau );

} // namespace blaze

The following functions provide an interface for the LAPACK functions sormlq(), dormlq(), cunmlq(), and zunmlq(), which can be used to multiply a matrix with the Q matrix from a LQ decomposition:

namespace blaze {

void ormlq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const float* A, blas_int_t lda, const float* tau, float* C, blas_int_t ldc, float* work, blas_int_t lwork, blas_int_t* info );

void ormlq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const double* A, blas_int_t lda, const double* tau, double* C, blas_int_t ldc, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void ormlq( DenseMatrix<MT1,SO1>& C, const DenseMatrix<MT2,SO2>& A, char side, char trans, const ElementType_<MT2>* tau );


void unmlq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<float>* A, blas_int_t lda, const complex<float>* tau, complex<float>* C, blas_int_t ldc, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void unmlq( char side, char trans, blas_int_t m, blas_int_t n, blas_int_t k, const complex<double>* A, blas_int_t lda, const complex<double>* tau, complex<double>* C, blas_int_t ldc, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT1, bool SO, typename MT2 >
void unmlq( DenseMatrix<MT1,SO>& C, DenseMatrix<MT2,SO>& A, char side, char trans, ElementType_<MT2>* tau );

} // namespace blaze

Matrix Inversion

Given a matrix that has already been decomposed, the following functions can be used to invert the matrix in-place.

LU-based Inversion

The following functions provide an interface for the LAPACK functions sgetri(), dgetri(), cgetri(), and zgetri(), which invert a general matrix that has already been decomposed by an LU Decomposition:

namespace blaze {

void getri( blas_int_t n, float* A, blas_int_t lda, const blas_int_t* ipiv, float* work, blas_int_t lwork, blas_int_t* info );

void getri( blas_int_t n, double* A, blas_int_t lda, const blas_int_t* ipiv, double* work, blas_int_t lwork, blas_int_t* info );

void getri( blas_int_t n, complex<float>* A, blas_int_t lda, const blas_int_t* ipiv, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void getri( blas_int_t n, complex<double>* A, blas_int_t lda, const blas_int_t* ipiv, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO >
void getri( DenseMatrix<MT,SO>& A, const blas_int_t* ipiv );

} // namespace blaze

The functions fail if ...

  • ... the given matrix is not a square matrix;
  • ... the given matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

LDLT-based Inversion

The following functions provide an interface for the LAPACK functions ssytri(), dsytri(), csytri(), and zsytri(), which invert a symmetric indefinite matrix that has already been decomposed by an LDLT Decomposition:

namespace blaze {

void sytri( char uplo, blas_int_t n, float* A, blas_int_t lda, const blas_int_t* ipiv, float* work, blas_int_t* info );

void sytri( char uplo, blas_int_t n, double* A, blas_int_t lda, const blas_int_t* ipiv, double* work, blas_int_t* info );

void sytri( char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, const blas_int_t* ipiv, complex<float>* work, blas_int_t* info );

void sytri( char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, const blas_int_t* ipiv, complex<double>* work, blas_int_t* info );

template< typename MT, bool SO >
void sytri( DenseMatrix<MT,SO>& A, char uplo, const blas_int_t* ipiv );

} // namespace blaze

The functions fail if ...

  • ... the given matrix is not a square matrix;
  • ... the given matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

LDLH-based Inversion

The following functions provide an interface for the LAPACK functions chetri() and zhetri(), which invert an Hermitian indefinite matrix that has already been decomposed by an LDLH Decomposition:

namespace blaze {

void hetri( char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, const blas_int_t* ipiv, complex<float>* work, blas_int_t* info );

void hetri( char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, const blas_int_t* ipiv, complex<double>* work, blas_int_t* info );

template< typename MT, bool SO >
void hetri( DenseMatrix<MT,SO>& A, char uplo, const blas_int_t* ipiv );

} // namespace blaze

The functions fail if ...

  • ... the given matrix is not a square matrix;
  • ... the given matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

Cholesky-based Inversion

The following functions provide an interface for the LAPACK functions spotri(), dpotri(), cpotri(), and zpotri(), which invert a positive definite matrix that has already been decomposed by an Cholesky Decomposition:

namespace blaze {

void potri( char uplo, blas_int_t n, float* A, blas_int_t lda, blas_int_t* info );

void potri( char uplo, blas_int_t n, double* A, blas_int_t lda, blas_int_t* info );

void potri( char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, blas_int_t* info );

void potri( char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, blas_int_t* info );

template< typename MT, bool SO >
void potri( DenseMatrix<MT,SO>& A, char uplo );

} // namespace blaze

The functions fail if ...

  • ... the given matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the given matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

Inversion of Triangular Matrices

The following functions provide an interface for the LAPACK functions strtri(), dtrtri(), ctrtri(), and ztrtri(), which invert the given triangular matrix in-place:

namespace blaze {

void trtri( char uplo, char diag, blas_int_t n, float* A, blas_int_t lda, blas_int_t* info );

void trtri( char uplo, char diag, blas_int_t n, double* A, blas_int_t lda, blas_int_t* info );

void trtri( char uplo, char diag, blas_int_t n, complex<float>* A, blas_int_t lda, blas_int_t* info );

void trtri( char uplo, char diag, blas_int_t n, complex<double>* A, blas_int_t lda, blas_int_t* info );

template< typename MT, bool SO >
void trtri( DenseMatrix<MT,SO>& A, char uplo, char diag );

} // namespace blaze

The functions fail if ...

  • ... the given matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the given diag argument is neither 'U' nor 'N';
  • ... the given matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

Substitution

Given a matrix that has already been decomposed the following functions can be used to perform the forward/backward substitution step to compute the solution to a system of linear equations. Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems:

Single right-hand side:

  • A *x=b if A is column-major
  • A^T*x=b if A is row-major

Multiple right-hand sides:

  • A *X =B if both A and B are column-major
  • A^T*X =B if A is row-major and B is column-major
  • A *X^T=B^T if A is column-major and B is row-major
  • A^T*X^T=B^T if both A and B are row-major

In this context the general system matrix A is a n-by-n matrix that has already been factorized by the according decomposition function, x and b are n-dimensional vectors and X and B are either row-major m-by-n matrices or column-major n-by-m matrices.

LU-based Substitution

The following functions provide an interface for the LAPACK functions sgetrs(), dgetrs(), cgetrs(), and zgetrs(), which perform the substitution step for a general matrix that has already been decomposed by an LU Decomposition:

namespace blaze {

void getrs( char trans, blas_int_t n, blas_int_t nrhs, const float* A, blas_int_t lda, const blas_int_t* ipiv, float* B, blas_int_t ldb, blas_int_t* info );

void getrs( char trans, blas_int_t n, blas_int_t nrhs, const double* A, blas_int_t lda, const blas_int_t* ipiv, double* B, blas_int_t ldb, blas_int_t* info );

void getrs( char trans, blas_int_t n, const complex<float>* A, blas_int_t lda, const blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void getrs( char trans, blas_int_t n, const complex<double>* A, blas_int_t lda, const blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void getrs( const DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char trans, const blas_int_t* ipiv );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void getrs( const DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char trans, const blas_int_t* ipiv );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Substitution). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations. The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given trans argument is neither 'N' nor 'T' nor 'C';
  • ... the sizes of the two given matrices do not match.

The first four functions report failure via the info argument, the last two functions throw a std::invalid_argument exception in case of an error.

LDLT-based Substitution

The following functions provide an interface for the LAPACK functions ssytrs(), dsytrs(), csytrs(), and zsytrs(), which perform the substitution step for a symmetric indefinite matrix that has already been decomposed by an LDLT Decomposition:

namespace blaze {

void sytrs( char uplo, blas_int_t n, blas_int_t nrhs, const float* A, blas_int_t lda, const blas_int_t* ipiv, float* B, blas_int_t ldb, blas_int_t* info );

void sytrs( char uplo, blas_int_t n, blas_int_t nrhs, const double* A, blas_int_t lda, const blas_int_t* ipiv, double* B, blas_int_t ldb, blas_int_t* info );

void sytrs( char uplo, blas_int_t n, blas_int_t nrhs, const complex<float>* A, blas_int_t lda, const blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void sytrs( char uplo, blas_int_t n, blas_int_t nrhs, const complex<double>* A, blas_int_t lda, const blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void sytrs( const DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo, const blas_int_t* ipiv );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void sytrs( const DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo, const blas_int_t* ipiv );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Substitution). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations. The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the sizes of the two given matrices do not match.

The first four functions report failure via the info argument, the last two functions throw a std::invalid_argument exception in case of an error.

LDLH-based Substitution

The following functions provide an interface for the LAPACK functions chetrs(), and zhetrs(), which perform the substitution step for an Hermitian indefinite matrix that has already been decomposed by an LDLH Decomposition:

namespace blaze {

void hetrs( char uplo, blas_int_t n, blas_int_t nrhs, const complex<float>* A, blas_int_t lda, const blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void hetrs( char uplo, blas_int_t n, blas_int_t nrhs, const complex<double>* A, blas_int_t lda, const blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void hetrs( const DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo, const blas_int_t* ipiv );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void hetrs( const DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo, const blas_int_t* ipiv );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Substitution). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations. The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the sizes of the two given matrices do not match.

The first two functions report failure via the info argument, the last two functions throw a std::invalid_argument exception in case of an error.

Cholesky-based Substitution

The following functions provide an interface for the LAPACK functions spotrs(), dpotrs(), cpotrs(), and zpotrs(), which perform the substitution step for a positive definite matrix that has already been decomposed by an Cholesky Decomposition:

namespace blaze {

void potrs( char uplo, blas_int_t n, blas_int_t nrhs, const float* A, blas_int_t lda, float* B, blas_int_t ldb, blas_int_t* info );

void potrs( char uplo, blas_int_t n, blas_int_t nrhs, const double* A, blas_int_t lda, double* B, blas_int_t ldb, blas_int_t* info );

void potrs( char uplo, blas_int_t n, blas_int_t nrhs, const complex<float>* A, blas_int_t lda, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void potrs( char uplo, blas_int_t n, blas_int_t nrhs, const complex<double>* A, blas_int_t lda, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void potrs( const DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void potrs( const DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Substitution). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations. The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the sizes of the two given matrices do not match.

The first two functions report failure via the info argument, the last two functions throw a std::invalid_argument exception in case of an error.

Substitution for Triangular Matrices

The following functions provide an interface for the LAPACK functions strtrs(), dtrtrs(), ctrtrs(), and ztrtrs(), which perform the substitution step for a triangular matrix:

namespace blaze {

void trtrs( char uplo, char trans, char diag, blas_int_t n, blas_int_t nrhs, const float* A, blas_int_t lda, float* B, blas_int_t ldb, blas_int_t* info );

void trtrs( char uplo, char trans, char diag, blas_int_t n, blas_int_t nrhs, const double* A, blas_int_t lda, double* B, blas_int_t ldb, blas_int_t* info );

void trtrs( char uplo, char trans, char diag, blas_int_t n, blas_int_t nrhs, const complex<float>* A, blas_int_t lda, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void trtrs( char uplo, char trans, char diag, blas_int_t n, blas_int_t nrhs, const complex<double>* A, blas_int_t lda, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void trtrs( const DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo, char trans, char diag );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void trtrs( const DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo, char trans, char diag );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Substitution). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations. The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the given trans argument is neither 'N' nor 'T' nor 'C';
  • ... the given diag argument is neither 'U' nor 'N';
  • ... the sizes of the two given matrices do not match.

The first four functions report failure via the info argument, the last two functions throw a std::invalid_argument exception in case of an error.

Linear System Solver

The following functions represent compound functions that perform both the decomposition step as well as the substitution step to compute the solution to a system of linear equations. Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems:

Single right-hand side:

  • A *x=b if A is column-major
  • A^T*x=b if A is row-major

Multiple right-hand sides:

  • A *X =B if both A and B are column-major
  • A^T*X =B if A is row-major and B is column-major
  • A *X^T=B^T if A is column-major and B is row-major
  • A^T*X^T=B^T if both A and B are row-major

In this context the general system matrix A is a n-by-n matrix that has already been factorized by the according decomposition function, x and b are n-dimensional vectors and X and B are either row-major m-by-n matrices or column-major n-by-m matrices.

LU-based Linear System Solver

The following functions provide an interface for the LAPACK functions sgesv(), dgesv(), cgesv(), and zgesv(), which combine an LU Decomposition and the according LU-based Substitution:

namespace blaze {

void gesv( blas_int_t n, blas_int_t nrhs, float* A, blas_int_t lda, blas_int_t* ipiv, float* B, blas_int_t ldb, blas_int_t* info );

void gesv( blas_int_t n, blas_int_t nrhs, double* A, blas_int_t lda, blas_int_t* ipiv, double* B, blas_int_t ldb, blas_int_t* info );

void gesv( blas_int_t n, blas_int_t nrhs, complex<float>* A, blas_int_t lda, blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void gesv( blas_int_t n, blas_int_t nrhs, complex<double>* A, blas_int_t lda, blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void gesv( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, blas_int_t* ipiv );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void gesv( DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, blas_int_t* ipiv );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Linear System Solver). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations and A has been decomposed by means of an LU Decomposition.

The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given system matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

LDLT-based Linear System Solver

The following functions provide an interface for the LAPACK functions ssysv(), dsysv(), csysv(), and zsysv(), which combine an LDLT Decomposition and the according LDLT-based Substitution:

namespace blaze {

void sysv( char uplo, blas_int_t n, blas_int_t nrhs, float* A, blas_int_t lda, blas_int_t* ipiv, float* B, blas_int_t ldb, float* work, blas_int_t lwork, blas_int_t* info );

void sysv( char uplo, blas_int_t n, blas_int_t nrhs, double* A, blas_int_t lda, blas_int_t* ipiv, double* B, blas_int_t ldb, double* work, blas_int_t lwork, blas_int_t* info );

void sysv( char uplo, blas_int_t n, blas_int_t nrhs, complex<float>* A, blas_int_t lda, blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void sysv( char uplo, blas_int_t n, blas_int_t nrhs, complex<double>* A, blas_int_t lda, blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void sysv( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo, blas_int_t* ipiv );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void sysv( DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo, blas_int_t* ipiv );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Linear System Solver). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations and A has been decomposed by means of an LDLT Decomposition.

The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the sizes of the two given matrices do not match;
  • ... the given system matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

LDLH-based Linear System Solver

The following functions provide an interface for the LAPACK functions shesv(), dhesv(), chesv(), and zhesv(), which combine an LDLH Decomposition and the according LDLH-based Substitution:

namespace blaze {

void hesv( char uplo, blas_int_t n, blas_int_t nrhs, complex<float>* A, blas_int_t lda, blas_int_t* ipiv, complex<float>* B, blas_int_t ldb, complex<float>* work, blas_int_t lwork, blas_int_t* info );

void hesv( char uplo, blas_int_t n, blas_int_t nrhs, complex<double>* A, blas_int_t lda, blas_int_t* ipiv, complex<double>* B, blas_int_t ldb, complex<double>* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void hesv( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo, blas_int_t* ipiv );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void hesv( DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo, blas_int_t* ipiv );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Linear System Solver). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations and A has been decomposed by means of an LDLH Decomposition.

The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the sizes of the two given matrices do not match;
  • ... the given system matrix is singular and not invertible.

The first two functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

Cholesky-based Linear System Solver

The following functions provide an interface for the LAPACK functions sposv(), dposv(), cposv(), and zposv(), which combine an Cholesky Decomposition and the according Cholesky-based Substitution:

namespace blaze {

void posv( char uplo, blas_int_t n, blas_int_t nrhs, float* A, blas_int_t lda, float* B, blas_int_t ldb, blas_int_t* info );

void posv( char uplo, blas_int_t n, blas_int_t nrhs, double* A, blas_int_t lda, double* B, blas_int_t ldb, blas_int_t* info );

void posv( char uplo, blas_int_t n, blas_int_t nrhs, complex<float>* A, blas_int_t lda, complex<float>* B, blas_int_t ldb, blas_int_t* info );

void posv( char uplo, blas_int_t n, blas_int_t nrhs, complex<double>* A, blas_int_t lda, complex<double>* B, blas_int_t ldb, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void posv( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo );

template< typename MT1, bool SO1, typename MT2, bool SO2 >
void posv( DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& B, char uplo );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Linear System Solver). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations and A has been decomposed by means of an Cholesky Decomposition.

The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the sizes of the two given matrices do not match;
  • ... the given system matrix is singular and not invertible.

The first four functions report failure via the info argument, the fifth function throws a std::invalid_argument exception in case of an error.

Linear System Solver for Triangular Matrices

The following functions provide an interface for the LAPACK functions strsv(), dtrsv(), ctrsv(), and ztrsv():

namespace blaze {

void trsv( char uplo, char trans, char diag, blas_int_t n, const float* A, blas_int_t lda, float* x, blas_int_t incX );

void trsv( char uplo, char trans, char diag, blas_int_t n, const double* A, blas_int_t lda, double* x, blas_int_t incX );

void trsv( char uplo, char trans, char diag, blas_int_t n, const complex<float>* A, blas_int_t lda, complex<float>* x, blas_int_t incX );

void trsv( char uplo, char trans, char diag, blas_int_t n, const complex<double>* A, blas_int_t lda, complex<double>* x, blas_int_t incX );

template< typename MT, bool SO, typename VT, bool TF >
void trsv( const DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& b, char uplo, char trans, char diag );

} // namespace blaze

Note that depending on the storage order of the system matrix and the given right-hand side the functions solve different equation systems (see Linear System Solver). If the function exits successfully, the vector b or the matrix B contain the solution(s) of the linear system of equations.

The functions fail if ...

  • ... the given system matrix is not a square matrix;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the given trans argument is neither 'N' nor 'T' nor 'C';
  • ... the given diag argument is neither 'U' nor 'N'.

The last function throws a std::invalid_argument exception in case of an error. Note that none of the functions does perform any test for singularity or near-singularity. Such tests must be performed prior to calling this function!

Eigenvalues/Eigenvectors

General Matrices

The following functions provide an interface for the LAPACK functions sgeev(), dgeev(), cgeev(), and zgeev(), which compute the eigenvalues and optionally the eigenvectors of the given general matrix:

namespace blaze {

void geev( char jobvl, char jobvr, blas_int_t n, float* A, blas_int_t lda, float* wr, float* wi, float* VL, blas_int_t ldvl, float* VR, blas_int_t ldvr, float* work, blas_int_t lwork, blas_int_t* info );

void geev( char jobvl, char jobvr, blas_int_t n, double* A, blas_int_t lda, double* wr, double* wi, double* VL, blas_int_t ldvl, double* VR, blas_int_t ldvr, double* work, blas_int_t lwork, blas_int_t* info );

void geev( char jobvl, char jobvr, blas_int_t n, complex<float>* A, blas_int_t lda, complex<float>* w, complex<float>* VL, blas_int_t ldvl, complex<float>* VR, blas_int_t ldvr, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* info );

void geev( char jobvl, char jobvr, blas_int_t n, complex<double>* A, blas_int_t lda, complex<double>* w, complex<double>* VL, blas_int_t ldvl, complex<double>* VR, blas_int_t ldvr, complex<double>* work, blas_int_t lwork, double* rwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void geev( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w );

template< typename MT1, bool SO1, typename MT2, bool SO2, typename VT, bool TF >
void geev( DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& VL, DenseVector<VT,TF>& w );

template< typename MT1, bool SO1, typename VT, bool TF, typename MT2, bool SO2 >
void geev( DenseMatrix<MT1,SO1>& A, DenseVector<VT,TF>& w, DenseMatrix<MT2,SO2>& VR );

template< typename MT1, bool SO1, typename MT2, bool SO2, typename VT, bool TF, typename MT3, bool SO3 >
void geev( DenseMatrix<MT1,SO1>& A, DenseMatrix<MT2,SO2>& VL, DenseVector<VT,TF>& w, DenseMatrix<MT3,SO3>& VR );

} // namespace blaze

The complex eigenvalues of the given matrix A are returned in the given vector w. Please note that no order of eigenvalues can be assumed, except that complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first.

If VR is provided as an argument, the right eigenvectors are returned in the rows of VR in case VR is a row-major matrix and in the columns of VR in case VR is a column-major matrix. The right eigenvector v[j] of A satisfies

A * v[j] = lambda[j] * v[j],

where lambda[j] is its eigenvalue.

If VL is provided as an argument, the left eigenvectors are returned in the rows of VL in case VL is a row-major matrix and in the columns of VL in case VL is a column-major matrix. The left eigenvector u[j] of A satisfies

u[j]^{H} * A = lambda[j] * u[j]^{H},

where u[j]^{H} denotes the conjugate transpose of u[j].

w, VL, and VR are resized to the correct dimensions (if possible and necessary). The functions fail if ...

  • ... the given matrix A is not a square matrix;
  • ... the given matrix VL is a fixed size matrix and the dimensions don't match;
  • ... the given vector w is a fixed size vector and the size doesn't match;
  • ... the given matrix VR is a fixed size matrix and the dimensions don't match;
  • ... the eigenvalue computation fails.

The first four functions report failure via the info argument, the last four functions throw an exception in case of an error.

Symmetric Matrices

The following functions provide an interface for the LAPACK functions ssyev() and dsyev(), which compute the eigenvalues and eigenvectors of the given symmetric matrix:

namespace blaze {

void syev( char jobz, char uplo, blas_int_t n, float* A, blas_int_t lda, float* w, float* work, blas_int_t lwork, blas_int_t* info );

void syev( char jobz, char uplo, blas_int_t n, double* A, blas_int_t lda, double* w, double* work, blas_int_t lwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void syev( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char jobz, char uplo );

} // namespace blaze

Alternatively, the following functions can be used, which provide an interface to the LAPACK functions ssyevd() and dsyevd(). In contrast to the syev() functions they use a divide-and-conquer strategy for the computation of the left and right eigenvectors:

namespace blaze {

void syevd( char jobz, char uplo, blas_int_t n, float* A, blas_int_t lda, float* w, float* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t liwork, blas_int_t* info );

void syevd( char jobz, char uplo, blas_int_t n, double* A, blas_int_t lda, double* w, double* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t liwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void syevd( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char jobz, char uplo );

} // namespace blaze

The real eigenvalues are returned in ascending order in the given vector w. w is resized to the correct size (if possible and necessary). In case A is a row-major matrix, the left eigenvectors are returned in the rows of A, in case A is a column-major matrix, the right eigenvectors are returned in the columns of A.

The functions fail if ...

  • ... the given matrix A is not a square matrix;
  • ... the given vector w is a fixed size vector and the size doesn't match;
  • ... the given jobz argument is neither 'V' nor 'N';
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the eigenvalue computation fails.

The first two functions report failure via the info argument, the last function throws an exception in case of an error.

Via the following functions, which wrap the LAPACK functions ssyevx() and dsyevx(), it is possible to compute a subset of eigenvalues and/or eigenvectors of a symmetric matrix:

namespace blaze {

void syevx( char jobz, char range, char uplo, blas_int_t n, float* A, blas_int_t lda, float vl, float vu, blas_int_t il, blas_int_t iu, float abstol, blas_int_t* m, float* w, float* Z, blas_int_t ldz, float* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t* ifail, blas_int_t* info );

void syevx( char jobz, char range, char uplo, blas_int_t n, double* A, blas_int_t lda, double vl, double vu, blas_int_t il, blas_int_t iu, double abstol, blas_int_t* m, double* w, double* Z, blas_int_t ldz, double* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t* ifail, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
size_t syevx( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char uplo );

template< typename MT, bool SO, typename VT, bool TF, typename ST >
size_t syevx( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char uplo, ST low, ST upp );

template< typename MT1, bool SO1, typename VT, bool TF, typename MT2, bool SO2 >
size_t syevx( DenseMatrix<MT1,SO1>& A, DenseVector<VT,TF>& w, DenseMatrix<MT2,SO2>& Z, char uplo );

template< typename MT1, bool SO1, typename VT, bool TF, typename MT2, bool SO2, typename ST >
size_t syevx( DenseMatrix<MT1,SO1>& A, DenseVector<VT,TF>& w, DenseMatrix<MT2,SO2>& Z, char uplo, ST low, ST upp );

} // namespace blaze

The number of eigenvalues to be computed is specified by the lower bound low and the upper bound upp, which either form an integral or a floating point range.

In case low and upp are of integral type, the function computes all eigenvalues in the index range [low..upp]. The num resulting real eigenvalues are stored in ascending order in the given vector w, which is either resized (if possible) or expected to be a num-dimensional vector. The eigenvectors are returned in the rows of Z in case Z is row-major matrix and in the columns of Z in case Z is a column-major matrix. Z is resized (if possible) or expected to be a num-by-n row-major matrix or a n-by-num column-major matrix.

In case low and upp are of floating point type, the function computes all eigenvalues in the half-open interval (low..upp]. The resulting real eigenvalues are stored in ascending order in the given vector w, which is either resized (if possible) or expected to be an n-dimensional vector. The eigenvectors are returned in the rows of Z in case Z is a row-major matrix and in the columns of Z in case Z is a column-major matrix. Z is resized (if possible) or expected to be a n-by-n matrix.

The functions fail if ...

  • ... the given matrix A is not a square matrix;
  • ... the given vector w is a fixed size vector and the size doesn't match;
  • ... the given matrix Z is a fixed size matrix and the dimensions don't match;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the eigenvalue computation fails.

The first two functions report failure via the info argument, the last four functions throw an exception in case of an error.

Hermitian Matrices

The following functions provide an interface for the LAPACK functions cheev() and zheev(), which compute the eigenvalues and eigenvectors of the given Hermitian matrix:

namespace blaze {

void heev( char jobz, char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, float* w, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* info );

void heev( char jobz, char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, double* w, complex<double>* work, blas_int_t lwork, float* rwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void heev( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char jobz, char uplo );

} // namespace blaze

Alternatively, the following functions can be used, which provide an interface to the LAPACK functions cheevd() and zheevd(). In contrast to the heev() functions they use a divide-and-conquer strategy for the computation of the left and right eigenvectors:

namespace blaze {

void heevd( char jobz, char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, float* w, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* lrwork, blas_int_t* iwork, blas_int_t* liwork, blas_int_t* info );

void heevd( char jobz, char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, double* w, complex<double>* work, blas_int_t lwork, double* rwork, blas_int_t lrwork, blas_int_t* iwork, blas_int_t* liwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void heevd( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char jobz, char uplo );

} // namespace blaze

The real eigenvalues are returned in ascending order in the given vector w. w is resized to the correct size (if possible and necessary). In case A is a row-major matrix, the left eigenvectors are returned in the rows of A, in case A is a column-major matrix, the right eigenvectors are returned in the columns of A.

The functions fail if ...

  • ... the given matrix A is not a square matrix;
  • ... the given vector w is a fixed size vector and the size doesn't match;
  • ... the given jobz argument is neither 'V' nor 'N';
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the eigenvalue computation fails.

The first two functions report failure via the info argument, the last function throws an exception in case of an error.

Via the following functions, which wrap the LAPACK functions cheevx() and zheevx(), it is possible to compute a subset of eigenvalues and/or eigenvectors of an Hermitian matrix:

namespace blaze {

void heevx( char jobz, char range, char uplo, blas_int_t n, complex<float>* A, blas_int_t lda, float vl, float vu, blas_int_t il, blas_int_t iu, float abstol, blas_int_t* m, float* w, complex<float>* Z, blas_int_t ldz, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* iwork, blas_int_t* ifail, blas_int_t* info );

void heevx( char jobz, char range, char uplo, blas_int_t n, complex<double>* A, blas_int_t lda, double vl, double vu, blas_int_t il, blas_int_t iu, double abstol, blas_int_t* m, double* w, complex<double>* Z, blas_int_t ldz, complex<double>* work, blas_int_t lwork, double* rwork, blas_int_t* iwork, blas_int_t* ifail, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
size_t heevx( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char uplo );

template< typename MT, bool SO, typename VT, bool TF, typename ST >
size_t heevx( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& w, char uplo, ST low, ST upp );

template< typename MT1, bool SO1, typename VT, bool TF, typename MT2, bool SO2 >
size_t heevx( DenseMatrix<MT1,SO1>& A, DenseVector<VT,TF>& w, DenseMatrix<MT2,SO2>& Z, char uplo );

template< typename MT1, bool SO1, typename VT, bool TF, typename MT2, bool SO2, typename ST >
size_t heevx( DenseMatrix<MT1,SO1>& A, DenseVector<VT,TF>& w, DenseMatrix<MT2,SO2>& Z, char uplo, ST low, ST upp );

} // namespace blaze

The number of eigenvalues to be computed is specified by the lower bound low and the upper bound upp, which either form an integral or a floating point range.

In case low and upp are of integral type, the function computes all eigenvalues in the index range [low..upp]. The num resulting real eigenvalues are stored in ascending order in the given vector w, which is either resized (if possible) or expected to be a num-dimensional vector. The eigenvectors are returned in the rows of Z in case Z is row-major matrix and in the columns of Z in case Z is a column-major matrix. Z is resized (if possible) or expected to be a num-by-n row-major matrix or a n-by-num column-major matrix.

In case low and upp are of floating point type, the function computes all eigenvalues in the half-open interval (low..upp]. The resulting real eigenvalues are stored in ascending order in the given vector w, which is either resized (if possible) or expected to be an n-dimensional vector. The eigenvectors are returned in the rows of Z in case Z is a row-major matrix and in the columns of Z in case Z is a column-major matrix. Z is resized (if possible) or expected to be a n-by-n matrix.

The functions fail if ...

  • ... the given matrix A is not a square matrix;
  • ... the given vector w is a fixed size vector and the size doesn't match;
  • ... the given matrix Z is a fixed size matrix and the dimensions don't match;
  • ... the given uplo argument is neither 'L' nor 'U';
  • ... the eigenvalue computation fails.

The first two functions report failure via the info argument, the last four functions throw an exception in case of an error.

Singular Values/Singular Vectors

The following functions provide an interface for the LAPACK functions sgesvd(), dgesvd(), cgesvd(), and zgesvd(), which perform a singular value decomposition (SVD) on the given general matrix:

namespace blaze {

void gesvd( char jobu, char jobv, blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float* s, float* U, blas_int_t ldu, float* V, blas_int_t ldv, float* work, blas_int_t lwork, blas_int_t* info );

void gesvd( char jobu, char jobv, blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double* s, double* U, blas_int_t ldu, double* V, blas_int_t ldv, double* work, blas_int_t lwork, blas_int_t* info );

void gesvd( char jobu, char jobv, blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, float* s, complex<float>* U, blas_int_t ldu, complex<float>* V, blas_int_t ldv, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* info );

void gesvd( char jobu, char jobv, blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, double* s, complex<double>* U, blas_int_t ldu, complex<double>* V, blas_int_t ldv, complex<double>* work, blas_int_t lwork, double* rwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void gesvd( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& s, char jobu, char jobv );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF >
void gesvd( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, char jobu, char jobv );

template< typename MT1, bool SO, typename VT, bool TF, typename MT2 >
void gesvd( DenseMatrix<MT1,SO>& A, DenseVector<VT,TF>& s, DenseMatrix<MT2,SO>& V, char jobu, char jobv );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF, typename MT3 >
void gesvd( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, DenseMatrix<MT3,SO>& V, char jobu, char jobv );

} // namespace blaze

Alternatively, the following functions can be used, which provide an interface to the LAPACK functions sgesdd(), dgesdd(), cgesdd(), and zgesdd(). In contrast to the gesvd() functions they compute the singular value decomposition (SVD) of the given general matrix by applying a divide-and-conquer strategy for the computation of the left and right singular vectors:

namespace blaze {

void gesdd( char jobz, blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float* s, float* U, blas_int_t ldu, float* V, blas_int_t ldv, float* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t* info );

void gesdd( char jobz, blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double* s, double* U, blas_int_t ldu, double* V, blas_int_t ldv, double* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t* info );

void gesdd( char jobz, blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, float* s, complex<float>* U, blas_int_t ldu, complex<float>* V, blas_int_t ldv, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* iwork, blas_int_t* info );

void gesdd( char jobz, blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, double* s, complex<double>* U, blas_int_t ldu, complex<double>* V, blas_int_t ldv, complex<double>* work, blas_int_t lwork, double* rwork, blas_int_t* iwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
void gesdd( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& s );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF >
void gesdd( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, char jobz );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF >
void gesdd( DenseMatrix<MT1,SO>& A, DenseVector<VT,TF>& s, DenseMatrix<MT2,SO>& V, char jobz );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF, typename MT3 >
void gesdd( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, DenseMatrix<MT3,SO>& V, char jobz );

} // namespace blaze

The resulting decomposition has the form

A = U S V,

where S is a m-by-n matrix, which is zero except for its min(m,n) diagonal elements, U is an m-by-m orthogonal matrix, and V is a n-by-n orthogonal matrix. The diagonal elements of S are the singular values of A, the first min(m,n) columns of U and rows of V are the left and right singular vectors of A, respectively.

The resulting min(m,n) real and non-negative singular values are returned in descending order in the vector s, which is resized to the correct size (if possible and necessary).

Via the following functions, which wrap the LAPACK functions sgesvdx(), dgesvdx(), cgesvdx(), and zgesvdx(), it is possible to compute a subset of singular values and/or vectors:

namespace blaze {

void gesvdx( char jobu, char jobv, char range, blas_int_t m, blas_int_t n, float* A, blas_int_t lda, float vl, float vu, blas_int_t il, blas_int_t iu, blas_int_t* ns, float* s, float* U, blas_int_t ldu, float* V, blas_int_t ldv, float* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t* info );

void gesvdx( char jobu, char jobv, char range, blas_int_t m, blas_int_t n, double* A, blas_int_t lda, double vl, double vu, blas_int_t il, blas_int_t iu, blas_int_t* ns, double* s, double* U, blas_int_t ldu, double* V, blas_int_t ldv, double* work, blas_int_t lwork, blas_int_t* iwork, blas_int_t* info );

void gesvdx( char jobu, char jobv, char range, blas_int_t m, blas_int_t n, complex<float>* A, blas_int_t lda, float vl, float vu, blas_int_t il, blas_int_t iu, blas_int_t* ns, float* s, complex<float>* U, blas_int_t ldu, complex<float>* V, blas_int_t ldv, complex<float>* work, blas_int_t lwork, float* rwork, blas_int_t* iwork, blas_int_t* info );

void gesvdx( char jobu, char jobv, char range, blas_int_t m, blas_int_t n, complex<double>* A, blas_int_t lda, double vl, double vu, blas_int_t il, blas_int_t iu, blas_int_t* ns, double* s, complex<double>* U, blas_int_t ldu, complex<double>* V, blas_int_t ldv, complex<double>* work, blas_int_t lwork, double* rwork, blas_int_t* iwork, blas_int_t* info );

template< typename MT, bool SO, typename VT, bool TF >
size_t gesvdx( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& s );

template< typename MT, bool SO, typename VT, bool TF, typename ST >
size_t gesvdx( DenseMatrix<MT,SO>& A, DenseVector<VT,TF>& s, ST low, ST upp );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF >
size_t gesvdx( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF, typename ST >
size_t gesvdx( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, ST low, ST upp );

template< typename MT1, bool SO, typename VT, bool TF, typename MT2 >
size_t gesvdx( DenseMatrix<MT1,SO>& A, DenseVector<VT,TF>& s, DenseMatrix<MT2,SO>& V );

template< typename MT1, bool SO, typename VT, bool TF, typename MT2, typename ST >
size_t gesvdx( DenseMatrix<MT1,SO>& A, DenseVector<VT,TF>& s, DenseMatrix<MT2,SO>& V, ST low, ST upp );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF, typename MT3 >
size_t gesvdx( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, DenseMatrix<MT3,SO>& V );

template< typename MT1, bool SO, typename MT2, typename VT, bool TF, typename MT3, typename ST >
size_t gesvdx( DenseMatrix<MT1,SO>& A, DenseMatrix<MT2,SO>& U, DenseVector<VT,TF>& s, DenseMatrix<MT3,SO>& V, ST low, ST upp );

} // namespace blaze

The number of singular values to be computed is specified by the lower bound low and the upper bound upp, which either form an integral or a floating point range.

In case low and upp form are of integral type, the function computes all singular values in the index range [low..upp]. The num resulting real and non-negative singular values are stored in descending order in the given vector s, which is either resized (if possible) or expected to be a num-dimensional vector. The resulting left singular vectors are stored in the given matrix U, which is either resized (if possible) or expected to be a m-by-num matrix. The resulting right singular vectors are stored in the given matrix V, which is either resized (if possible) or expected to be a num-by-n matrix.

In case low and upp are of floating point type, the function computes all singular values in the half-open interval (low..upp]. The resulting real and non-negative singular values are stored in descending order in the given vector s, which is either resized (if possible) or expected to be a min(m,n)-dimensional vector. The resulting left singular vectors are stored in the given matrix U, which is either resized (if possible) or expected to be a m-by-min(m,n) matrix. The resulting right singular vectors are stored in the given matrix V, which is either resized (if possible) or expected to be a min(m,n)-by-n matrix.

The functions fail if ...

  • ... the given matrix U is a fixed size matrix and the dimensions don't match;
  • ... the given vector s is a fixed size vector and the size doesn't match;
  • ... the given matrix V is a fixed size matrix and the dimensions don't match;
  • ... the given scalar values don't form a proper range;
  • ... the singular value decomposition fails.

The first four functions report failure via the info argument, the remaining functions throw an exception in case of an error.

Previous: BLAS Functions ---- Next: Block Vectors and Matrices

Updated