# PlayaKrylovSolver.hpp

```00001 /* @HEADER@ */
00002 //
00004
00005 #ifndef PLAYA_KRYLOVSOLVER_HPP
00006 #define PLAYA_KRYLOVSOLVER_HPP
00007
00009 #include "PlayaIterativeSolver.hpp"
00010 #include "PlayaPreconditionerFactory.hpp"
00011 #include "PlayaILUKPreconditionerFactory.hpp"
00012 #include "PlayaSimpleComposedOpDecl.hpp"
00013
00014 namespace Playa
00015 {
00016 using namespace Teuchos;
00017
00021 template <class Scalar>
00022 class KrylovSolver : public IterativeSolver<Scalar>
00023 {
00024 public:
00026   KrylovSolver(const ParameterList& params);
00028   KrylovSolver(const ParameterList& params,
00029     const PreconditionerFactory<Scalar>& precond);
00030
00032   virtual ~KrylovSolver(){;}
00033
00035   virtual SolverState<Scalar> solve(const LinearOperator<Scalar>& op,
00036     const Vector<Scalar>& rhs,
00037     Vector<Scalar>& soln) const ;
00038 protected:
00039   virtual SolverState<Scalar> solveUnprec(const LinearOperator<Scalar>& op,
00040     const Vector<Scalar>& rhs,
00041     Vector<Scalar>& soln) const = 0 ;
00042
00043   const PreconditionerFactory<Scalar>& precond() const {return precond_;}
00044
00045 private:
00046   PreconditionerFactory<Scalar> precond_;
00047 };
00048
00049
00050 template <class Scalar> inline
00051 KrylovSolver<Scalar>::KrylovSolver(const ParameterList& params)
00052   : IterativeSolver<Scalar>(params), precond_()
00053 {
00054   if (!params.isParameter("Precond")) return;
00055
00056   const std::string& precondType = params.template get<string>("Precond");
00057
00058   if (precondType=="ILUK")
00059   {
00060     precond_ = new ILUKPreconditionerFactory<Scalar>(params);
00061   }
00062 }
00063
00064 template <class Scalar> inline
00065 KrylovSolver<Scalar>::KrylovSolver(const ParameterList& params,
00066   const PreconditionerFactory<Scalar>& precond)
00067   : IterativeSolver<Scalar>(params), precond_(precond)
00068 {
00069   TEUCHOS_TEST_FOR_EXCEPTION(params.isParameter("Precond"), std::runtime_error,
00070     "ambiguous preconditioner specification in "
00071     "KrylovSolver ctor: parameters specify "
00072     << params.template get<string>("Precond")
00073     << " but preconditioner argument is "
00074     << precond);
00075 }
00076
00077 template <class Scalar> inline
00078 SolverState<Scalar> KrylovSolver<Scalar>
00079 ::solve(const LinearOperator<Scalar>& op,
00080   const Vector<Scalar>& rhs,
00081   Vector<Scalar>& soln) const
00082 {
00083   if (precond_.ptr().get()==0)
00084   {
00085     return solveUnprec(op, rhs, soln);
00086   }
00087
00088
00089   Preconditioner<Scalar> p = precond_.createPreconditioner(op);
00090
00091   if (!p.hasRight())
00092   {
00093     LinearOperator<Scalar> A = p.left()*op;
00094     Vector<Scalar> newRHS = rhs.space().createMember();
00095     p.left().apply(rhs, newRHS);
00096     return solveUnprec(A, newRHS, soln);
00097   }
00098   else if (!p.hasLeft())
00099   {
00100     LinearOperator<Scalar> A = op * p.right();
00101     Vector<Scalar> intermediateSoln;
00102     SolverState<Scalar> rtn
00103       = solveUnprec(A, rhs, intermediateSoln);
00104     if (rtn.finalState()==SolveConverged)
00105     {
00106       p.right().apply(intermediateSoln, soln);
00107     }
00108     return rtn;
00109   }
00110   else
00111   {
00112     LinearOperator<Scalar> A = p.left() * op * p.right();
00113     Vector<Scalar> newRHS;
00114     p.left().apply(rhs, newRHS);
00115     Vector<Scalar> intermediateSoln;
00116     SolverState<Scalar> rtn
00117       = solveUnprec(A, newRHS, intermediateSoln);
00118     if (rtn.finalState()==SolveConverged)
00119     {
00120       p.right().apply(intermediateSoln, soln);
00121     }
00122     return rtn;
00123   }
00124 }
00125
00126 }
00127
00128 #endif
00129
```