[petsc-dev] MATNEST MATSUBMATRIX MATLOCALREF MATBLOCKMAT and friends

[Barry Smith]

On Jan 15, 2011, at 2:32 PM, Jed Brown wrote:

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

  I didn't mean a setter that allows setting one block at a time. I meant something like

    MatSetType(mat,MATNEST)
    MatNestSetBlocks(mat,IS isrows[], IS iscols[], Mat mats[]); 

   I see my name MatNestAddBlocks() was probably the cause of the confusion, since one would expect that one would call it multiple times.

   I agree there is no reason to make this dynamic, that is able to add blocks later.

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