[petsc-users] Efficient Block Matrix implementation

[tibo at berkeley.edu]

> On Fri, Apr 27, 2012 at 3:50 PM, <tibo at berkeley.edu> wrote:
>
>> 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 ?
>>
>
> What exactly do you want to do with the blocks? Its easy to extract them
> with
> MatGetSubmatrix() and they can be used in the PCFIELDSPLIT preconditioner.
>
>    Matt
>
>

Thank you for your answer (and previous one).
I don't need to use block decomposition, my goal is just to solve as
efficiently as possible systems of the form Ax=b at each time step. So
what I do now as a first approach is that I declare A as SEQAIJ, then
fills it at each time step (since A changes at each time step), and then
solves many times a system of the form Ax = b (Runge Kutta steps with
inner Newton iterations) using KSPsolve and go to the next time step.
Since my matrix is block-tridiagonal, I was just wondering if I could take
advantage of this using a SEQBAIJ declaration, which seems not to be the
case since the blocks are sparse.
The next step will be to implement this sequence in parallel, but I wanted
to make sure that I was using the appropriate sequential implementation
first to base myself upon

>> 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
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>