00001
00002
00003
00004
00005
00006 #include "PlayaAmesosSolver.hpp"
00007 #include "PlayaEpetraVector.hpp"
00008 #include "PlayaEpetraMatrix.hpp"
00009
00010
00011
00012 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00013 #include "PlayaVectorImpl.hpp"
00014 #include "PlayaLinearOperatorImpl.hpp"
00015 #include "PlayaLinearSolverImpl.hpp"
00016 #endif
00017
00018 #include "Amesos.h"
00019 #include "Amesos_BaseSolver.h"
00020
00021
00022 using namespace Teuchos;
00023
00024 namespace Playa
00025 {
00026
00027 AmesosSolver::AmesosSolver(const ParameterList& params)
00028 : LinearSolverBase<double>(params),
00029 kernel_()
00030 {
00031 if (parameters().isParameter("Kernel"))
00032 {
00033 kernel_ = getParameter<string>(parameters(), "Kernel");
00034 }
00035 else
00036 {
00037 kernel_ = "Klu";
00038 }
00039 }
00040
00041
00042
00043 SolverState<double> AmesosSolver::solve(const LinearOperator<double>& op,
00044 const Vector<double>& rhs,
00045 Vector<double>& soln) const
00046 {
00047 Playa::Vector<double> bCopy = rhs.copy();
00048 Playa::Vector<double> xCopy = rhs.copy();
00049
00050 Epetra_Vector* b = EpetraVector::getConcretePtr(bCopy);
00051 Epetra_Vector* x = EpetraVector::getConcretePtr(xCopy);
00052
00053 Epetra_CrsMatrix& A = EpetraMatrix::getConcrete(op);
00054
00055 Epetra_LinearProblem prob(&A, x, b);
00056
00057 Amesos amFactory;
00058 RCP<Amesos_BaseSolver> solver
00059 = rcp(amFactory.Create("Amesos_" + kernel_, prob));
00060 TEUCHOS_TEST_FOR_EXCEPTION(solver.get()==0, std::runtime_error,
00061 "AmesosSolver::solve() failed to instantiate "
00062 << kernel_ << "solver kernel");
00063
00064 int ierr = solver->Solve();
00065
00066 soln = xCopy;
00067
00068 SolverStatusCode state;
00069 std::string msg;
00070
00071 switch(ierr)
00072 {
00073 case 0:
00074 state = SolveConverged;
00075 msg = "converged";
00076 break;
00077 default:
00078 state = SolveCrashed;
00079 msg = "amesos failed: ierr=" + Teuchos::toString(ierr);
00080 }
00081
00082 SolverState<double> rtn(state, "Amesos solver " + msg, 0, 0);
00083 return rtn;
00084 }
00085
00086 }