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

[Barry Smith]

  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> 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> wrote:
>> 
>> Olivier Jamond <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/20201008/a50c74b7/attachment.html>