Derek Gaston wrote:
> let's look at what the function parameters would be for  
> preconditioning a vector.  You wouldn't simply call  
> Preconditioner.matvec(q, y) as you might believe.

Yes I would.

> The trouble is that  
> you have to provide a Matrix to be preconditioned....

Yes, in the preconditioner constructor.  Many preconditioners would 
royally suck if they had to be able to deal with a different matrix 
every time you applied them.


> which starts to  
> look more like: Preconditioner.solve(Matrix, q, y).

Or: Preconditioner p(Matrix); p.matvec(q, y).

> Essentially a  
> Preconditioner doesn't resemble a matvec... because it needs to take a  
> matrix as an argument.... which doesn't make any sense for a matvec...  

Sure it does: A^-1*q ~= invP(A)*q
Here your matrix is a function of another matrix.  And since it's a 
function in the mathematical "same input => same output" sense rather 
than just the lazy C++ "take some stuff and do some stuff to it based on 
possible global state with possible side effects" sense, it's a good 
idea to cache the evaluation of that function by, for example doing it 
in the Preconditioner constructor.

> Further... if the Preconditioner is just a matvec (maybe you really do  
> have some factorization of the matrix stored away) there is nothing  
> awkward about calling Preconditioner.solve(Matrix, q, y)

Except when you have to repeatedly check that the Matrix hasn't changed 
to see if you have to rebuild your internal data.

> if  
> you treat a preconditioner that really does a solve as a matvec it  
> feels more awkward to do Preconditioner.matvec(q, y) where internally  
> it's not doing a matvec at all but really doing a solve using a matrix  
> that you must have setup beforehand.

Perhaps it's just the name "matvec" that's throwing you?  Think of it as 
a shell matrix - it's just an operator that takes one vector and applies 
a linear transformation to produce another vector.

Of course, by this expansive definition a LinearSolver is also a shell 
matrix (assuming we ignore the nonlinearities introduced by the 
iterative algorithm), so the questions then are just: is every 
Preconditioner a LinearSolver, and/or is every LinearSolver a 
Preconditioner?

I say "no" to the first, "yes" to the second.  Not all preconditioners 
need to explicitly store or even implicitly know what "P" is in order to 
apply the action of P^-1.  If you happen to have a nicely defined "P", 
you can make a SparseMatrix or a shell matrix out of it, hand that to a 
LinearSolver, and wrap a shell matrix interface around that solver to 
get your matvec... but you can't always (efficiently!) go the other way 
around.
---
Roy

------------------------------------------------------------------------------
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
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-devel

Reply via email to