53 #include "Thyra_BelosLinearOpWithSolveFactory.hpp" 54 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" 55 #include "Thyra_TpetraLinearOp.hpp" 58 #include "Tpetra_ElementSpace.hpp" 59 #include "Tpetra_VectorSpace.hpp" 60 #include "Tpetra_CisMatrix.hpp" 63 #ifdef HAVE_BELOS_TRIUTILS 68 #include "Teuchos_ParameterList.hpp" 69 #include "Teuchos_VerboseObject.hpp" 70 #include "Teuchos_GlobalMPISession.hpp" 73 # include "Tpetra_MpiPlatform.hpp" 75 # include "Tpetra_SerialPlatform.hpp" 81 int main(
int argc,
char* argv[])
86 Teuchos::RCP<Teuchos::FancyOStream>
87 out = Teuchos::VerboseObjectBase::getDefaultOStream();
89 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
92 typedef std::complex<double>
ST;
94 typedef std::complex<double>
ST;
99 typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;
101 ST one = Teuchos::ScalarTraits<ST>::one();
102 ST zero = Teuchos::ScalarTraits<ST>::zero();
105 MPI_Comm mpiComm = MPI_COMM_WORLD;
106 const Tpetra::MpiPlatform<OT,OT> ordinalPlatform(mpiComm);
107 const Tpetra::MpiPlatform<OT,ST> scalarPlatform(mpiComm);
109 const Tpetra::SerialPlatform<OT,OT> ordinalPlatform;
110 const Tpetra::SerialPlatform<OT,ST> scalarPlatform;
118 std::string matrixFile =
"mhd1280b.cua";
126 info = readHB_newmat_double(matrixFile.c_str(),&dim,&dim2,&nnz,
127 &colptr,&rowind,&dvals);
129 if (info == 0 || nnz < 0) {
130 *out <<
"Error reading '" << matrixFile <<
"'" << std::endl;
138 for (
int ii=0; ii<nnz; ii++) {
139 cvals[ii] =
ST(dvals[ii*2],dvals[ii*2+1]);
146 const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
147 const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
150 RCP<Tpetra::Operator<OT,ST> >
151 tpetra_A = rcp(
new MyOperator<OT,ST>(vectorSpace,dim,colptr,nnz,rowind,cvals) );
154 RCP<Thyra::LinearOpBase<ST> >
155 A = Teuchos::rcp(
new Thyra::TpetraLinearOp<OT,ST>(tpetra_A) );
161 int maxIterations = globalDim;
162 int maxRestarts = 15;
163 int gmresKrylovLength = 50;
164 int outputFrequency = 100;
165 bool outputMaxResOnly =
true;
168 Teuchos::RCP<Teuchos::ParameterList>
169 belosLOWSFPL = Teuchos::rcp(
new Teuchos::ParameterList() );
171 belosLOWSFPL->set(
"Solver Type",
"Block GMRES");
173 Teuchos::ParameterList& belosLOWSFPL_solver =
174 belosLOWSFPL->sublist(
"Solver Types");
176 Teuchos::ParameterList& belosLOWSFPL_gmres =
177 belosLOWSFPL_solver.sublist(
"Block GMRES");
179 belosLOWSFPL_gmres.set(
"Maximum Iterations",
int(maxIterations));
180 belosLOWSFPL_gmres.set(
"Convergence Tolerance",MT(maxResid));
181 belosLOWSFPL_gmres.set(
"Maximum Restarts",
int(maxRestarts));
182 belosLOWSFPL_gmres.set(
"Block Size",
int(blockSize));
183 belosLOWSFPL_gmres.set(
"Num Blocks",
int(gmresKrylovLength));
184 belosLOWSFPL_gmres.set(
"Output Frequency",
int(outputFrequency));
185 belosLOWSFPL_gmres.set(
"Show Maximum Residual Norm Only",
bool(outputMaxResOnly));
195 Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
198 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
202 belosLOWSFactory->setParameterList( belosLOWSFPL );
206 belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
209 Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
210 nsA = belosLOWSFactory->createOp();
213 Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
216 Teuchos::RCP< Thyra::MultiVectorBase<ST> >
217 b = Thyra::createMembers(domain, numRhs);
220 Teuchos::RCP< Thyra::MultiVectorBase<ST> >
221 x = Thyra::createMembers(domain, numRhs);
222 Thyra::assign(&*x, one);
225 A->apply( Thyra::NONCONJ_ELE, *x, &*b );
226 Thyra::assign(&*x, zero);
230 Thyra::SolveStatus<ST> solveStatus;
231 solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
234 *out <<
"\nBelos LOWS Status: "<< solveStatus << std::endl;
239 std::vector<MT> norm_b(numRhs), norm_res(numRhs);
240 Teuchos::RCP< Thyra::MultiVectorBase<ST> >
241 y = Thyra::createMembers(domain, numRhs);
244 Thyra::norms_2( *b, &norm_b[0] );
247 A->apply( Thyra::NONCONJ_ELE, *x, &*y );
250 Thyra::update( -one, *b, &*y );
253 Thyra::norms_2( *y, &norm_res[0] );
257 *out <<
"Final relative residual norms" << std::endl;
258 for (
int i=0; i<numRhs; ++i) {
259 rel_res = norm_res[i]/norm_b[i];
260 if (rel_res > maxResid)
262 *out <<
"RHS " << i+1 <<
" : " 263 << std::setw(16) << std::right << rel_res << std::endl;
266 return ( success ? 0 : 1 );
int main(int argc, char *argv[])
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
Simple example of a user's defined Tpetra::Operator class.