PETSc and block matrices

[Uday K]

Dear all, I have been trying to use deal.ii with PETSc to solve my own set of FE classes. In other words, I am presently only using the linear algebra functionality of the library and have been facing some trouble. 

I would like to define a block matrix, say of 2x2 size, whose individual blocks are themselves sparse matrices. Similarly, I would like to compose a 2x1 block vector for the solution and RHS vector corresponding to a matrix equation. To start with, I would like to solve the resulting system using a direct solver such as MatrixMUMPS. This is what I have done so far: I create the four individual sparse matrices using the PETScWrappers::SparseMatrix class, where I also specify the number of non-zero entries per row. My own function then goes in and sets each of the non-zero entries. 

1. Now, how do I efficiently setup the block matrix? Ideally, I would like to use a functionality where I can define the number of (block) rows and columns and then place pointers to the existing sparse matrices. I set out by trying to use the PETScWrappers::BlockSparseMatrix class, but it isn't clear from the documentation if I can achieve what I want. For starters, I tried to declare an empty object and used the reinit function to initialize a 2x2 block matrix, only to get a compilation error.

error: request for member ‘reinit’ in ‘A_TM’, which is of non-class type ‘dealii::PETScWrappers::BlockSparseMatrix()’

2. A simple code snippet which demonstrates how to implement these features, as provided in the BlockMatrixArray documentation would be great (http://www.dealii.org/8.1.0/doxygen/deal.II/classBlockMatrixArray.html). Of course, BlockMatrixArray can't be used to compose objects of the class PETScWrappers::SparseMatrix.

 Thanks and regards,

Uday 

[Toby Young]

> 1. Now, how do I efficiently setup the block matrix? Ideally, I would like to use a functionality where I can define the number of (block) rows and columns and then place pointers to the existing sparse matrices. I set out by trying to use the PETScWrappers::BlockSparseMatrix class, but it isn't clear from the documentation if I can achieve what I want.

It does not quite work like this.PETScWrappers::BlockSparseMatrix *is* a PETScWrappers::SparseMatrix.

The block renumbering makes it look like a block matrix is some internal wizardry (if I remember correctly this is covered by operators in a BaseClass functionality).

I do not believe it is possible to *efficiently* do what you want to do (but I like your idea!),  The only method I know of is to copy your four individual matrices into a BlockSparseMatrix element-wise. That, unfortunately is very not efficient! Additionally, you have no solvers for the problem you want set up..

I recommend you assemble a block matrix directly (see section "Block solvers and preconditioners" in the examples).

Can that be done for your case? Is there a specific reason you want to do this?

Best,

   Toby

[Uday Khankhoje]

On Thu, Mar 6, 2014 at 7:50 PM, Toby Young <oneli...@gmail.com> wrote:

I recommend you assemble a block matrix directly (see section "Block solvers and preconditioners" in the examples).

Can that be done for your case? Is there a specific reason you want to do this?

​Thanks for the clarifications, Toby. 

It should be possible for me to directly assemble a block matrix, though it isn't the most intuitive way of thinking about the problem. My motivation ​for sectioning the matrix into blocks comes from the fact that I am trying to solve a complex valued problem by formulating its equivalent in real space (i.e. combining two real equations into one matrix equation). So, the block elements turn out very naturally to be the real and imaginary components of the original complex matrix equation. Further, I would like to run some optimization routines using these matrices, and am led to computing gradients etc., and once again, it is intuitive to formulate in terms of block matrices and vectors.

Before I attempt to assemble the block matrix directly, I was wondering if I could 

(1) define/assemble each sparse matrix using dealii::SparseMatrix, and then

(2) define a block matrix using dealii::BlockMatrixArray (where I can use its member function 'enter ' to assign the blocks to the matrices assembled in (1) above, and finally

(3) use SparseDirectUMFPACK to solve the resulting equation? 

I know that UMFPACK is designed for unsymmetric matrices.. is its performance undermined if it is applied to anti-symmetric matrices? I ask because ideally I would like to use MatrixMUMPS, but I had issues with MUMPS in my installation, detailed below;

Per Matthias's suggestion, I had to compile deal.ii with MUMPS support switched off. Here is the reason according to him: 

"The problem is simply that the mumps configuration module does not pick up scotch at the moment. I'll fix that later.

For the moment, I suggest that you just reconfigure deal.II with

disabled, direct mumps support: -DDEAL_II_WITH_MUMPS=OFF. Please note,

that you can use mumps (if you wish to do so) over the appropriate petsc

wrapper as well."

(See https://groups.google.com/d/msg/dealii/12kW9ifjENo/FbI9pNrmXucJ)

So now, I am in a bit of a pickle if UMFPACK won't work for anti-symmetric matrices AND I can't use MUMPS via PETSc!

Thanks and regards, 

 Uday

[Timo Heister]

> It does not quite work like this.PETScWrappers::BlockSparseMatrix *is* a
> PETScWrappers::SparseMatrix.

No, that is not true.

PETScWrappers::BlockSparseMatrix (and
PETScWrappers::MPI::BlockSparseMatrix) are based on BlockMatrixBase<>.
This means that BlockSparseMatrix is a table (basically a 2d array) of
PETScWrappers::SparseMatrix. So each block is a SparseMatrix, which is
a thin wrapper around a PETSc MAT object.
As a consequence, you can not apply any PETSc preconditioner/solver to
a BlockSparseMatrix, but only to individual blocks.

That means for you:
./ You have to assemble everything into a single
PETScWrappers::SparseMatrix to use a solver or preconditioner from
PETSc like MUMPS.
./ The only kind of solvers working on block matrices are deal.II
solvers, but you can of course build a preconditioner in deal.II that
uses PETSc preconditioners for individual blocks.
./ One thing I sometimes do is have a setting in my code to switch
from say a 2x2 BlockSparseMatrix to a 1x1 BlockSparseMatrix (the
construction of the matrix is the only place that needs to change).
Then I can use a PETSc direct solver on matrix.block(0,0) for testing.
./ PETSc has support for block systems, that we are currently not
using. A long-term goal is to implement a way to create a single
PETScWrappers::SparseMatrix and tell PETSc "hey, this is actually a
block matrix, these are the different blocks", so that you can use
PETSc's block solver features. I might have a student working on that
in the summer.

And yes, you should create a BlockSparseMatrix and fill the individual
matrices instead of creating them separately and trying to copying
them in.


--
Timo Heister
http://www.math.clemson.edu/~heister/

[Timo Heister]

> Before I attempt to assemble the block matrix directly, I was wondering if I
> could
> (1) define/assemble each sparse matrix using dealii::SparseMatrix, and then
> (2) define a block matrix using dealii::BlockMatrixArray (where I can use
> its member function 'enter ' to assign the blocks to the matrices assembled
> in (1) above, and finally
> (3) use SparseDirectUMFPACK to solve the resulting equation?

Why don't you create the BlockSparseMatrix and them assemble into
matrix.block(i,j) (which is of type SparseMatrix)?

SparseDirectUMFPACK only works with SparseMatrix and BlockSparseMatrix
(and SparseMatrixEZ).

> I know that UMFPACK is designed for unsymmetric matrices.. is its
> performance undermined if it is applied to anti-symmetric matrices?

Of course it will work. You can throw any non-singular matrix at
UMFPACK. The only limit is that you will run out of memory on large 3d
problems.

> I ask
> because ideally I would like to use MatrixMUMPS, but I had issues with MUMPS
> in my installation, detailed below;

MUMPS is a big pain to install, unfortunately. The way I deal with
this is that I let PETSc install MUMPS and use it through PETSc.

[Uday Khankhoje]

On Thu, Mar 6, 2014 at 9:08 PM, Timo Heister <hei...@clemson.edu> wrote:

Why don't you create the BlockSparseMatrix and them assemble into
matrix.block(i,j) (which is of type SparseMatrix)?

​Thanks! I think this is easy enough to do and is a better solution than placing everything inside a PETScWrappers::SparseMatrix object.

For now my problem is 2D, so the direct solvers are fine, but in the long term it will become essential to have block matrix support through PETSc, particularly for the preconditioners and solvers to work. I hope you get that summer student!

 Best,
Uday

[Uday Khankhoje]

On Thu, Mar 6, 2014 at 9:08 PM, Timo Heister <hei...@clemson.edu> wrote:

Why don't you create the BlockSparseMatrix and them assemble into
matrix.block(i,j) (which is of type SparseMatrix)?

I had a follow up question. Are there any situations where a BlockMatrixArray might be more efficient? For instance, if the constitutive blocks are not distinct matrices, but rather, are scaled/transposed versions of a single matrix. In such a case, the memory requirements of a BlockMatrixArray would be lesser than that of a BlockSparseMatrix, right? ​However, one has lost the functionality of preconditioners/solvers, I guess.

 Thanks,
Uday

[Timo Heister]

> I had a follow up question. Are there any situations where a
> BlockMatrixArray might be more efficient? For instance, if the constitutive
> blocks are not distinct matrices, but rather, are scaled/transposed versions
> of a single matrix. In such a case, the memory requirements of a
> BlockMatrixArray would be lesser than that of a BlockSparseMatrix, right?
> However, one has lost the functionality of preconditioners/solvers, I guess.

Correct. As a consequence it is not a class that is used very often.

[Uday K]

On Thursday, 6 March 2014 21:08:19 UTC+5:30, Timo Heister wrote:

Why don't you create the BlockSparseMatrix and them assemble into
matrix.block(i,j) (which is of type SparseMatrix)?

​Following up on this suggestion, I have a potential bug to report. I cooked up a sparsity pattern and assembled a sparse matrix based on it like so:

unsigned int i,j,k,Ns=3, Nn=2;

//Ns is the number of rows/columns and Nn is the number of non-zero entries per row

std::vector< std::vector<unsigned int> > indices(Ns,std::vector<unsigned int>(Nn));

for(i=0; i<Ns; ++i){

//square matrix, so store the diagonal entry first

indices[i][0] = i;

//Store the rest of the coln indices in increasing order

for(j=1; j<Nn; ++j){

//some placeholder here to ensure sanity, won't work if Nn>2

k = (i+1)%Ns;

//k is the actual coln number that is non-zero

indices[i][j] = k;

}

}

//Copy this structure into the sparsity pattern

SparsityPattern sp;

sp.copy_from(Ns,Ns,indices.begin(),indices.end());

//Assemble the sparse matrix

SparseMatrix<double> mat;

//Using reinit sets all the sparse matrix entries to zero so we can use add later

mat.reinit(sp);

for(i=0; i<Ns; ++i){

for(j=0; j<Nn; j++){

mat.add(i,indices[i][j],i+j);

}

}

I verify that the sparsity pattern and assembled matrix is what I expect. Then, I build a 2x2 BlockSparsityPattern based on the existing SparsityPattern like so;

BlockSparsityPattern bsp(2,2);

for(i=0; i<2; i++){

for(j=0; j<2; j++){

bsp.block(i,j).copy_from(sp);

}

}

I verify that the number of block rows and columns is actually 2.  Then I attempt to assemble a BlockSparseMatrix like so;

BlockSparseMatrix<double> bmat;

bmat.reinit(bsp);

Here is the problem. If I look at the number of blocks in the matrix, it is actually 0, as revealed by examining bmat.n_block_rows()
or bmat.n_block_cols(). Which is why if I try to copy in one of the blocks by bmat.block(1,1).copy_from(mat), I get a core dump. Shouldn't bmat report the correct number of block rows and cols once it has been initialized with bsp? Or am I doing something incorrect here?

 Thanks,
Uday

[Timo Heister]

> BlockSparsityPattern bsp(2,2);
> for(i=0; i<2; i++){
> for(j=0; j<2; j++){
> bsp.block(i,j).copy_from(sp);
> }
> }

you need to call collect_sizes() after changing the size of any block.

[Toby Young]

> Thanks for the clarifications, Toby.

Actually, it seems I had misled you (and myself). Sorry.
Thanks for the clarifications Timo.

> My motivation
> for sectioning the matrix into blocks comes from the fact that I am trying
> to solve a complex valued problem by formulating its equivalent in real
> space (i.e. combining two real equations into one matrix equation).

Aha! The complex problem. Double up on the real-space size and solve
as a 2x2 real matrix problem.
This can easily be done without BlockSparseMatrix. Right now, I can
not remember how I do this.
If you are interested I can send you a code sample when I get back to school.

Best,
Toby

[Uday K]

On Saturday, March 8, 2014 12:44:23 AM UTC+5:30, Toby Young wrote:

Aha! The complex problem. Double up on the real-space size and solve
as a 2x2 real matrix problem.
This can easily be done without BlockSparseMatrix. Right now, I can
not remember how I do this.
If you are interested I can send you a code sample when I get back to school.

Toby, that would be great! I think posting the sample code will benefit a lot of other people on the mailing list as well. The BlockSparseMatrix solution would end up taking four times the memory requirement for the LHS system matrix, whereas a more efficient solution should take about two times the memory (as compared to a single real valued matrix equation).

[Uday K]

you need to call collect_sizes() after changing the size of any block.

Thanks, Timo, that did it. I'm still getting used to deal.ii and sometimes get lost in the sea of documentation!

Meanwhile, after calling collect_sizes(), I can set the matrix entries of each block using the 'add' function, but if I try to copy an entire matrix into the block (i.e. from the code above, 'mat' into any of the blocks of 'bmat' using 'copy_from'), I get a core dump due to an exception being raised on account of different sparsity patterns. I note that that sparsity patterns are two identical, but distinct objects. Since copy_from requires the two matrices to point to the same sparsity object, I guess this is the cause of the exception?

[Timo Heister]

> Since copy_from requires
> the two matrices to point to the same sparsity object, I guess this is the
> cause of the exception?

Yes, we require the SparsityPattern to be identical (not only
containing the same elements). Is there a reason you are creating the
same pattern twice?

If you really need to do this you can do:
dst = 0;
SparseMatrix<double>::iterator it = src.begin(), endit = src.end();
for (; it!=endit; ++it)
dst.add(it->row(), it->column(), it->value());

Note that this won't be very efficient, though.

[Uday Khankhoje]

On Mon, Mar 10, 2014 at 3:16 AM, Timo Heister <hei...@clemson.edu> wrote:

Is there a reason you are creating the
same pattern twice?

​Simply because in my case the constitutive matrices of the BlockMatrix have the same sparsity pattern. Ideally of course, it would be great to be able to place four pointers to the same sparsity pattern, but I don't suppose the class offers that functionality.​

[Timo Heister]

does the following for you?

SparsityPattern sp;
[initialize sp here ...]
BlockSparseMatrix mat(2,1);
mat.block(0,0).reinit(sp);
mat.block(1,0).reinit(sp);
mat.collect_sizes();

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Uday Khankhoje]

On Mon, Mar 10, 2014 at 6:30 PM, Timo Heister <hei...@clemson.edu> wrote:

does the following for you?

SparsityPattern sp;
[initialize sp here ...]
BlockSparseMatrix mat(2,1);
mat.block(0,0).reinit(sp);
mat.block(1,0).reinit(sp);
mat.collect_sizes();

​Nope. It cribs about not finding a constructor of the form (int, int) for the BlockSparseMatrix class (neither could I find one in the documentation of the latest stable version, 8.1.0). 

The reason I have been trying to make these block sparse matrices to work is that I want to solve the complex matrix equation Ax = b by recasting it as A' x' = b', where A'_{00} = A'_{11} = A_r and A'_{10} = -A'{01} = A_i and x', b' are column vectors of the form b'_0 = b_r, b'_1 = b_i, where 'r' and 'i' denote the real and imaginary parts, respectively. This is why all the components of A' have the same sparsity pattern. If you can point me in the direction of an efficient method of implementing this in deal.ii, or can think of something altogether different (as Toby seemed to have suggested), that would be really very helpful.

 Thanks,
Uday

​

[Timo Heister]

> Nope. It cribs about not finding a constructor of the form (int, int) for
> the BlockSparseMatrix class (neither could I find one in the documentation
> of the latest stable version, 8.1.0).

My bad, of course that won't work. Maybe this will though:

BlockSparsityPattern bsp;
[initialize bsp here ...]
BlockSparseMatrix mat(bsp);
mat.block(0,0).reinit(bsp.block(0,0));
mat.block(1,0).reinit(bsp.block(0,0));
// now block0,0 and block1,0 have the identical SP.

> If you can point me in the
> direction of an efficient method of implementing this in deal.ii, or can
> think of something altogether different (as Toby seemed to have suggested),
> that would be really very helpful.

I think you are on the right track and you probably shouldn't worry
about the loss of efficiency when storing a SP twice (premature
optimization is bad).

[Uday Khankhoje]

Great, that did it. Thanks!

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---

You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/tMZLz9u7Qac/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.