[petsc-users] VecGetSubVector

Matthew Knepley knepley at gmail.com
Sat Jan 15 09:43:07 CST 2011 

On Sat, Jan 15, 2011 at 8:23 AM, Gianluca Meneghello <gianmail at gmail.com>wrote:

> Dear Barry
>
> In order to use matrix-free methods with PCFIELDSPLIT, is it correct
> to start from
>
> MatShellSetOperation(mat,MATOP_MATGETSUBMATRIX, (void(*)(void))
> PetscErrorCode (*UserMatGetSubMatrix)(......));
>
> or does it require to dive deeper inside PETSc code?
>

That is fine. You must return a Mat object, probably a MatShell in your
case,
which applys the requested block of the operator.

   Matt


> Thanks
>
> Gianluca
>
> On 13 January 2011 22:16, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > On Jan 13, 2011, at 2:09 PM, Gianluca Meneghello wrote:
> >
> >> Barry,
> >>
> >> I'm afraid I didn't explain myself well in my first question.
> >>
> >> If I use PCFIELDSPLIT with
> PCFieldSplitSetType(pc,PC_COMPOSITE_ADDITIVE), is it the same as using
> PCASM?
> >
> >   Same in what sense? It solves a bunch of subproblems independently and
> adds together all the solutions. There can be overlapping in the fields or
> not depending how what you choose. The decomposition in ASM is by "geometry"
> while the decomposition in the PCFIELDSPLIT is between different "fields" or
> "types of variables". So yes they have many similarities.
> >
> >
> >>
> >> Concerning the MATSHELL to be used with PCFIELDSPLIT, is there an
> example from where to start from? I  guess I'm one of the people for which
> it is non-trivial :-)
> >
> >  Not really.
> >
> >   Barry
> >
> >>
> >> Thanks
> >>
> >> Gianluca
> >>
> >> Il giorno 13/gen/2011, alle ore 20.14, Barry Smith ha scritto:
> >>
> >>>
> >>> On Jan 13, 2011, at 12:00 PM, Gianluca Meneghello wrote:
> >>>
> >>>> Dear Barry and Jed,
> >>>>
> >>>> thanks again for your answers. I'm at the moment trying to understand
> more about how ASM and FieldSplit works. I've started reading Barry's book
> "Domain Decomposition, Parallel Multilevel Methods for Elliptic Partial
> Differential Equations". I hope it's the right starting point.
> >>>>
> >>>> Please let me know if you have further suggested readings, taking into
> account that I know nothing on domain decomposition in general --- but
> something in multigrid.
> >>>>
> >>>> Let me ask you a couple of questions for the moment:
> >>>> Is PCFIELDSPLIT additive the same as PCASM (I guess no as it does not
> use overlap)?
> >>>
> >>>  It can be additive or multiplicative depending on what you set with
> PCFieldSplitSetType()
> >>>
> >>>> Or is it a Substructuring Method (I haven't yet arrived to that
> chapter of the book!).
> >>>
> >>>  No, nothing to do with substructuring.
> >>>
> >>>> What's the difference between multiplicative and
> symmetric-multiplicative for PCFIELDSPLIT?  Does symmetric-multiplicative
> refers to eq 1.15 of "Domain Decomposition"?
> >>>
> >>> Yes
> >>>
> >>>>
> >>>> Is it possible use ASM and/or FieldSplit with a matrix-free method?
> >>>
> >>>  Yes, BUT the algorithms are coded around MatGetSubMatrix() and or
> MatGetSubMatrices() so to do matrix free you need to have code that applies
> "part" of the operator at a time (that is you cannot just have a matrix
> vector product that applies the entire operator to the entire vector. Once
> you have the ability to apply "part" of the operator at a time you need to
> code up a MATSHELL that responds appropriately to MatGetSubMatrix() and or
> MatGetSubMatrices() and returns new matrix-free shell matrices that apply
> only "their" part of the operator. This is non-trivial for many people but
> possible.
> >>>
> >>>  Barry
> >>>
> >>>
> >>>>
> >>>> Thanks
> >>>>
> >>>> Gianluca
> >>>>
> >>>> Il giorno 06/gen/2011, alle ore 21.55, Barry Smith ha scritto:
> >>>>
> >>>>>
> >>>>> On Jan 6, 2011, at 6:01 AM, Gianluca Meneghello wrote:
> >>>>>
> >>>>>> Dear Barry,
> >>>>>>
> >>>>>> thanks a lot for your answer.
> >>>>>>
> >>>>>> I tried to do some experiments with MatGetSubMatrix, but I guess I'm
> >>>>>> doing something wrong as it looks like being painfully slow.
> >>>>>
> >>>>> Hmm, we've always found the getsubmatrix takes a few percent of the
> time. Perhaps you are calling it repeatedly for the each domain, rather than
> once and reusing it? Also use MatGetSubMatrices() and get all the
> submatrices in one call rather than one at a time.
> >>>>>
> >>>>>>
> >>>>>> I changed approach and now I'm using the ASM preconditioner. What
> I'm
> >>>>>> actually trying to do is to split the domain in different parts ---
> >>>>>> like interior and boundaries --- and relax (solve) each one with a
> >>>>>> different smoother (solver). In your opinion, is this the right
> >>>>>> approach?
> >>>>>
> >>>>> Worth trying since it is easy. You can experiment with different
> smoothers on the subdomains using the -sub_pc_type etc options and set
> different prefixes for different subdomains.
> >>>>>
> >>>>>> So far it looks much faster than my previous approach of
> >>>>>> extracting each submatrix.
> >>>>>
> >>>>> ASM just uses MatGetSubMatrices() so shouldn't be faster or slower
> than a custom code that does the same thing.
> >>>>>
> >>>>>>
> >>>>>> Also, please let me ask you one more thing. When using ASM with
> >>>>>> different subdomains on the same process, is the order in which the
> >>>>>> domains are solved the same as the one in which they are stored in
> the
> >>>>>> IS array passed to PCASMSetLocalSubdomains()? I would be interested
> in
> >>>>>> controlling this in order to build a downstream marching smoother.
> >>>>>
> >>>>> It is only additive, there is no order as Jed noted. Doing
> multiplicative in general is tricky because you want to just update the
> parts of the residual that need to be updated.
> >>>>>
> >>>>>>
> >>>>>> Looking at the references, I've noticed you have worked on
> multigrid.
> >>>>>> What I'm trying to do is close to what is described in Diskin,
> Thomas,
> >>>>>> Mineck, "Textbook Multigrid Efficiency for Leading Edge Stagnation",
> >>>>>> in case you already know the paper.
> >>>>>>
> http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20040081104_2004082284.pdf
> >>>>>>
> >>>>>> Again, thanks a lot.
> >>>>>>
> >>>>>> Gianluca
> >>>>>>
> >>>>>> On 3 January 2011 17:43, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >>>>>>>
> >>>>>>> Gianluca,
> >>>>>>>
> >>>>>>> The expected use is with the VecScatter object. First you create a
> VecScatter object with VecScatterCreate() then each time you need the
> "subvector" you call VecScatterBegin() followed by VecScatterEnd() Note that
> usually the VecScatter object is retained and used many times.
> >>>>>>>
> >>>>>>> Barry
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>> On Jan 3, 2011, at 5:22 AM, Gianluca Meneghello wrote:
> >>>>>>>
> >>>>>>>> Hi,
> >>>>>>>>
> >>>>>>>> I'm new to PETSc, so that this can be a very simple question:
> >>>>>>>>
> >>>>>>>> I'm looking for something like VecGetSubVector, which I've seen it
> >>>>>>>> exists in the dev version but not in the released one.
> >>>>>>>>
> >>>>>>>> I need to write a smoother for a multigrid algorithm (something
> like a
> >>>>>>>> block Gauss Seidel) which can be written in matlab as
> >>>>>>>>
> >>>>>>>> for j = 1:ny
> >>>>>>>> P = <some matrix indices as function of j>;
> >>>>>>>> du(P) = L(P,P) \ (  rhs(P) - L(P,:)*du + D2(P,P)*du(P) );
> >>>>>>>> end
> >>>>>>>>
> >>>>>>>> where L is a matrix (in my case the linearized Navier Stokes).
> >>>>>>>>
> >>>>>>>> I was thinking about using IS for declaring P, so that D2(P,P) can
> be
> >>>>>>>> obtained using MatGetSubMatrix. I would need the same for the
> vector
> >>>>>>>> du.
> >>>>>>>>
> >>>>>>>> Is there a way to do that without using the developer version? (I
> >>>>>>>> really don't feel like being "experienced with building, using and
> >>>>>>>> debugging PETSc).
> >>>>>>>>
> >>>>>>>> Thanks in advance
> >>>>>>>>
> >>>>>>>> Gianluca
> >>>>>>>
> >>>>>>>
> >>>>>
> >>>>
> >>>
> >>
> >
> >
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110115/b809ebcb/attachment-0001.htm>

More information about the petsc-users mailing list