## Re: [Libmesh-devel] Preconditioner Class

 Re: [Libmesh-devel] Preconditioner Class From: Derek Gaston - 2009-01-23 18:02:41 ```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 speciﬁed 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 > 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 ```

 [Libmesh-devel] Preconditioner Class From: Derek Gaston - 2009-01-23 02:38:45 Attachments: Message as HTML ```Ok... so I have a huge email already written that has a big long justification of a Preconditioner Class. But after writing it I realized that it was long winded and no one would really care about the discovery process I went through during the last few days.... so I've decided to just skip to the meat. If you want the justification later let me know.... We need a Preconditioner class. It will allow you to setup a preconditioner and use it multiple times with different RHS vectors. It should inherit from LinearSolver, using the matrix and rhs vector to produce the preconditioned vector in solution whenever solve() is called. Concrete implementations of this class would be: PetscPreconditioner MLPreconditioner (from Trilinos) ShellPreconditioner etc.. I also wouldn't mind a PetscHyprePreconditioner where both MLPreconditioner and it inherit from MultigridPreconditioner. This would allow simple manipulation of common multigrid preconditioner paramters (like smoothing, cycle numbers, etc.). But this can be negotiated. LinearImplicitSystem should allow for attaching a preconditioner class that would then be used by the solver. LinearImplicitSystem should also allow you to specify a Preconditioner _as_ the solver itself... or maybe there should be a separate PreconditionerSystem... I can't decide. 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). 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... Derek ```
 Re: [Libmesh-devel] Preconditioner Class From: John Peterson - 2009-01-23 03:42:17 ```On Thu, Jan 22, 2009 at 7:40 PM, Derek Gaston 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 ```
 Re: [Libmesh-devel] Preconditioner Class From: Derek Gaston - 2009-01-23 18:02:41 ```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 speciﬁed 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 > 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 ```
 Re: [Libmesh-devel] Preconditioner Class From: Roy Stogner - 2009-01-23 18:22:27 ```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 ```
 Re: [Libmesh-devel] Preconditioner Class From: Kirk, Benjamin (JSC-EG) - 2009-01-23 19:05:09 ```>> which starts to >> look more like: Preconditioner.solve(Matrix, q, y). > > Or: p(Matrix); p.matvec(q, y). p(Matrix); p.matvec(q, y). Or both... Preconditioner p(Matrix); p.matvec(q, y); ... (keep using Matrix and whatever you have done to it...) p.matvec(q2,y2); p.matvec(Matrix,q,y); ...(reset the matrix, do whatever you have to do...) Can't we all just get along?? -Ben ```
 Re: [Libmesh-devel] Preconditioner Class From: John Peterson - 2009-01-23 19:18:25 ```On Fri, Jan 23, 2009 at 12:01 PM, Derek Gaston 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 ```
 Re: [Libmesh-devel] Preconditioner Class From: Derek Gaston - 2009-01-23 21:20:57 ```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 > 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 ```
 Re: [Libmesh-devel] Preconditioner Class From: Derek Gaston - 2009-01-23 21:25:51 ```I think I'm going to go ahead and start on a Preconditioner base class... that doesn't inherit from LinearSolver for now. I'll get it up and running with a PetscPreconditioner class. This will pretty much meet my requirements for now. Derek On Jan 23, 2009, at 2:19 PM, Derek Gaston wrote: > 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 >> 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 > ```
 Re: [Libmesh-devel] Preconditioner Class From: John Peterson - 2009-01-23 22:31:27 ```On Fri, Jan 23, 2009 at 3:19 PM, Derek Gaston wrote: > 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. I know that's a simple example but it really doesn't look too bad. It seems that x, and y are only temporary, so if they became members of Preconditioner and you called the ApplyPreconditioner function solve() instead, I think it would be very close to what you want... Also, I wouldn't worry too much about unused PC/KSP objects. They can't possibly take up much memory or take much time to build up/destroy, can they?. In any event, an ImplicitSystem may be more applicable in your example (assuming it's not a virtual base class, I didn't check!) if System::solve() is *never* to be called. The System hierarchy doesn't get a LinearSolver until the LinearImplicitSystem level. -- John ```
 Re: [Libmesh-devel] Preconditioner Class From: Derek Gaston - 2009-01-23 22:27:52 ```On Jan 23, 2009, at 3:02 PM, John Peterson wrote: > In any event, an ImplicitSystem may be more > applicable in your example (assuming it's not a virtual base class, I > didn't check!) if System::solve() is *never* to be called. The System > hierarchy doesn't get a LinearSolver until the LinearImplicitSystem > level. If ImplicitSystem doesn't have a pure virtual hanging around... then you are absolutely right.... that's what I should be using. Ok... well I'm going to make this completely detached Preconditioner object for now. We'll see later if we want to tie it to something. I've actually got most of the class already put together... but it looks like I'm not going to get to finish it until Monday. Thanks for the discussion! Derek ```