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 Where P~A and (more importantly) P^-1 ~ A^-1 In GMRES the things you need to calculate with each type of preconditioning are: Left: P^-1 Ar Right: AP^-1r (which as far as we're concerned is just: P^1r) Where r is the current residual. Let's make the substitution of q = Ar (Left) or q = r (Right)... therefore we need to calculate P^-1 q. At first glance this would appear to substantiate your claims that a preconditioner is a matvec. But (for nontrivial preconditioners) it's not that simple. In reality you never have or form P^-1. Instead you (usually approximately) _solve_ the system: Py = q Which of course means y ~= P^-1 q. Essentially you approximate the action of the inverse of the preconditioning matrix.... by _doing a linear solve_. Indeed in the Petsc documentation they have this to say: (In section 16.4: Unimportant Details of PC) "In particular, for a preconditioner matrix, B, that has been set via PCSetOperators(pc,A,B,flag), the routine PCApply(pc,x,y) computes y = B^−1 x by solving the linear system By = x with the specified preconditioner method." Where their "x" matches the "q" I was using above. 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. Now 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. The trouble is that you have to provide a Matrix to be preconditioned.... which starts to look more like: Preconditioner.solve(Matrix, 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... but makes perfect sense for a solve. 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) and internally doing y = P^-1q. However the opposite is not true... 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. Finally, I think that the polymorphism works out nicely to make a Preconditioner a LinearSolver... that way you can actually _solve_ (or approximately solve) systems with it if you please (which is exactly what I'm going to be doing... because I setup separate Systems for each piece of my preconditioner and "solve" them using a (possibly different for each piece) preconditioner). 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. Derek On Jan 22, 2009, at 8:42 PM, John Peterson wrote: > On Thu, Jan 22, 2009 at 7:40 PM, Derek Gaston <fried...@gmail.com> > wrote: >> >> We need a Preconditioner class. It will allow you to setup a >> preconditioner >> and use it multiple times with different RHS vectors. > > The linear solver interface seems to be already general enough to do > this for solves; you can specify a matrix for preconditioning and you > can re-use it multiple times. Granted we don't have a separate way of > just applying a preconditioner outside of a solve scenario, but this > is essentially always a matrix-vector product, or whatever shell > routine you have specified to perform that action, and we have generic > interfaces for matrix-vector products already. > >> It should inherit from LinearSolver, using the matrix and rhs >> vector to > > This I don't agree with at all, it conflates the idea of an object > whose responsibility is to solve a linear system and an object which > is a component part of a linear solve. The LinearSolver class has > always been intended as a generic interface class for hiding the > implementation details of concrete sparse linear solver instantiations > (petsc, trilinos, laspack) and I don't think what you are describing > "is a" LinearSolver in that sense. > > >> Just to give an example of a Preconditioner: the >> PetscPreconditioner would >> create a PC struct calling PCCreate(), PCSetType() and >> PCSetOperators() for >> you. Then when solve is called it would do a PCApply(pc, rhs, >> solution). > > The PCSetType part duplicates existing functionality in the > PetscLinearSolver class, but it should be possible to unify that and > other duplicated things in a nice way. > >> Looking over ML (in Trilinos) it would have similar behavior. >> To sum up the justification of these classes: I have to do a _lot_ >> of manual >> petsc and libmesh manipulation to setup separate preconditioning >> systems and >> _efficiently_ solve them.... and a lot of that manipulation is >> going to have >> to be duplicated for using ML in Trilinos. >> So what do you guys think? If I get the thumbs up I could probably >> have >> something in the repository by the end of the day tomorrow... > > The functionality you are describing does not appear to me to be too > far outside the ShellMatrix interface class we had an extensive thread > on previously. Even if you don't want your Preconditioner class to > derive from a NumericMatrix of some kind, don't you agree that it will > have a lot of "matrix-like" interfaces? In my opinion you should make > Preconditioner an abstract base class, and make it usable within the > existing LinearSolver interface, not derive from it. > > > -- > 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
