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 <[email protected]> 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 > >
