Feature request: Block Sparse matrices

[Lakshman]

Most physics based simulations (continuum mechanics, classical mechanics, etc.) involves block sparse matrices. Though one can use blaze::CompressedMatrix, which has block size 1, higher performance is possible with native support for block sparse matrices. It would be great native block sparse matrices can be included in Blaze.

Due to the lack of libraries with block sparse support many researchers are developing their own codes. Availability of block sparse matrix support in a high performing library like Blaze will surely be a good news for many.

[Klaus Iglberger]

Hi Lakshman!

Thanks a lot for your proposal. We are aware that block structured matrices are very important for many scientific applications. This is why Blaze already supports this feature. Please take a look at the wiki page about block structured vectors and matrices. The examples only use dense matrices and vectors, but it also works as well for sparse vectors and matrices, which enables you to create block structured sparse matrices.

In case you have seen the wiki page before: What functionality is missing from your point of view? If you didn't see it before: Do you think the existing functionality is sufficient for your application?

Best regards,

Klaus!

[Lakshman]

Dear Klaus,

Thanks for the prompt reply.

I am sorry, I did not notice that I can use block structured matrices can be used for sparse systems. I will try small example with block structured sparse matrices and write you if I have any doubts.

A small suggestion. Wouldn't it be better to mention in the documentation about compressed block sparse matrices. Some users may not notice that block structured can be used with sparse matrices, just like me.

Just for information, let me write about the problems I am working with. I am working with two kind of problems. One has sparse matrices with 3x3 block size; quite standard I guess. The second involves mixed block sizes. The diagonals are either 3x3 or 6x6, and hence the off-diagonals can be 3x3, 3x6 or 6x3. Surely I can use 3x3 block sizes for the second problem. However, it is more convenient if I can use the original block sizes.

Thank you very much for the great library.

[Klaus Iglberger]

Hi Lakshman!

This is a reasonable suggestion. We will use this issue for the update of the wiki and the tutorial.

Your second problem with diagonal elements of different block size can be efficiently solved with HybridMatrix:

using BlockSparseMatrix = blaze::CompressedMatrix< blaze::HybridMatrix<double,6UL,6UL> >;

HybridMatrix uses static memory just as StaticMatrix, but can be resized within the specified bounds. Since you have at maximum 6 rows and 6 columns, a 6-by-6 HybridMatrix will do the trick. Alternatively you can of course also use DynamicMatrix, which uses dynamic memory:

using BlockSparseMatrix = blaze::CompressedMatrix< blaze::DynamicMatrix<double> >;

Best regards,

Klaus!

[Lakshman]

Dear Klaus,

Great! Blaze is more versatile than I thought. I will give it a try soon.

Thanks for the tip, and the great work.

Lalith

[Klaus Iglberger]

Summary

We have updated the tutorial and wiki with a second example for block structured vectors and matrices. The second example uses a compressed block structured matrix with blocks of varying size.

The updated tutorial is immediately available via cloning the Blaze repository. The updated wiki can be found here.

Compressed Block Structured Matrices

This example shows the multiplication between a compressed block matrix with blocks of varying size and a compressed block vector:

using namespace blaze;

// ( ( 1 -2  3 )         ( 5 -1 ) )   ( (  1 ) )   ( ( -3 ) )
// ( ( 4  1  0 )         ( 1  2 ) )   ( (  0 ) )   ( (  7 ) )
// ( ( 0  2  4 )         ( 3  1 ) )   ( (  1 ) )   ( (  3 ) )
// (                              )   (        )   (        )
// (              ( 1 )           ) * ( (  2 ) ) = ( (  2 ) )
// (                              )   (        )   (        )
// ( ( 0 -1  1 )         ( 1  0 ) )   ( ( -1 ) )   ( (  0 ) )
// ( ( 2 -1  2 )         ( 0  1 ) )   ( (  2 ) )   ( (  6 ) )

using M3x3 = HybridMatrix<int,3UL,3UL,rowMajor>;
using V3   = HybridVector<int,3UL,columnVector>;

CompressedMatrix<M3x3,rowMajor> A( 3UL, 3UL, 5UL );
A(0,0) = M3x3{ { 1, -2, 3 }, { 4, 1, 0 }, { 0, 2, 4 } };
A(0,2) = M3x3{ { 5, -1 }, { 1, 2 }, { 3, 1 } };
A(1,1) = M3x3{ { 1 } };
A(2,0) = M3x3{ { 0, -1, 1 }, { 2, -1, 2 } };
A(2,2) = M3x3{ { 1, 0 }, { 0, 1 } };

CompressedVector<V3,columnVector> x( 3UL, 3UL );
x[0] = V3{ 1, 0, 1 };
x[1] = V3{ 2 };
x[2] = V3{ -1, 2 };

CompressedVector<V3,columnVector> y( A * x );

[Lakshman]

Dear Klaus,

Thank you very much. It is really helpful.

With best regards Lalith

[Mikhail Katliar]

Dear Klaus,

With your approach to block matrices, inconsistent block sizes are possible. For example, the following code will run:

    using namespace blaze;

    using DynM = DynamicMatrix<double>;
    DynamicMatrix<DynM> A{ { DynM(1, 2), DynM(3, 4) },
                           { DynM(5, 6), DynM(7, 8) } };

A can be seen as a matrix of matrices, but not as a matrix of numbers. The main problem with this is that the assignment to a numeric matrix will not work:

    DynamicMatrix<int> A1;
    A1 = A; // compiler error

I have cases when it is desirable to have a block view of a numeric matrix, and at the same time to be able to use it as a normal numeric matrix.

Therefore I believe a BlockMatrixView<> class would be helpful. The BlockMatrixView<> ctor would take 3 arguments: a reference to a matrix m, and two integer arrays si and sj defining the row and column sizes of the blocks. The operator()(size_t i, size_t j) would return a Submatrix representing the block of m with si[i] rows and sj[j] columns.

Best regards,

Mikhail