On Wed, Mar 18, 2015 at 2:10 PM, Luc Berger-Vergiat <[email protected]> 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 <[email protected]> >>> 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 <[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 >>>>> >>>>> >>>>> >>>>> > > -- 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
