[petsc-users] Fast LU solver for block diagonal matrix

[Matthew Knepley]

On Wed, Mar 18, 2015 at 2:10 PM, Luc Berger-Vergiat <lb2653 at columbia.edu>
wrote:

> I will have to do a MatMatMatMult or two MatMatMult since J is unsymmetric
> (Jeo != Joe^T).
> I think two MatMatMult could be more efficient if I use the
> MatMatMultSymbolic/MatMatMultNumeric correctly.
>

Test it first using

  -pc_fieldsplit_schur_precondition full

which does what you want, although maybe not in an optimal way:


https://bitbucket.org/petsc/petsc/src/f412d3e7cd17c738b1d11aea53d16ae27a4a5abb/src/ksp/pc/impls/fieldsplit/fieldsplit.c?at=master#cl-777

  Thanks,

    Matt

Best,
> Luc
>
>
> On 03/18/2015 04:04 PM, Barry Smith wrote:
>
>> On Mar 18, 2015, at 2:54 PM, Luc Berger-Vergiat <lb2653 at columbia.edu>
>>> wrote:
>>>
>>> Hi Barry,
>>> I would like to compute S explicitly, I have a few good option to
>>> precondition S but they are based on S not an approximation of S (i.e. I
>>> don't want to compute my preconditioner using Sp different from S).
>>>
>>> Also Jee is obtained using MatGetSubmatrix(J,isrow_e,iscol_e) and J is
>>> an AIJ matrix so I assume that Jee is AIJ too.
>>> I can convert Jee to BAIJ to take advantage of the block structure but
>>> from your conversation with Chung-Kan it might require to reassemble?
>>> In that case what about the following:
>>> MatCreateBAIJ(comm, 5, Jee_inv)
>>> for(i=0,i<nblocks,i++)
>>> {
>>>      MatGetValues(Jee,5,[5*i+0,5*i+1,5*i+2,5*i+3,5*i+4],5,[5*i+0,5*i+1,5*i+2,5*i+3,5*i+4],
>>> block_values)
>>>      some_package_inverts( block_values )
>>>      MatSetValuesBlocked(Jee_inv,5,[5*i+0,5*i+1,5*i+2,5*i+3,5*i+
>>> 4],5,[5*i+0,5*i+1,5*i+2,5*i+3,5*i+4], block_values)
>>> }
>>> MatAssemblyBegin()
>>> MatAssemblyEnd()
>>> With this I could then just go on and do matrix multiplications to
>>> finish my Schur complement calculations.
>>> Would this be decently efficient?
>>>
>>     This would be fine. Then you can use MatPtAP or MatRARt to compute
>> the triple matrix product.
>>
>>    Barry
>>
>>  Best,
>>> Luc
>>>
>>> On 03/18/2015 03:22 PM, Barry Smith wrote:
>>>
>>>>    Do you want to explicitly compute (as a matrix) S = Joo - joe *
>>>> inv(Jee) jeo   or do you want to just have an efficient computation of
>>>> S y  for any y vector?
>>>>
>>>>    Here is some possibly useful information. If you create Jee as a
>>>> BAIJ matrix of block size 5 and use MatILUFactor() it will efficiently
>>>> factor this matrix (each 5 by 5 block factorization is done with custom
>>>> code) then you can use MatSolve() efficiently with the result (note that
>>>> internally when factoring a BAIJ matrix PETSc actually stores the inverse
>>>> of the diagonal blocks so in your case the MatSolve() actually ends up
>>>> doing little matrix-vector products (and there are no triangular solves).
>>>>
>>>>    To use this with the MatCreateSchurComplement() object you can do
>>>>
>>>>     MatCreateSchurComplement(...,&S)
>>>>     MatSchurComplementGetKSP(S,&ksp)
>>>>     KSPSetType(ksp,KSPPREONLY);
>>>>
>>>>     now MatMult(S,y,z) will be efficient.
>>>>
>>>>     Of course you still have the question, how do you plan to solve S?
>>>> This depends on its structure and if you have a good way of preconditioning
>>>> it.
>>>>
>>>>     If you want to explicitly form S you can use MatMatSolve( fact,jeo)
>>>> but this requires making jeo dense which seems to defeat the purpose.
>>>>
>>>>    Barry
>>>>
>>>>
>>>>
>>>>
>>>>  On Mar 18, 2015, at 1:41 PM, Luc Berger-Vergiat <lb2653 at columbia.edu>
>>>>>   wrote:
>>>>>
>>>>> Hi all,
>>>>> I am solving multi-physics problem that leads to a jacobian of the
>>>>> form:
>>>>>
>>>>> [ Jee  Jeo ]
>>>>> [ Joe  Joo ]
>>>>>
>>>>> where Jee is 5by5 block diagonal. This feature makes it a very good
>>>>> candidate for a Schur complement.
>>>>> Indeed, Jee could be inverted in parallel with no inter-nodes
>>>>> communication.
>>>>> My only issue is the fact that the Schur complement is not accessible
>>>>> directly with PETSC, only an approximation is available, for direct solvers
>>>>> (usually S~Joo or S~Joo-Joe* diag(Jee)^-1 *Jeo).
>>>>>
>>>>> Any advice on how I could efficiently compute Jee^-1 for the given
>>>>> structure?
>>>>> I am currently thinking about hard coding the formula for the inverse
>>>>> of a 5by5 and sweeping through Jee (with some threading) and storing the
>>>>> inverse in-place. Instead of hard coding the formula for a 5by5 I could
>>>>> also do a MatLUFactorSym on a 5by5 matrix but it would not give me an
>>>>> inverse, only a factorization...
>>>>>
>>>>> Thanks in advance for your suggestions!
>>>>>
>>>>> --
>>>>> Best,
>>>>> Luc
>>>>>
>>>>>
>>>>>
>>>>>
>
>


-- 
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/20150318/60f33fb1/attachment-0001.html>