On Tue, Jan 20, 2009 at 09:58, John Peterson <jwpeterson@...> wrote:
>> Also... on a related question.... what do you guys think about
>> providing a Preconditioner base class in libMesh? Concrete
>> implementations of this class would be HyprePreconditioner,
>> MLPreconditioner, LinearSolvePreconditioner, etc... I guess they are
>> similar to linear solvers... but not exactly the same....
> I guess I don't see the need for this, but that's not to say I can't
> be convinced. There's nothing like that in PETSc AFAIK.
There sure is, but the linear solve is actually a method for Mat. In
OO language, all PETSc objects come from a factory which can be
influenced by command line parameters. When you use -pc_type hypre,
you make the factory create an instance of the derived class PC_HYPRE.
At least that is how a similar API would be implemented in C++ so
that the user never hard-codes concrete types. In the current
interface, the factories are not bona-fide objects, but it's a
reasonable way to think of PETSc objects before XXXSetFromOptions or
whatever is called. The trouble with implementing this on the Libmesh
side is that you duplicate a lot of work and lose significant
functionality unless you reimpliment a lot of PETSc. PCSHELL is
always there for when you have a clever way to form a preconditioner.
I don't know if it is useful for anyone on this list, but I like to
think of forming a preconditioner as
PCApply = PC(A,Ap)
where A is an implementation of MatMult which is frequently ignored
and Ap is ``extra information'' which is usually an assembled sparse
matrix, but could be any information allowing you to create a
(nominally) linear operation PCApply so that
PCApply \circ A
A \circ PCApply
has favorable properties for the Krylov solver being used ($\circ$ is
function composition). This view makes it clear that the
preconditioner is intimately coupled to the format of Ap. PETSc
handles this by implementing parts of many preconditioners as methods
of Mat (factor, solve, relax, etc). Many third party preconditioners
require their own matrix type, usually derived from AIJ.