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

[Barry Smith]

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