Provide an optimized kernel for the multiplication of a matrix with its own transpose

[Klaus Iglberger]

Description

The primary goal of the Blaze library is to provide maximum performance for all operations. However, it doesn't yet provide a special kernel for the multiplication of a matrix with its own transpose:

blaze::DynamicMatrix<double,rowMajor> A, B, C;
// ... Resizing and initialization
B = A * trans(A);
C = trans(A) * A;

In both cases, the resulting matrix is guaranteed to be a symmetric matrix, even in case the matrix A is not symmetric. Thus it is possible to speed up the computation by approx. a factor of 2. Blaze should provide the according optimized kernels that perform better than the kernel for the general matrix multiplication.

Tasks

• provide an optimized kernel for A * trans(A) for both dense and sparse matrices
• provide an optimized kernel for trans(A) * A for both dense and sparse matrices
• guarantee correctness and robustness for all new kernels
• ensure compatibility with all existing matrix classes
• ensure compatibility with all existing matrix expressions
• guarantee maximum performance for the operations
• add the necessary number of test cases for the functionality

[Klaus Iglberger]

Summary

The feature has been implemented, tested, optimized (including vectorization and parallelization) and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 3.1.

In case the matrix resulting from a matrix multiplication is known to be symmetric or Hermitian, the computation can be optimized by explicitly declaring the multiplication as symmetric or Hermitian via the declsym() or declherm() operations, respectively:

using blaze::DynamicMatrix;

DynamicMatrix<double> M1, M2, M3;

// ... Initialization of the matrices

M3 = declsym( M1 * M2 );  // Declare the result of the matrix multiplication as symmetric

using blaze::DynamicMatrix;

DynamicMatrix< complex<double> > M1, M2, M3

// ... Initialization of the matrices

M3 = declherm( M1 * M2 );  // Declare the result of the matrix multiplication as Hermitian

Declaring the multiplication as symmetric or Hermitian can speed up the computation by up to a factor of 2. Note however that the caller of the declsym() and declherm() operations takes full responsibility for the correctness of the declaration. Falsely declaring a multiplication as symmetric or Hermitian leads to undefined behavior!

declsym()

The declsym() operation can be used to explicitly declare any matrix or matrix expression as symmetric:

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

B = declsym( A );

Any matrix or matrix expression that has been declared as symmetric via declsym() will gain all the benefits of a symmetric matrix, which range from reduced runtime checking to a considerable speed-up in computations:

using blaze::DynamicMatrix;
using blaze::SymmetricMatrix;

DynamicMatrix<double> A, B, C;
SymmetricMatrix< DynamicMatrix<double> > S;
// ... Resizing and initialization

isSymmetric( declsym( A ) );  // Will always return true without runtime effort

S = declsym( A );  // Omit any runtime check for symmetry

C = declsym( A * B );  // Declare the result of the matrix multiplication as symmetric
                       // i.e. perform an optimized matrix multiplication

Warning: The declsym() operation has the semantics of a cast: The caller is completely responsible and the system trusts the given information. Declaring a non-symmetric matrix or matrix expression as symmetric via the declsym() operation leads to undefined behavior (which can be violated invariants or wrong computation results)!

declherm()

The declherm() operation can be used to explicitly declare any matrix or matrix expression as Hermitian:

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

B = declherm( A );

Any matrix or matrix expression that has been declared as Hermitian via declherm() will gain all the benefits of an Hermitian matrix, which range from reduced runtime checking to a considerable speed-up in computations:

using blaze::DynamicMatrix;
using blaze::HermitianMatrix;

DynamicMatrix<double> A, B, C;
HermitianMatrix< DynamicMatrix<double> > S;
// ... Resizing and initialization

isHermitian( declherm( A ) );  // Will always return true without runtime effort

S = declherm( A );  // Omit any runtime check for Hermitian symmetry

C = declherm( A * B );  // Declare the result of the matrix multiplication as Hermitian
                        // i.e. perform an optimized matrix multiplication

Warning: The declherm() operation has the semantics of a cast: The caller is completely responsible and the system trusts the given information. Declaring a non-Hermitian matrix or matrix expression as Hermitian via the declherm() operation leads to undefined behavior (which can be violated invariants or wrong computation results)!

[Klaus Iglberger]

Issue #241 was marked as a duplicate of this issue.