Provide support for LU decomposition

[Klaus Iglberger]

Description

The LU decomposition is an essential step for solving systems of linear equations. Unfortunately, this feature is still missing in Blaze. As it is an often required algorithm, Blaze should provide support for it.

Tasks

• design the interface for the LU decomposition
• provide an implementation of the LU decomposition that ...
  • ... is robust
  • ... is as efficient as possible
  • ... can deal with errors in a graceful way
• provide a shared-memory parallel implementation (if possible)
• provide a full documentation of the feature, including ...
  • ... an explanation of the algorithm
  • ... limitations of the algorithms
  • ... examples of use
  • ... example computations
• ensure compatibility to all existing matrix types
• add appropriate test cases for the feature

[Klaus Iglberger]

The Blaze library already provides LAPACK wrapper functions for the (P)LU decomposition of a dense matrix:

namespace blaze {

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

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

} // namespace blaze

These functions immediately available via cloning the Blaze repository and will be officially released in Blaze 2.6.

[Klaus Iglberger]

Summary

The feature has been implemented, tested, optimized for different kinds of matrices, and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 2.6.

LU Decomposition

The LU decomposition of a dense matrix can be computed via the lu() function:

blaze::DynamicMatrix<double,blaze::rowMajor> A;
// ... Resizing and initialization

blaze::DynamicMatrix<double,blaze::rowMajor> L, U, P;

lu( A, L, U, P );  // LU decomposition of a row-major matrix

assert( A == L * U * P );

blaze::DynamicMatrix<double,blaze::columnMajor> A;
// ... Resizing and initialization

blaze::DynamicMatrix<double,blaze::columnMajor> L, U, P;

lu( A, L, U, P );  // LU decomposition of a column-major matrix

assert( A == P * L * U );

The function works for both rowMajor and columnMajor matrices. Note, however, that the three matrices A, L and U are required to have the same storage order. Also, please note that the way the permutation matrix P needs to be applied differs between row-major and column-major matrices, since the algorithm uses column interchanges for row-major matrices and row interchanges for column-major matrices.

Furthermore, lu() can be used with adaptors. For instance, the following example demonstrates the LU decomposition of a symmetric matrix into a lower and upper triangular matrix:

blaze::SymmetricMatrix< blaze::DynamicMatrix<double,blaze::columnMajor> > A;
// ... Resizing and initialization

blaze::LowerMatrix< blaze::DynamicMatrix<double,blaze::columnMajor> > L;
blaze::UpperMatrix< blaze::DynamicMatrix<double,blaze::columnMajor> > U;
blaze::DynamicMatrix<double,blaze::columnMajor> P;

lu( A, L, U, P );  // LU decomposition of A

Note that the LU decomposition can only be used for dense matrices with float, double, complex<float> or complex<double> element type. The attempt to call the function with matrices of any other element type or with a sparse matrix results in a compile time error!

Also note that the lu() function decomposes a dense matrix by means of LAPACK kernels. Thus the function can only be used if the fitting LAPACK library is available and linked to the executable. Otherwise a linker error will be created.