mixed matrix type?

[Matthew Knepley]

On Wed, May 21, 2008 at 11:14 AM, Sean Dettrick <sdettrick at gmail.com> wrote:
> Hi,
>
> I have a sparse N*N matrix generated from a DA and a 5 point stencil, with a
> total of approx 5*N non-zero entries.  Now I would like to extend this
> matrix by adding a smaller M*M dense matrix to the bottom right hand side,
> i.e. so that there is a dense square in the bottom right hand corner of the
> otherwise sparse matrix.  The total number of new non-zero entries, M*M, is
> comparable to the total number of old entries, 5*N.  On top of this, there
> would be a small number of non-zero entries in the new upper-right and
> lower-left rectangular portions of the matrix, due to coupling of the two
> systems.  The new total matrix size (including zeroes) would be (N+M)*(N+M).
>
> Can anybody recommend a Mat type to store the new matrix?
>
> One possibility I was thinking of was to establish the original sparse Mat
> with a DA in a sub-communicator (with half the CPUs), and get the ownership
> range with MatGetOwnershipRange.  Then in the petsc_comm_world communicator,
> the complete matrix could be constructed (by element-wise copying the old
> one I suppose), and the ownership range could be maintained manually.
>
> Does this sound like a reasonable strategy?

Not really. If you only want parallelism suitable for the M matrix, I
would use a Schur
complement strategy:

  1) Make the DA matrix

  2) Make M, and the coupling matrices B, B^T

  3) Make a MatShell that has a MatMult() that applies

       B^T A^{-1} B + M

   with A^{-1} done using a KSPSolve().

This seems easier to program and solve to me than the fully assembled one. We
hope later to have something built to do the full assembly.

  Matt

> I would very much appreciate any suggestions or advice.
>
> Thanks,
>
> Sean Dettrick
>
>



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