Provide lower and upper views for dense and sparse matrices

[Klaus Iglberger]

Description

Since version 2.3 the Blaze library provides the LowerMatrix and UpperMatrix adaptors, which have the compile time guarantee to be square triangular matrices. Although these adaptors in general are very useful and beneficial for performance, they are not the right tool for all situations. A very valuable addition to the tool set of the Blaze library would be lower and upper triangular views.

Blaze should provide an API which enables to "cut out" the lower and/or upper part of any type of matrix. The views should act as regular dense or sparse matrix, respectively, and thus it should be possible to assign to them, assign them to any other matrix, and it should be possible to use them in arithmetic operations. During these operations they should provide the maximum possible performance.

Conceptual Example

// Setup of matrix A
//
//        ( 0 0 0 )
//        ( 0 0 0 )
//    A = ( 0 0 0 )
//        ( 0 0 0 )
//        ( 0 0 0 )
//
blaze::DynamicMatrix<int,blaze::rowMajor> A( 5UL, 3UL, 0 );

// Setup of matrix B
//
//        ( 9 9 9 9 9 )
//    B = ( 9 9 9 9 9 )
//        ( 9 9 9 9 9 )
//
blaze::DynamicMatrix<int,blaze::rowMajor> B( 3UL, 5UL, 9 );

// Setting the lower part of A to 2
//
//        ( 2 0 0 )
//        ( 2 2 0 )
//    A = ( 2 2 2 )
//        ( 2 2 2 )
//        ( 2 2 2 )
//
lower( A ) = 2;

// Setting the upper part of B to 5
//
//        ( 5 5 5 5 5 )
//    B = ( 9 5 5 5 5 )
//        ( 9 9 5 5 5 )
//
upper( B ) = 5;

// Subtracting the transpose of the lower part of matrix A from the upper part of B
//
//        ( 3 3 3 3 3 )
//    B = ( 9 3 3 3 3 )
//        ( 9 9 3 3 3 )
//
upper( B ) -= trans( lower( A ) );

Tasks

• design an API for lower and upper triangular views
• implement the lower triangular view for dense and sparse matrices
• implement the upper triangular view for dense and sparse matrices
• provide a full documentation of the feature
• ensure compatibility with all existing matrix classes
• ensure compatibility with all existing matrix expressions/operations
• guarantee maximum performance during use of views
• add the necessary number of test cases for the entire functionality

[Klaus Iglberger]

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

[Mikhail Katliar]

Hello Klaus,

What you propose would be a very useful feature. In my case, I need to multiply a matrix with a lower-triangular result of Cholesky decomposition, calculated by potrf():

DynamicMatrix<double> L(m, m);
DynamicMatrix<double> A(m, n);
// Init L, A ...

// Compute Cholesky decomposition of L. Only the lower part of L is referenced.
potrf(L, 'L');

// Multiply a lower-triangular matrix with a full matrix (similar to BLAS dtrmm())
DynamicMatrix<double> LA = lower(L) * A;

As long as the lower() function is not implemented, how would you suggest to code the matrix multiplication in line 9 for maximum performance?

Best regards,

Mikhail

[Klaus Iglberger]

Hi Mikhail!

If you would make sure the matrix L is indeed a lower triangular matrix (i.e. all elements in the upper part are reset to 0) then you can speed up the multiplication by declaring the L matrix as lower triangular via the decllow() function:

DynamicMatrix<double> LA = decllow(L) * A;

From my experience this is the only way how you can effectively gain performance from exploiting a lower matrix during a matrix multiplication. Although the lower() function sound like a great idea and is very convenient to use, it would incur way too many conditionals during the computation and therefore would probably result in decreased performance. The lower() function would therefore best serve to prepare a lower matrix. Note, however, that due to vectorization the optimization via the decllow() function does not pay off for tiny and small matrices. I hope this helps,

Best regards,

Klaus!

[Mikhail Katliar]

Thank you Klaus!

[GK Palem]

Came across this while looking for a way to assign a general matrix (symmetric) to a StrictlyUpperMatrix; The non-upper portion of this general matrix is not zero (being symmetric). This makes it difficult to extract the upper triangular portion with assignment.

DynamicMatrix<bool, blaze::rowMajor> M(5, 5, 0);
StrictlyUpperMatrix<DynamicMatrix<bool, blaze::rowMajor> > U;

M(3, 4) = 1;
M(4, 3) = 1; // symmetric

U = M; // assign only the upper portion of the matrix
// throws error "Invalid setup of strictly upper matrix"

1. Is there a way to make the assignment as “extract and assign only the upper portion of the source matrix” ignoring the lower portion?
2. Alternately, is there a way to “force the non-upper triangular” portion of the general matrix M to be zero (without having to iterate over each element) so that it can be used as the source of assignment ?

The only way I see to use a `StrictlyUpperMatrix` is either you already have a strictly upper matrix or create one from scratch. You can NOT extract one from a general matrix. Is this correct? Would really a love a way to:

1. “extract” the triangular portion of the general matrix without having to iterate over.
2. “override” the non-upper (non-lower/diagonal) portions of a general matrix with default values (say 0) to “convert it in-place” into an upper / lower etc.

‌

[Klaus Iglberger]

Hi GK!

What you would need is exactly what is proposed in this issue. However, this feature is not implemented yet, hence currently there is no better (good?) way to deal with only the lower or upper part. We suggest to implement the according functionality yourself in some form of functions: assignLower(), assignUpper(), etc.

Best regards,

Klaus!