52#include "Thyra_BelosLinearOpWithSolveFactory.hpp"
53#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
54#include "Thyra_EpetraLinearOp.hpp"
57#include "Epetra_SerialComm.h"
58#include "EpetraExt_readEpetraLinearSystem.h"
61#ifdef HAVE_BELOS_IFPACK
66#include "Teuchos_ParameterList.hpp"
67#include "Teuchos_VerboseObject.hpp"
69int main(
int argc,
char* argv[])
74 Teuchos::RCP<Teuchos::FancyOStream>
75 out = Teuchos::VerboseObjectBase::getDefaultOStream();
80 int maxIterations = 400;
82 int gmresKrylovLength = 25;
83 int outputFrequency = 1;
84 bool outputMaxResOnly =
true;
85 double maxResid = 1e-6;
87 Teuchos::RCP<Teuchos::ParameterList>
88 belosLOWSFPL = Teuchos::rcp(
new Teuchos::ParameterList() );
90 belosLOWSFPL->set(
"Solver Type",
"Block GMRES");
92 Teuchos::ParameterList& belosLOWSFPL_solver =
93 belosLOWSFPL->sublist(
"Solver Types");
95 Teuchos::ParameterList& belosLOWSFPL_gmres =
96 belosLOWSFPL_solver.sublist(
"Block GMRES");
98 belosLOWSFPL_gmres.set(
"Maximum Iterations",
int(maxIterations));
99 belosLOWSFPL_gmres.set(
"Convergence Tolerance",
double(maxResid));
100 belosLOWSFPL_gmres.set(
"Maximum Restarts",
int(maxRestarts));
101 belosLOWSFPL_gmres.set(
"Block Size",
int(blockSize));
102 belosLOWSFPL_gmres.set(
"Num Blocks",
int(gmresKrylovLength));
103 belosLOWSFPL_gmres.set(
"Output Frequency",
int(outputFrequency));
104 belosLOWSFPL_gmres.set(
"Show Maximum Residual Norm Only",
bool(outputMaxResOnly));
106#ifdef HAVE_BELOS_IFPACK
110 Teuchos::ParameterList &ifpackPFSL = belosLOWSFPL->sublist(
"IfpackPreconditionerFactory");
112 ifpackPFSL.set(
"Overlap",
int(2));
113 ifpackPFSL.set(
"Prec Type",
"ILUT");
124 std::string matrixFile =
"orsirr1.hb";
127 Epetra_SerialComm comm;
128 Teuchos::RCP<Epetra_CrsMatrix> epetra_A;
129 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
132 Teuchos::RCP<const Thyra::LinearOpBase<double> >
133 A = Thyra::epetraLinearOp(epetra_A);
136 Teuchos::RCP<const Thyra::VectorSpaceBase<double> > domain = A->domain();
139 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
142#ifdef HAVE_BELOS_IFPACK
145 belosLOWSFactory->setPreconditionerFactory(
147 ,
"IfpackPreconditionerFactory"
152 belosLOWSFactory->setParameterList( belosLOWSFPL );
155 belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
158 Teuchos::RCP<Thyra::LinearOpWithSolveBase<double> >
159 nsA = belosLOWSFactory->createOp();
162 Thyra::initializeOp<double>( *belosLOWSFactory, A, &*nsA );
165 Teuchos::RCP< Thyra::MultiVectorBase<double> >
166 b = Thyra::createMembers(domain, numRhs);
167 Thyra::seed_randomize<double>(0);
168 Thyra::randomize(-1.0, 1.0, &*b);
171 Teuchos::RCP< Thyra::MultiVectorBase<double> >
172 x = Thyra::createMembers(domain, numRhs);
173 Thyra::assign(&*x, 0.0);
177 Thyra::SolveStatus<double> solveStatus;
178 solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
181 *out <<
"\nBelos LOWS Status: "<< solveStatus << std::endl;
186 std::vector<double> norm_b(numRhs), norm_res(numRhs);
187 Teuchos::RCP< Thyra::MultiVectorBase<double> >
188 y = Thyra::createMembers(domain, numRhs);
191 Thyra::norms_2( *b, &norm_b[0] );
194 A->apply( Thyra::NONCONJ_ELE, *x, &*y );
197 Thyra::update( -1.0, *b, &*y );
200 Thyra::norms_2( *y, &norm_res[0] );
203 double rel_res = 0.0;
204 *out <<
"Final relative residual norms" << std::endl;
205 for (
int i=0; i<numRhs; ++i) {
206 rel_res = norm_res[i]/norm_b[i];
207 if (rel_res > maxResid)
209 *out <<
"RHS " << i+1 <<
" : "
210 << std::setw(16) << std::right << rel_res << std::endl;
213 return ( success ? 0 : 1 );
int main(int argc, char *argv[])
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
Concrete preconditioner factory subclass based on Ifpack.