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 #include "Teuchos_ParameterList.hpp" 64 #include "Teuchos_VerboseObject.hpp" 65 #include "Teuchos_GlobalMPISession.hpp" 68 # include "Tpetra_MpiPlatform.hpp" 70 # include "Tpetra_SerialPlatform.hpp" 73 int main(
int argc,
char* argv[])
78 Teuchos::RCP<Teuchos::FancyOStream>
79 out = Teuchos::VerboseObjectBase::getDefaultOStream();
81 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
90 typedef std::complex<double>
ST;
92 typedef std::complex<double>
ST;
96 typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;
106 MPI_Comm mpiComm = MPI_COMM_WORLD;
107 const Tpetra::MpiPlatform<OT,OT> ordinalPlatform(mpiComm);
108 const Tpetra::MpiPlatform<OT,ST> scalarPlatform(mpiComm);
110 const Tpetra::SerialPlatform<OT,OT> ordinalPlatform;
111 const Tpetra::SerialPlatform<OT,ST> scalarPlatform;
118 const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
119 const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
122 RCP<Tpetra::CisMatrix<OT,ST> >
123 tpetra_A = rcp(
new Tpetra::CisMatrix<OT,ST>(vectorSpace));
126 const int numMyElements = vectorSpace.getNumMyEntries();
127 const std::vector<int> &myGlobalElements = vectorSpace.elementSpace().getMyGlobalElements();
130 const ST offDiag = -1.0, diag = 2.0;
131 int numEntries;
ST values[3];
int indexes[3];
132 for(
int k = 0; k < numMyElements; ++k ) {
133 const int rowIndex = myGlobalElements[k];
134 if( rowIndex == 0 ) {
136 values[0] = diag; values[1] = offDiag;
137 indexes[0] = 0; indexes[1] = 1;
139 else if( rowIndex == globalDim - 1 ) {
141 values[0] = offDiag; values[1] = diag;
142 indexes[0] = globalDim-2; indexes[1] = globalDim-1;
146 values[0] = offDiag; values[1] = diag; values[2] = offDiag;
147 indexes[0] = rowIndex-1; indexes[1] = rowIndex; indexes[2] = rowIndex+1;
149 tpetra_A->submitEntries(Tpetra::Insert,rowIndex,numEntries,values,indexes);
153 tpetra_A->fillComplete();
156 RCP<Thyra::LinearOpBase<ST> >
157 A = Teuchos::rcp(
new Thyra::TpetraLinearOp<OT,ST>(
158 Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_A) ) );
164 int maxIterations = globalDim;
166 int gmresKrylovLength = globalDim;
167 int outputFrequency = 100;
168 bool outputMaxResOnly =
true;
171 ST one = Teuchos::ScalarTraits<ST>::one();
172 ST zero = Teuchos::ScalarTraits<ST>::zero();
174 Teuchos::RCP<Teuchos::ParameterList>
175 belosLOWSFPL = Teuchos::rcp(
new Teuchos::ParameterList() );
177 belosLOWSFPL->set(
"Solver Type",
"Block GMRES");
179 Teuchos::ParameterList& belosLOWSFPL_solver =
180 belosLOWSFPL->sublist(
"Solver Types");
182 Teuchos::ParameterList& belosLOWSFPL_gmres =
183 belosLOWSFPL_solver.sublist(
"Block GMRES");
185 belosLOWSFPL_gmres.set(
"Maximum Iterations",
int(maxIterations));
186 belosLOWSFPL_gmres.set(
"Convergence Tolerance",
double(maxResid));
187 belosLOWSFPL_gmres.set(
"Maximum Restarts",
int(maxRestarts));
188 belosLOWSFPL_gmres.set(
"Block Size",
int(blockSize));
189 belosLOWSFPL_gmres.set(
"Num Blocks",
int(gmresKrylovLength));
190 belosLOWSFPL_gmres.set(
"Output Frequency",
int(outputFrequency));
191 belosLOWSFPL_gmres.set(
"Show Maximum Residual Norm Only",
bool(outputMaxResOnly));
207 Teuchos::RCP<Thyra::MultiVectorBase<ST> > x, b;
215 Teuchos::RCP<Tpetra::Vector<OT,ST> > tpetra_b =
216 Teuchos::rcp(
new Tpetra::Vector<OT,ST>(vectorSpace) );
219 tpetra_b->setAllToRandom();
222 b = Thyra::create_Vector(tpetra_b);
225 Teuchos::RCP<Tpetra::Vector<OT,ST> > tpetra_x =
226 Teuchos::rcp(
new Tpetra::Vector<OT,ST>(vectorSpace) );
229 tpetra_x->setAllToScalar( zero );
232 x = Thyra::create_Vector(tpetra_x);
243 Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
246 x = Thyra::createMembers(domain, numRhs);
247 b = Thyra::createMembers(domain, numRhs);
252 for (
int j=0; j<numRhs; ++j ) {
256 Teuchos::RCP<Tpetra::Vector<OT,ST> >
257 tpetra_b_j = Thyra::get_Tpetra_Vector(vectorSpace,b->col(j));
258 tpetra_b_j->setAllToRandom();
262 Teuchos::RCP<Tpetra::Vector<OT,ST> >
263 tpetra_x_j = Thyra::get_Tpetra_Vector(vectorSpace,x->col(j));
264 tpetra_x_j->setAllToScalar( zero );
269 for (
int i=0; i<numMyElements; ++i ) {
273 const int rowIndex = myGlobalElements[ i ];
276 (*tpetra_x_j)[ rowIndex ] = zero;
291 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
296 belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
299 belosLOWSFactory->setParameterList( belosLOWSFPL );
302 Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
303 nsA = belosLOWSFactory->createOp();
306 Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
310 Thyra::SolveStatus<ST> solveStatus;
311 solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
314 *out <<
"\nBelos LOWS Status: "<< solveStatus << std::endl;
323 std::vector<ST> norm_b(numRhs), norm_res(numRhs);
325 for (
int j=0; j<numRhs; ++j ) {
329 Teuchos::RCP<Tpetra::Vector<OT,ST> >
330 tpetra_x_j = Thyra::get_Tpetra_Vector(vectorSpace,x->col(j));
332 Teuchos::RCP<Tpetra::Vector<OT,ST> >
333 tpetra_b_j = Thyra::get_Tpetra_Vector(vectorSpace,b->col(j));
336 norm_b[j] = tpetra_b_j->norm2();
339 Tpetra::Vector<OT,ST> y(vectorSpace);
340 tpetra_A->apply( *tpetra_x_j, y );
343 y.update( one, *tpetra_b_j, -one );
346 norm_res[j] = y.norm2();
351 *out <<
"Final relative residual norms" << std::endl;
352 for (
int i=0; i<numRhs; ++i) {
353 rel_res = Teuchos::ScalarTraits<ST>::real(norm_res[i]/norm_b[i]);
354 if (rel_res > maxResid)
356 *out <<
"RHS " << i+1 <<
" : " 357 << std::setw(16) << std::right << rel_res << std::endl;
360 return ( success ? 0 : 1 );
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
int main(int argc, char *argv[])