[petsc-users] Ainsworth formula to solve saddle point problems / preconditioner for shell matrices

[Olivier Jamond]

Hi Barry,

Thanks for this work! I tried this branch with my code and sequential 
matrices on a small case: it does work!

Thanks a lot,
Olivier

On 09/10/2020 03:50, Barry Smith wrote:
>
>   Olivier,
>
>     The branch *barry/2020-10-08/invert-block-diagonal-aij* contains 
> an example src/mat/tests/ex178.c that shows how to compute inv(CC'). 
> It works for SeqAIJ matrices.
>
>     Please let us know if it works for you and then I will implement 
> the parallel version.
>
>   Barry
>
>
>> On Oct 8, 2020, at 3:59 PM, Barry Smith <bsmith at petsc.dev 
>> <mailto:bsmith at petsc.dev>> wrote:
>>
>>
>>  Olivier
>>
>>  I am working on extending the routines now and hopefully push a 
>> branch you can try fairly soon.
>>
>>  Barry
>>
>>
>>> On Oct 8, 2020, at 3:07 PM, Jed Brown <jed at jedbrown.org 
>>> <mailto:jed at jedbrown.org>> wrote:
>>>
>>> Olivier Jamond <olivier.jamond at cea.fr 
>>> <mailto:olivier.jamond at cea.fr>> writes:
>>>
>>>>>   Given the structure of C it seems you should just explicitly 
>>>>> construct Sp and use GAMG (or other preconditioners, even a direct 
>>>>> solver) directly on Sp. Trying to avoid explicitly forming Sp will 
>>>>> give you a much slower performing solving for what benefit? If C 
>>>>> was just some generic monster than forming Sp might be unrealistic 
>>>>> but in your case CCt is is block diagonal with tiny blocks which 
>>>>> means (C*Ct)^(-1) is block diagonal with tiny blocks (the blocks 
>>>>> are the inverses of the blocks of (C*Ct)).
>>>>>
>>>>>    Sp = Ct*C  + Qt * S * Q = Ct*C  +  [I - Ct * (C*Ct)^(-1)*C] S 
>>>>> [I - Ct * (C*Ct)^(-1)*C]
>>>>>
>>>>> [Ct * (C*Ct)^(-1)*C] will again be block diagonal with slightly 
>>>>> larger blocks.
>>>>>
>>>>> You can do D = (C*Ct) with MatMatMult() then write custom code 
>>>>> that zips through the diagonal blocks of D inverting all of them 
>>>>> to get iD then use MatPtAP applied to C and iD to get Ct * 
>>>>> (C*Ct)^(-1)*C then MatShift() to include the I then MatPtAP or 
>>>>> MatRAR to get [I - Ct * (C*Ct)^(-1)*C] S [I - Ct * (C*Ct)^(-1)*C] 
>>>>>  then finally MatAXPY() to get Sp. The complexity of each of the 
>>>>> Mat operations is very low because of the absurdly simple 
>>>>> structure of C and its descendants.   You might even be able to 
>>>>> just use MUMPS to give you the explicit inv(C*Ct) without writing 
>>>>> custom code to get iD.
>>>>
>>>> At this time, I didn't manage to compute iD=inv(C*Ct) without using
>>>> dense matrices, what may be a shame because all matrices are sparse 
>>>> . Is
>>>> it possible?
>>>>
>>>> And I get no idea of how to write code to manually zip through the
>>>> diagonal blocks of D to invert them...
>>>
>>> You could use MatInvertVariableBlockDiagonal(), which should perhaps 
>>> return a Mat instead of a raw array.
>>>
>>> If you have constant block sizes, MatInvertBlockDiagonalMat will 
>>> return a Mat.
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20201012/fd395c29/attachment.html>