[petsc-dev] MATNEST MATSUBMATRIX MATLOCALREF MATBLOCKMAT and friends

[Jed Brown]

On Thu, Jan 13, 2011 at 19:01, Barry Smith <bsmith at mcs.anl.gov> wrote:

> Since Mike McCourt is in dire need of flexible multiphysics infrastructure
> I've started to look at some of the newly implement Mat stuff in PETSc (and
> my head is spinning).
>
>   Could the author's of these things take a look at my comments and respond
> (hopefully with changes, but I can change things also).
>

I am responsible for many of these, but I don't know how much time I'll have
to make changes in the next two weeks.


>
>   Stream of consciousness writing because I can't get my head around it
> all.
>
> ----
>
>   MatCreateNest() (and VecCreateNexs()).  Does not have the the proper
> construction process that always goes through MatCreate_Nest() (and
> VecCreate_Nest()). This needs to be fixed. Should add something like
> MatNestAddBlocks() or similar name for putting the blocks in after the call
> to MatCreate_Nest().
>

We used MatCreateNest() instead of lots of setters because in all cases I
could think of, the user has all the nested Mats available up-front and
adding them one at a time took more effort for the user.  Furthermore, the
internal data structure would need to be more dynamic and we wanted to make
something correct and functional as quickly as possible.  The construction
process could be made more dynamic, but it would make the API a bit bigger,
make the implementation more complex, and I don't think anyone would want to
use that interface anyway.  If you have a use case where it's important to
construct it dynamically, I'll prioritize changing it.


>  MatCreateNest() seems to assume that all Mat's live on all processes in
> the outer Mat communicator. How hard would it be to change the code so that
> each block could live on some subcommunicator?
>

I think it's similar difficulty to making subproblems in PCFieldSplit live
on subcommunicators.  It would also provide by far the most benefit if
PCFieldSplit could support that, so both changes should probably be made at
a similar time.


>  MatCreateSubMatrix() will it eventually have MatSetValues() that directs
> requests down to the appropriate block?
>

It could, I think it's not too hard to add, but I'm not sure how useful it
would be since MatGetLocalSubMatrix() offers, in my opinion, a better API
for this purpose.


>
>
> ----
>
>   MatCreateSubMatrix(). Again not the proper construction process built on
> MatCreate_SubMatrix() (when the heck was this written 10 years ago before we
> enforced that rule.)
>

Not nearly that old, this was added to simplify PCFieldSplit and make more
things work with MFFD and other matrix-free operators.  I had looked at your
MatCreateSchurComplement() just before writing this.

This needs to be fixed. What's with the manual page description "Creates a
> composite matrix that acts as a submatrix" what the heck does 'composite'
> mean here and what does "acts as a submatrix" mean?
>

A_sub = R^T A R, this just holds A and R.


>  MatSubMatrixUpdate() what does the documentation "full matrix in the
> submatrix" mean? The full matrix is  bigger than the submatrix so how can it
> be in the submatrix. Terrible phrase.
>

"original matrix"?  Please suggest better wording.


>
>   MatSubMatrixUpdate() why have this routine, why not model on
> MatGetSubMatrix() and have a single interface function
> MatCreateSubMatrix(Mat,IS,IS,MatRuse,Mat *)?
>
>   MatCreateSubMatrix() terrible name, why not MatGetSubMatrixImplicit()?
>

I don't think the user should ever call these, that's why it is "developer"
level.  You should call MatGetSubMatrix in essentially all cases.  This is
used as the fall-back if the matrix format doesn't have a better way to
produce a submatrix.  Only in some strange case where the user does NOT want
to use the "better way" would they call this function.


> ------
>
>  MatCreateLocalRef() again with no proper construction with
> MatCreate_LocalRef() (come on guys get with the program)
>

Again, users should only need MatGetLocalSubMatrix().


>  MatCreateLocalRef() 'Gets a logical reference to a local submatrix, for
> use in assembly' what does this mean?
>

It allows for multi-physics matrix-assembly assembly using "per-physics"
numbering.  Look at snes/examples/tutorials/ex28.c.  Otherwise the user
would have to translate to global numbering themselves, which sucks and then
wouldn't work with MatNest.  MatLocalRef is pretty much just translating
indices.


> Assembly in PETSc is reserved for inside MatAssembly,VecAssembly it should
> not be used for the process of setting values into the matrix. Thus this
> sentence is meaningless gibberish.  Could possible be something like 'Gets a
> logical reference to a local submatrix, for use in putting values into that
> block with MatSetValues() or MatSetValuesLocal()?
>

Yes, my use of "assembly" referred to the larger process, not
MatAssemblyBegin/End.  Your wording sounds better.  Note that
MatSetValuesBlocked{,Local} also works if the *submatrix* has block
structure (regardless of whether the global matrix does).


>
> MatCreateLocalRef()  How/why is this different than MatCreateSubMatrix()?
> Could they be merged into one construct? Why is this a 'local' ref, is it a
> block only only on this process?
>

Perhaps they could be merged, but the semantics are pretty different.
 MatCreateLocalRef takes "local index sets" and is not collective (and
returns a matrix on COMM_SELF) while MatCreateSubMatrix is collective and
returns a matrix with global semantics (and where MatMult is actually
defined).


> -------
>  How is MatCreateBlockMat() which fortunately has the proper constructor
> related to MatCreateNest()? Do they serve the same purpose? Different
> purposes? Can they be merged together?
> --------
>

MatCreateBlockMat() is serial only and almost unused.


>
>   In the end we want to be able to run 'multphysics'  codes where it is a
> command line switch between 'block matrix' storage of the global Jacobian
> where each block is a submatrix of the entire matrix and stored itself for
> example in AIJ format; or one big honking AIJ matrix. In both cases the user
> calls MatSetValues() or MatSetValuesLocal() to put values in and has no
> switches in the code they write that depends on the matrix format, things
> like MatGetSubMatrix() and MatGetSubMatrices() work transparently on them.
> How do these four Matrix classes take us in that direction?
>

See snes/examples/tutorials/ex28.c.  It's a command line switch (read the
comments at the top of that file), subphysics modules get called with the
result of MatGetLocalSubMatrix() so identical sub-physics user-code is used
for single problems and for the coupled problem, when backed by any matrix
type (AIJ or Nest, Nest can have (S)BAIJ blocks).
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20110115/b9e0bec9/attachment.html>