Sundance 2.1

00001 /* @HEADER@ */
00002 //
00003  /* @HEADER@ */
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()) soln = rhs.copy();
00046 
00047   RCP<OP> APtr = rcp(new LinearOperator<double>(A));
00048   RCP<MV> bPtr = rcp(new MV(1));
00049   (*bPtr)[0] = rhs;
00050   RCP<MV> ansPtr = rcp(new MV(1));
00051   (*ansPtr)[0] = soln;
00052   
00053   
00054   RCP<LP> prob = rcp(new LP(APtr, ansPtr, bPtr));
00055 
00056   TEUCHOS_TEST_FOR_EXCEPT(!prob->setProblem());
00057 
00058   
00059   if (pf_.ptr().get())
00060   {
00061     Preconditioner<double> P = pf_.createPreconditioner(A);
00062     if (P.hasLeft())
00063     {
00064       prob->setLeftPrec(rcp(new OP(P.left())));
00065     }
00066   
00067     if (P.hasRight())
00068     {
00069       prob->setRightPrec(rcp(new OP(P.right())));
00070     }
00071   }
00072 
00073   if (!hasSolver_)
00074     {
00075 
00076       ParameterList plist = parameters();
00077 
00078       RCP<ParameterList> belosList = rcp(&plist, false);
00079 
00080       std::string solverType = parameters().get<string>("Method");
00081       
00082       if (solverType=="GMRES")
00083   {
00084     solver_=rcp(new Belos::BlockGmresSolMgr<double, MV, OP>(prob, belosList));
00085   }
00086       else if (solverType=="CG")
00087   {
00088     solver_=rcp(new Belos::BlockCGSolMgr<double, MV, OP>(prob, belosList));
00089   }
00090       else if (solverType=="TFQMR")
00091   {
00092     solver_=rcp(new Belos::TFQMRSolMgr<double, MV, OP>(prob, belosList));
00093   }
00094       else if (solverType=="GCRODR")
00095   {
00096     solver_=rcp(new Belos::GCRODRSolMgr<double, MV, OP>(prob, belosList));
00097     hasSolver_ = true; // only cache recycling solvers
00098   }
00099       else if (solverType=="RCG")
00100   {
00101     solver_=rcp(new Belos::RCGSolMgr<double, MV, OP>(prob, belosList));
00102     hasSolver_ = true; // only cache recycling solvers
00103   }
00104       else
00105   {
00106     TEUCHOS_TEST_FOR_EXCEPT(!(solverType=="GMRES" || solverType=="CG"));
00107   }
00108     }
00109   else // reset problem
00110     {
00111       solver_->setProblem( prob );
00112     }
00113   
00114   Belos::ReturnType rtn = solver_->solve();
00115 
00116   int numIters = solver_->getNumIters();
00117   double resid = -1.0;
00118   
00119   SolverStatusCode code = SolveFailedToConverge;
00120   if (rtn==Belos::Converged) code = SolveConverged;
00121   SolverState<double> state(code, "Belos solver completed", numIters, resid);
00122   
00123   return state;
00124 }
00125 
00126 
00127
PlayaBelosSolver.cpp
