00001
00002
00003
00004
00005 #ifndef PLAYA_KRYLOVSOLVER_HPP
00006 #define PLAYA_KRYLOVSOLVER_HPP
00007
00008 #include "PlayaDefs.hpp"
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
00018
00019
00020
00021 template <class Scalar>
00022 class KrylovSolver : public IterativeSolver<Scalar>
00023 {
00024 public:
00025
00026 KrylovSolver(const ParameterList& params);
00027
00028 KrylovSolver(const ParameterList& params,
00029 const PreconditionerFactory<Scalar>& precond);
00030
00031
00032 virtual ~KrylovSolver(){;}
00033
00034
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