[petsc-users] Efficient Block Matrix implementation

[tibo at berkeley.edu]

Hi,

I am using petsc in a Fortran code to solve large systems of linear
equations Ax = b where A is a matrix resulting from a discretization of a
reaction diffusion problem: A is block-tridiagonal, and each of the three
blocks is sparse. It turns out that I know exactly the sparsity structure
of each block (ie the number and position of non-zeros per row).

For now, I am only writing a sequential code. From the petsc manual, it
seems there are two ways to create my matrix A: Either using
MatCreateSeqAIJ, or using MatCreateSeqBAIJ. The second approach seems
recommended since A is indeed a Block matrix.

However, the command MatCreateSeqAIJ allows me to specify the number of
non zeros per row, therefore allowing me to take advantage of knowing the
sparsity structure of my blocks (but losing the advantage of using the
fact that the matrix is a block matrix). It seems that the command
MatCreateSeqBAIJ only allows me to specify the number of non zeros blocks
per block rows, therefore losing the available information on the sparsity
of each block.

My question is then: Is there a way to declare the matrix to take
advantage of both the fact that my matrix is a block matrix (or even
better, that it is a block tridiagonal matrix) and that I know the
sparsity structure of the blocks ?

I am asking this question because If I use the block matrix approach as
described above (just specifying the number of non zeros blocks per block
rows) versus the AIJ approach with specifying the number of non zeros for
each row, the computing time increases by a factor of 5, every other
parameters being unchanged.

Thank you