> 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 >
