Ok - all good arguments... and I'm seeing some of your reasoning  
now.... but let me show you the use case I'm thinking of which might  
better explain why I think the LinearSolver link is useful...

I'm needing to solve individual systems _just_ for preconditioning  
purposes... and they use a preconditioner _as a solver_.  If a  
Preconditioner was _not_ a LinearSolver this would be what I would  
have to do (and is close to what I'm doing now... although there is  
even more manual manipulation).

ApplyPreconditioner(Vector x, Vector & y)
{
LinearImplicitSystem sys;
sys.add_variable("u");

AssemblePrecMatrix(sys.matrix);

copyCorrectPiecesOfX(x, sys.rhs);

Preconditioner pre(sys.matrix);
pre.apply(sys.rhs, sys.solution);

copyCorrectPiecesOfSolution(sys.solution, y);
}

This really doesn't look too bad... the problem is that there's  
datastructures that aren't used... or are redundant.  For instance  
there's a KSP object inside of that LinearImplicitSystem that never  
gets used... and it has it's own PC object... again that never gets  
used.  Then there's the fact that you have to call apply() and pass in  
two things that are both held within the system.

If you could just create the LinearSystem, passing in the  
Preconditioner then just call sys.solve().... it would simplify  
things.  I suppose that there's nothing stopping us from adding a  
precondition() method to LinearSystem.... that just calls the  
Preconditioner object (that you can pass in at construction time or  
attach).

What do you guys think?

Derek

On Jan 23, 2009, at 12:18 PM, John Peterson wrote:

> On Fri, Jan 23, 2009 at 12:01 PM, Derek Gaston <fried...@gmail.com>  
> wrote:
>> Alright!  Now we have a debate!  Let me see if I can persuade  
>> you... and if
>> I can't... I'm probably going to code it up anyway to give you a  
>> chance to
>> see what I mean.
>>
>> First, let's look at the two different forms of preconditioning  
>> (I'll be
>> using GMRES the whole time because it's easy to think about)...
>>
>> Solving: Ax=b
>>
>> Left Preconditioning: Find P^-1 Ax = P^-1b
>> Right Preconditioning: Find A P^-1 P x = b
>
> You're right, I wasn't taking right-preconditioning into account in my
> previous email.
>
>
>> Further, from Axelsson: Iterative Solution Methods (2000) pg. 472  
>> when
>> talking about left preconditioning GMRES:
>>
>> "The statement h:=C^-1 r is to be interpreted as solving the system  
>> Ch = r."
>>
>> Where h = y and C = P.
>
> What you've perfectly described is a Preconditioner object that *uses*
> a LinearSolver during its application, not something which "is a"
> LinearSolver.  Of course, this doesn't preclude you from also having a
> Preconditioner::solve() member, which just applies the preconditioner.
> It doesn't, in my opinion, make it a LinearSolver, in the sense in
> which we use that name in Libmesh.
>
>
>> You say that a Preconditioner doesn't satisfy the "is a"  
>> relationship for
>> inheritance... but I hope I've shown you that in fact it does.  A
>> Preconditioner "is a" LinearSolver that solves a specific Py=q  
>> system... the
>> solution of which is used internally by another LinearSolver.
>
> Perhaps, but in LibMesh a LinearSolver isn't a just thing that "solves
> a specific Py=q system," although it does do that.  It's an abstract
> interface to 3rd-party concrete linear algebra libraries as I
> mentioned before.  So, it's a little unfortunate that we have already
> taken the name that you would like to use but I still argue strongly
> against adding your proposed Preconditioner to this hierarchy.  We
> might need to have a discussion about expanding the meaning of
> LinearSolver in libmesh, but I think the current system is elegant and
> relatively simple to understand for new users.
>
> -- 
> John


------------------------------------------------------------------------------
This SF.net email is sponsored by:
SourcForge Community
SourceForge wants to tell your story.
http://p.sf.net/sfu/sf-spreadtheword
_______________________________________________
Libmesh-devel mailing list
Libmesh-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-devel

Reply via email to