Provide support for Cholesky decomposition

[Klaus Iglberger]

Description

The Cholesky decomposition is the method of choice for the decomposition of positive definite matrices and serves as basis for the inversion of such matrices. 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 Cholesky decomposition
• provide an implementation of the Cholesky 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]

Summary

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

Cholesky (LLH) Decomposition

The Cholesky decomposition of a dense matrix can be computed via the llh() function:

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

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

llh( A, L );  // LLH decomposition of a row-major matrix

assert( A == L * ctrans( L ) );

The function works for both rowMajor and columnMajor matrices and the two matrices A and L can have any storage order.

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

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

blaze::LowerMatrix< blaze::DynamicMatrix<double,blaze::columnMajor> > L;

llh( A, L );  // Cholesky decomposition of A

Note that the LLH 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 llh() 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.