Barry Smith wrote: > > On May 20, 2009, at 3:47 PM, Jed Brown wrote: > >> Do you have a preferred interface for specifying a preconditioning >> matrix for A in Mat_SchurComplement > > I would add a new second argument to MatCreateSchurComplement() Ap > that would be the pmat passed into the KSPSetOperators(). To get the > current behavior one would use MatCreateSchurComplement(A,A,....) > >> and for the Schur complement in PC_FieldSplit? > > PCFieldSplitSchurPrecondition() could be modified to take an > additional argument that is used as the preconditioner matrix.
This is what I had in mind so I'll do it, but maybe not right away. >> It looks like I must implement MatGetSubMatrix for my matrix-free >> operators. This is of course only possible for very special IS. Is >> there any chance of a higher-level description of the submatrix? (I'm >> not sure ther is a better way, but as is, both PC_FieldSplit and my >> Mat_Shell already need to have duplicate representations of the >> submatrices but we'll be talking by complementing each IS.) >> > Not likely. The interface is to provide an IS. It will always be an > IS, but of course IS is an abstract base class and perhaps additional > ISs could be introduced that are "higher level" but in the end they > still have to be an abstract representation of a bunch of integers :-( That's fine, and it's only done once so even the naive representation isn't a performance issue. >> Do you have an interface in mind for specifying reduced forms for block >> preconditioners? For example, Galerkin discretization of incompressible >> flow gives >> >> J = [A B'] >> [B 0 ] >> >> The reduced forms >> >> [A ] or [A B'] >> [B S] [ S ] >> >> are appropriate for left- and right-preconditioned GMRES respectively >> (the preconditioned operator has minimal degree 2 and all unit >> eigenvalues). Also, >> >> [A ] >> [ S] >> >> is popular for symmetric problems. >> >> > Not sure what you mean "specifying". See the manual page for > Block_Preconditioners (bottom of fieldsplit.c) where it gives the > cases. Sadly it does not tell you directly how to run each different > case (It should) but I think you just need to use the proper > -pc_fieldsplit_type additive,multiplicitive,symmetric_multiplicative I've looked at that. This is different than what is described there because all of my examples are variants of a Schur-complement preconditioner where as those examples are block relaxation of some sort. As I see it, there are two sources of preconditioners for these systems, relaxation and factorization. So additive, multiplicative, and symmetric_multiplicative are all block relaxation where as Schur is the canonical example of factorization. What you have set up is nice and flexible for relaxation, but I'm having trouble seeing how to extend the factorization choices. Some more thoughts on this: So factorization always produces Schur complements and it is rare to get a scalable Schur solver without finding a way to interpret the Schur complement physically, or at least use some approximate commutators, to build a preconditioner. Since factorization of larger systems produces very complex Schur complements they will be very hard to understand, so I think it is useful to do both relaxation and factorization. There is a good chance that the right way to do this is to nest a factorization preconditioner as one block inside a larger split based on relaxation. Unfortunately, factorization preconditioners don't necessarily require preconditioners for diagonal blocks. For instance, consider the Jacobian in LNKS. This is a 3x3 (state, adjoint, design) system which is indefinite about the adjoint. The factorization you end up with needs preconditioners for the forward and adjoint models which sit at the (adjoint,state) and (state,adjoint) blocks. This case is not reducible to a 2x2 form and the factorization requires pivoting. The reduced Hessian (the block that ends up on the diagonal in the design space) is of course a nasty thing and you might not have a preconditioner available for it. Within the factorization variants, it is very common that you want to drop either the upper or lower factor. That is, suppose there is no pivoting and you have factors J = L D U. If the system is nonsymmetric, and sometimes if it is, you usually want to solve with (LD) or (DU), but not with the whole (LDU). This is because the effect of the dropping the factor is well-controlled by GMRES (usually the preconditioned operator has a low-degree minimal polynomial) and you need one more inner preconditioner application if you use the full factorization. I think it would be very useful to come up with an interface for this so that the user only needs to specify a factorization and provide preconditioners for the resulting Schur complements. Maybe such problems never really get bigger than 3x3 or 4x4 so they can be special-cased, but I think it ought to be possible to support any size and arbitrary factorizations. A related issue is telling the application (UserJacobian) what auxiliary systems are going to be needed by the preconditioner. For example, if you switch from factorization to relaxation, the application no longer needs to assemble an auxiliary system used to precondition the Schur complement. I can think of conventions like attaching a Mat with a special name to the (global, packed) preconditioning matrix (can be wrapped in an API call), but that would mean that the PC needs to have seen the matrix *before* SNESComputeFunction is called, and it just seems dirty and backwards. Probably the right thing is to have a PCFieldSplit API, so UserJacobian would unwrap the PC from SNES and query for which pieces are needed. (My motivation here is that we definitely don't want UserJacobian snooping -pc_fieldsplit_* to find out what will be needed.) For more general factorization, just identifying the blocks of the factorization seems tricky (at least coupled to the way you specify the factorization). Thanks for the discussion. 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/20090521/9264e44d/attachment.pgp>
