[petsc-dev] Jed!!!

[Eric Chamberland]

On 04/27/2014 12:32 PM, Jed Brown wrote:
 > Barry Smith <bsmith at mcs.anl.gov> writes:
 >
 >>    I am totally confused and likely users are to. Please point me to an
 >>    example that actually uses Matnest the way it is “supposed to be
 >>    used”.
 >
 > src/snes/examples/tutorials/ex28.c
 >
 > Note that this does not contain the word "nest".  That is as intended.
 >
 >>    I cannot even find a MatSetValues_Nest() so I do not see how one can
 >>    even interchange code between AIJ and Nest.
 >
 > You're supposed to use MatGetLocalSubMatrix().  Translating from global
 > indices to local indices is a disaster that we want to avoid.  So we go
 > the other way.  Speak the language of "split local" spaces during
 > assembly and the data structure itself can live wherever is most
 > efficient for the solver.

We have a question about this:

Where is the "disaster"?  In our mind, we think the  MatSetValues_Nest 
should:

1) do the assembly for local lines in the proper submatrices
2) communicate the non-local lines at MatAssemblyEnd() so that each 
process receives some lines which will be assembled by method 1)

So to do 1) (the assembly of local lines), you "just have to" convert 
column indices to:
a) retrieve which sub-matrices it belongs to
b) for each sub-matrix, do the assembly in that sub-matrix which 
involves converting lines/columns indices to -1 for elementary values 
not in that matrix.

We coded something (in sequential) in our code to do the work, and 
observed surprisingly that it is faster to do the assembly in a nested 
matrix:
3.97s vs 4.24s for a MatNest vs 357 369^2 CSR Matrix

(we suppose it could be faster because of the dichotomous search done in 
MatSetValues_SeqAIJ)

However, to do the correct work in parallel we think it should be 
implemented directly into something like MatSetValues_Nest... Have we 
missed something?

The use case scenario is the following: we have a P2-hierarchical (for 
example, a velocity field) which is made of a linear part and a 
quadratic "correction" part.  For a specific preconditioner 
(http://onlinelibrary.wiley.com/doi/10.1002/nla.757/abstract)
we want to do the "normal" assembly then work on the submatrices which 
are split by the linear vs quadratic parts of the velocity field.  The 
assembly is done "normally" in our code, meaning we don't want to 
duplicate the  computations done by the formulation (physical 
properties) but instead we want to have the elementary matrix to be 
split by the assembly functionality.

Thanks for your insights!

Eric


 >
 >>    The users manual has the word nest several times but no indication
 >>    of how to use the damn thing.
 >
 >  From the manual:
 >
 >    The key to format-independent assembly is the function
 >
 >    MatGetLocalSubMatrix(Mat A,IS isrow,IS iscol,Mat *submat);
 >