00001
00002
00003
00004
00005
00006 #include "PlayaBelosSolver.hpp"
00007 #include "PlayaPreconditioner.hpp"
00008 #include "PlayaPreconditionerFactory.hpp"
00009 #include "PlayaParameterListPreconditionerFactory.hpp"
00010
00011
00012 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00013 #include "PlayaVectorImpl.hpp"
00014 #include "PlayaLinearOperatorImpl.hpp"
00015 #include "PlayaLinearSolverImpl.hpp"
00016 #endif
00017
00018 using namespace Playa;
00019 using namespace Teuchos;
00020
00021
00022 BelosSolver::BelosSolver(const ParameterList& params)
00023 : LinearSolverBase<double>(params), pf_(), hasSolver_(false)
00024 {
00025 if (params.isSublist("Preconditioner"))
00026 {
00027 ParameterList precParams = params.sublist("Preconditioner");
00028 pf_ = new ParameterListPreconditionerFactory(precParams);
00029 }
00030 }
00031
00032
00033
00034 SolverState<double> BelosSolver::solve(const LinearOperator<double>& A,
00035 const Vector<double>& rhs,
00036 Vector<double>& soln) const
00037 {
00038 typedef Anasazi::SimpleMV MV;
00039 typedef LinearOperator<double> OP;
00040 typedef Belos::LinearProblem<double, MV, OP> LP;
00041
00042 TEUCHOS_TEST_FOR_EXCEPT(!A.ptr().get());
00043 TEUCHOS_TEST_FOR_EXCEPT(!rhs.ptr().get());
00044
00045 if (!soln.ptr().get())
00046 {
00047 soln = rhs.copy();
00048
00049 soln.zero();
00050 }
00051
00052 if (rhs.norm2()==0.0)
00053 {
00054 soln.zero();
00055 SolverStatusCode code = SolveConverged;
00056 SolverState<double> state(code, "Detected trivial solution", 0, 0.0);
00057
00058 return state;
00059 }
00060
00061
00062 RCP<OP> APtr = rcp(new LinearOperator<double>(A));
00063 RCP<MV> bPtr = rcp(new MV(1));
00064 (*bPtr)[0] = rhs;
00065 RCP<MV> ansPtr = rcp(new MV(1));
00066 (*ansPtr)[0] = soln;
00067
00068
00069 RCP<LP> prob = rcp(new LP(APtr, ansPtr, bPtr));
00070
00071 TEUCHOS_TEST_FOR_EXCEPT(!prob->setProblem());
00072
00073
00074 if (pf_.ptr().get())
00075 {
00076 Preconditioner<double> P = pf_.createPreconditioner(A);
00077 if (P.hasLeft())
00078 {
00079 prob->setLeftPrec(rcp(new OP(P.left())));
00080 }
00081
00082 if (P.hasRight())
00083 {
00084 prob->setRightPrec(rcp(new OP(P.right())));
00085 }
00086 }
00087
00088 if (!hasSolver_)
00089 {
00090
00091 ParameterList plist = parameters();
00092
00093 RCP<ParameterList> belosList = rcp(&plist, false);
00094
00095 std::string solverType = parameters().get<string>("Method");
00096
00097 if (solverType=="GMRES")
00098 {
00099 solver_=rcp(new Belos::BlockGmresSolMgr<double, MV, OP>(prob, belosList));
00100 }
00101 else if (solverType=="CG")
00102 {
00103 solver_=rcp(new Belos::BlockCGSolMgr<double, MV, OP>(prob, belosList));
00104 }
00105 else if (solverType=="TFQMR")
00106 {
00107 solver_=rcp(new Belos::TFQMRSolMgr<double, MV, OP>(prob, belosList));
00108 }
00109 else if (solverType=="GCRODR")
00110 {
00111 solver_=rcp(new Belos::GCRODRSolMgr<double, MV, OP>(prob, belosList));
00112 hasSolver_ = true;
00113 }
00114 else if (solverType=="RCG")
00115 {
00116 solver_=rcp(new Belos::RCGSolMgr<double, MV, OP>(prob, belosList));
00117 hasSolver_ = true;
00118 }
00119 else
00120 {
00121 TEUCHOS_TEST_FOR_EXCEPT(!(solverType=="GMRES" || solverType=="CG"));
00122 }
00123 }
00124 else
00125 {
00126 solver_->setProblem( prob );
00127 }
00128
00129 Belos::ReturnType rtn = solver_->solve();
00130
00131 int numIters = solver_->getNumIters();
00132 double resid = solver_->achievedTol();
00133
00134 SolverStatusCode code = SolveFailedToConverge;
00135 if (rtn==Belos::Converged) code = SolveConverged;
00136 SolverState<double> state(code, "Belos solver completed", numIters, resid);
00137
00138 return state;
00139 }
00140
00141
00142