Matthew Knepley wrote: > On Wed, May 27, 2009 at 2:50 PM, Barry Smith <bsmith at mcs.anl.gov> wrote: > >> On May 27, 2009, at 2:23 PM, Jed Brown wrote: >> >> >>> Gee, having vectors with a little extra space at the end is hard! >>> >>> Notation: J is the Jacobian, P is the associated preconditioning matrix >>> >>> PCFieldSplit gets the full row blocks (Afield) or the off-diagonal >>> blocks (B,C) from J. I implement these matrix-free and use VecGhost >>> internally. It also queries P for diagonal blocks (A,D). I provide >>> these as assembled AIJ and BAIJ matrices. It then obtains work vectors >>> by getting vecs from A and D which returns ordinary MPI vectors (no >>> ghosts). >>> >> You could provide the getvecs() method for the A and D matrices that you >> provide to give you the ghosted vectors that you need? > > > I agree with Barry. If you want to use ghosts, make everything use VecGhost.
Thanks guys, this occurred to me while riding home. Hooray for dynamic polymorphism. I'm now using PCFieldSplit_Schur for non-Newtonian Stokes, with everything completely matrix-free except the preconditioners for A and S, so I'm finally phasing out my homegrown stuff. I would like to add an option -pc_fieldsplit_real_diagonal, to use diagonal blocks of J preconditioned by diagonal blocks of Jp, instead of only using blocks of Jp. I implemented this using an array 'mat' obatained in the same way as 'pmat', but you might have something else in mind. I am running into an issue with order of events. I want UserJacobian to be able to determine what assembly is needed, but I can't do things like call PCFieldSplitGetSubKSP because PCSetUp has not been called the first time around. The issue here is that I want to be able to optionally replace the PC in kspschur with my own. This is tricky because we can't create the Schur complement matrix before PCSetUp, so there won't be a subksp[0] and PCFieldSplitGetSubKSP fails. Creating the KSP before the Schur matrix exists creates an ownership issue that would add complexity so I'm hesitant to do it unless it proves necessary. For the interface in these cases, I've been thinking about what information to put in the preconditioning matrix (Jp and submatrices) versus communicating directly to the PC. My inclination is to put all real data in Jp because it's messy for UserJacobian to talk directly to the PC, and just move metadata (setting options, usually prior to PCSetFromOptions) directly to the PC. That is, I like to think of Jp, less as a "matrix", and more as "data" that will be used by the preconditioner. It just so happens that an assembled matrix is a really useful format (understood by almost all PCs) to store that "data" so basing it on the Mat interface is appropriate. To illustrate, consider the usual problem J = [A B; C D], S = D-C*inv(A)*B. I don't want PCFieldSplit to need to do anything special based on how S gets preconditioned, and I don't want UserJacobian to need to unwrap sub-KSPs from PCFieldSplit (creation issues aside, this is a dependency that I think is avoidable). But we have this Jp hanging around and it's usual purpose is to route information from UserJacobian through the KSP to the PC so maybe we can also use it to route information through PCFieldSplit (via submatrices). The query for Ap is obvious, but what does a query for Dp mean? The matrix that actually sits there (often 0) is not actually useful. One choice is to return some other user-assembled matrix Sp (e.g. mass matrix for Stokes). I think this is a slight abuse of the interface since it specializes the definition of Jp to a particular preconditioner. In particular, it won't work correctly with FieldSplit relaxation (but that would fail for the class of problems we're discussing anyway). Let's assume that PCFieldSplit gets Sp somehow, perhaps by just by querying Dp. Now what if our preconditioner for S needs more information? For example, unscaled BFBt amounts to bfbt(S,Sp) = inv(Bp Cp) B A C inv(Bp Cp) where the product Bp*Cp is actually formed. Now "Sp" is no longer needed as an assembled matrix, but assembled variants of Bp and Cp are needed (B, A, and C are part of S). An alternative is to recognize that B*C corresponds to a Laplacian L, and actually assemble this as Lp. Let's call this variation the least squares commutator. It's normally accompanied by a couple scaling vectors, but the unscaled version looks like lsc(S,Sp) = ksp(L,Lp) B A C ksp(L,Lp) How should we get L, Lp, and scaling vectors through the layers of PCFieldSplit? Another preconditioner for S is from Kay, Loghin, and Wathen, klw(S,Sp) = inv(Mp) Ac ksp(L,Lp) where Mp is a (lumped) pressure mass matrix and Ac is a representation of A in the pressure space. Now we need to get Mp, Ac, L, Lp through. At this point, my preference is to compose this stuff with Sp. This avoids any need to unwrap internals of PCFieldSplit and PCFieldSplit does not need to concern itself with what auxiliary information the preconditioner for S might need. It would also make it easy to write alternate preconditioners for S, for example PCBFBt could be greatly simplified, and lots of other variants could be created easily. If we do indeed go this route, I see two remaining tricky bits. First, how to identify the various "Sp" showing up in a factorization? When pivoting is used, they are not uniquely identified by where they sit (coordinates) in the factored matrix. Second, how to avoid computing way more than needed in UserJacobian. For example, if I want to support default, bfbt, lsc, and klw, I should be able to provide Sp, Bp, Cp, L, Lp, Mp, Ac, and scaling, but not all of these will be used in any particular run. Before PCSetUp_FieldSplit has been called, the submatrices don't even exist and none of the wiring for the preconditioners has been set up. One option is to assemble everything the first time around, but this could could blow your memory. Another is to assemble lazily (e.g. in MatGetSubMatrix) at least the first time around (it's usually much faster to assemble several matrices at once). This goofs with profiling, but the first time is often nonsense anyway. Jed -------------- next part -------------- A non-text attachment was scrubbed... Name: signature.asc Type: application/pgp-signature Size: 260 bytes Desc: OpenPGP digital signature URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20090529/cdae0181/attachment.pgp>