[petsc-users] MatPtAP

[Samuel Lanthaler]

Hi,

I was wondering what the most efficient way to use MatPtAP would be in 
the following situation: I am discretizing a PDE system. The 
discretization yields a matrix A that has a band structure (with k upper 
and lower bands, say). In order to implement the boundary conditions, I 
use a transformation matrix P which is essentially the unit matrix, 
except for the entries P_{ij} where i,j<k and n-i,n-j<k, so

P =  [ B, 0, 0, 0, ..., 0, 0 ]
        [  0, 1, 0, 0, ..., 0, 0 ]
        [                              ]
        [                              ]
        [                  ..., 1, 0 ]
        [  0, 0, 0, 0, ..., 0, C ]

with B,C are (k-by-k) matrices.
Right now, I'm simply constructing A, P and calling

     CALL 
MatPtAP(petsc_matA,petsc_matP,MAT_INITIAL_MATRIX,PETSC_DEFAULT_REAL,petsc_matPtAP,ierr)

where I haven't done anything to pestc_matPtAP, prior to this call. Is 
this the way to do it?

I'm asking because, currently, setting up the matrices A and P takes 
very little time, whereas the operation MatPtAP is taking quite long, 
which seems very odd... The matrices are of type MPIAIJ. In my problem, 
the total matrix dimension is around 10'000 and the matrix blocks (B,C) 
are of size ~100.

Thanks in advance for any ideas.

Cheers,
Samuel