[petsc-users] VecGetSubVector

[Barry Smith]

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