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