[petsc-users] Use of MatRestrict/MatInterpolate with PCMG.

[Vijay S. Mahadevan]

>  Well if they happen to be equal then it will never apply the transpose thus giving a bad algorithm and garbage.

hmm, in the case when I provide a matrix for both
Interpolate/Restrict, I have (M==N==Nx==Ny),  this would always call
the MatMult routine. As long as I provide both these
operators/matrices explicitly, there is still no problem. A possible
issue is only when someone provides just the restriction or
prolongation. But I understand that when this happens, the other
operator is computed explicitly as its transpose. If this is actual
implementation, I still dont see a problem. Although, if the
restriction/prolongation operator are implicitly assumed to be
transpose of the other, then it will quite horribly fail since only
MatMult is called for both. I am not completely sure about the mode
currently used in petsc but it would be great if you can help me
understand.

Note: I probe more on this since my linear tesselation (most often)
results in the same number of dofs on the coarser level
(p-coarsened/h-refined) and I dont want a glitch to come back and bite
me later on..

Vijay

On Wed, Dec 15, 2010 at 8:36 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> On Dec 15, 2010, at 8:16 PM, Vijay S. Mahadevan wrote:
>
>> Barry,
>>
>> Thanks for the prompt change ! I do not work on the development
>> version but I can update these matrix routines alone.
>>
>>>  Note it can still glitch if the restricted size is exactly the original size. :-(
>>
>> Why would it glitch if the restricted size is the same as the original
>> size though ? I dont see a case where your check (M==Ny) would fail.
>> Can you please elaborate more on this ?
>
>  Well if they happen to be equal then it will never apply the transpose thus giving a bad algorithm and garbage.
>
>  Barry
>
>>
>> Vijay
>>
>> On Wed, Dec 15, 2010 at 8:04 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>
>>>  I have pushed this change to petsc-dev and it is ready for use.
>>>
>>>   Barry
>>>
>>>  Note it can still glitch if the restricted size is exactly the original size. :-(
>>>
>>>
>>> On Dec 15, 2010, at 7:53 PM, Barry Smith wrote:
>>>
>>>>
>>>>  Vijay,
>>>>
>>>>    The use of M>N in MatRestrict and MatInterpolate was always a bit cheesy since it has this broken case that you reported. I will change it to do as you suggest and use the size of the vectors in determining which way to apply. But note I will do this in petsc-dev http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html not petsc-3.1 so you'll need to switch if you are not using petsc-dev.
>>>>
>>>>   I'll try to get it down in the next few hours but it may take a little longer.
>>>>
>>>>
>>>>   Barry
>>>>
>>>> On Dec 15, 2010, at 6:06 PM, Vijay S. Mahadevan wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> I have an implementation issue with the MatRestrict/Interpolate
>>>>> functions. The problem is that one of my coarser levels (with PCMG)
>>>>> has higher dofs than the finest level. This does not always happen and
>>>>> requires a weird fine mesh system (in a sense) that uses multi-grid,
>>>>> but the idea is that the finest level problem has a high order (HO)
>>>>> discretization while the lower level mesh has a linear tesselation of
>>>>> the finest HO level (which I can optimize) and then adaptively
>>>>> coarsened levels beyond that. Since the number of columns in this case
>>>>> is larger than the number of rows, MatRestrict invariably calls
>>>>> MatMultTranspose to multiply instead of MatMult and vice-versa while
>>>>> calling  MatInterpolate. These result in assertion errors while
>>>>> comparing the length of Mat and Vec. The chosen method is based on
>>>>> whether (M>N) which seems to act against what I am doing here...
>>>>>
>>>>> I can always implement a shell matrix to replicate
>>>>> Restrict/Interpolate actions but my question is whether if such
>>>>> discretization will yield a consistent convergence in MG algorithm ?
>>>>> Is there a strong reason for checking if (M>N) rather than just doing
>>>>> (mat->rmap->N==y->map->N && mat->cmap->N==x->map->N) ? I would
>>>>> appreciate any detailed answer that you can provide for this and any
>>>>> suggestions to use the existing methods (without implementing the
>>>>> shell restriction) is very welcome.
>>>>>
>>>>> Thanks,
>>>>> vijay
>>>>
>>>
>>>
>
>