53 #include "Thyra_BelosLinearOpWithSolveFactory.hpp" 54 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp" 55 #include "Thyra_PreconditionerFactoryHelpers.hpp" 56 #include "Thyra_TpetraLinearOp.hpp" 57 #include "Thyra_DefaultPreconditioner.hpp" 60 #include "Tpetra_ElementSpace.hpp" 61 #include "Tpetra_VectorSpace.hpp" 62 #include "Tpetra_CisMatrix.hpp" 65 #include "Teuchos_ParameterList.hpp" 66 #include "Teuchos_VerboseObject.hpp" 67 #include "Teuchos_GlobalMPISession.hpp" 70 # include "Tpetra_MpiPlatform.hpp" 72 # include "Tpetra_SerialPlatform.hpp" 75 int main(
int argc,
char* argv[])
80 Teuchos::RCP<Teuchos::FancyOStream>
81 out = Teuchos::VerboseObjectBase::getDefaultOStream();
83 Teuchos::GlobalMPISession mpiSession(&argc,&argv);
90 typedef std::complex<double>
ST;
92 typedef std::complex<double>
ST;
96 typedef Teuchos::ScalarTraits<ST>::magnitudeType MT;
99 ST one = Teuchos::ScalarTraits<ST>::one();
100 ST zero = Teuchos::ScalarTraits<ST>::zero();
107 MPI_Comm mpiComm = MPI_COMM_WORLD;
108 const Tpetra::MpiPlatform<OT,OT> ordinalPlatform(mpiComm);
109 const Tpetra::MpiPlatform<OT,ST> scalarPlatform(mpiComm);
111 const Tpetra::SerialPlatform<OT,OT> ordinalPlatform;
112 const Tpetra::SerialPlatform<OT,ST> scalarPlatform;
119 const Tpetra::ElementSpace<OT> elementSpace(globalDim,0,ordinalPlatform);
120 const Tpetra::VectorSpace<OT,ST> vectorSpace(elementSpace,scalarPlatform);
123 RCP<Tpetra::CisMatrix<OT,ST> >
124 tpetra_A = rcp(
new Tpetra::CisMatrix<OT,ST>(vectorSpace));
127 RCP<Tpetra::CisMatrix<OT,ST> >
128 tpetra_Prec = rcp(
new Tpetra::CisMatrix<OT,ST>(vectorSpace));
131 const int numMyElements = vectorSpace.getNumMyEntries();
132 const std::vector<int> &myGlobalElements = vectorSpace.elementSpace().getMyGlobalElements();
135 const ST offDiag = -1.0, diag = 2.0;
136 int numEntries;
ST values[3];
int indexes[3];
137 ST prec_value[1];
int prec_index[1];
138 prec_value[0] = one/diag;
139 for(
int k = 0; k < numMyElements; ++k ) {
140 const int rowIndex = myGlobalElements[k];
141 prec_index[0] = rowIndex;
142 if( rowIndex == 0 ) {
144 values[0] = diag; values[1] = offDiag;
145 indexes[0] = 0; indexes[1] = 1;
147 else if( rowIndex == globalDim - 1 ) {
149 values[0] = offDiag; values[1] = diag;
150 indexes[0] = globalDim-2; indexes[1] = globalDim-1;
154 values[0] = offDiag; values[1] = diag; values[2] = offDiag;
155 indexes[0] = rowIndex-1; indexes[1] = rowIndex; indexes[2] = rowIndex+1;
157 tpetra_A->submitEntries(Tpetra::Insert,rowIndex,numEntries,values,indexes);
158 tpetra_Prec->submitEntries(Tpetra::Insert,rowIndex,1,prec_value,prec_index);
162 tpetra_A->fillComplete();
163 tpetra_Prec->fillComplete();
166 RCP<Thyra::LinearOpBase<ST> >
167 A = Teuchos::rcp(
new Thyra::TpetraLinearOp<OT,ST>(
168 Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_A) ) );
171 RCP<Thyra::LinearOpBase<ST> >
172 Prec = Teuchos::rcp(
new Thyra::TpetraLinearOp<OT,ST>(
173 Teuchos::rcp_implicit_cast<Tpetra::Operator<OT,ST> >(tpetra_Prec) ) );
176 RCP<const Thyra::DefaultPreconditioner<ST> >
177 DefPrec = Thyra::unspecifiedPrec<ST>(Prec);
184 int maxIterations = globalDim;
186 int gmresKrylovLength = globalDim;
187 int outputFrequency = 100;
188 bool outputMaxResOnly =
true;
191 Teuchos::RCP<Teuchos::ParameterList>
192 belosLOWSFPL = Teuchos::rcp(
new Teuchos::ParameterList() );
194 belosLOWSFPL->set(
"Solver Type",
"Block GMRES");
196 Teuchos::ParameterList& belosLOWSFPL_solver =
197 belosLOWSFPL->sublist(
"Solver Types");
199 Teuchos::ParameterList& belosLOWSFPL_gmres =
200 belosLOWSFPL_solver.sublist(
"Block GMRES");
202 belosLOWSFPL_gmres.set(
"Maximum Iterations",
int(maxIterations));
203 belosLOWSFPL_gmres.set(
"Convergence Tolerance",
double(maxResid));
204 belosLOWSFPL_gmres.set(
"Maximum Restarts",
int(maxRestarts));
205 belosLOWSFPL_gmres.set(
"Block Size",
int(blockSize));
206 belosLOWSFPL_gmres.set(
"Num Blocks",
int(gmresKrylovLength));
207 belosLOWSFPL_gmres.set(
"Output Frequency",
int(outputFrequency));
208 belosLOWSFPL_gmres.set(
"Show Maximum Residual Norm Only",
bool(outputMaxResOnly));
218 Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domain = A->domain();
221 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<ST> >
225 belosLOWSFactory->setParameterList( belosLOWSFPL );
229 belosLOWSFactory->setVerbLevel(Teuchos::VERB_LOW);
232 Teuchos::RCP<Thyra::LinearOpWithSolveBase<ST> >
233 nsA = belosLOWSFactory->createOp();
237 Thyra::initializeOp<ST>( *belosLOWSFactory, A, &*nsA );
240 Teuchos::RCP< Thyra::MultiVectorBase<ST> >
241 b = Thyra::createMembers(domain, numRhs);
242 Thyra::seed_randomize<ST>(0);
243 Thyra::randomize(-one, one, &*b);
246 Teuchos::RCP< Thyra::MultiVectorBase<ST> >
247 x = Thyra::createMembers(domain, numRhs);
248 Thyra::assign(&*x, zero);
252 Thyra::SolveStatus<ST> solveStatus;
253 solveStatus = Thyra::solve( *nsA, Thyra::NONCONJ_ELE, *b, &*x );
256 *out <<
"\nBelos LOWS Status: "<< solveStatus << std::endl;
261 std::vector<MT> norm_b(numRhs), norm_res(numRhs);
262 Teuchos::RCP< Thyra::MultiVectorBase<ST> >
263 y = Thyra::createMembers(domain, numRhs);
266 Thyra::norms_2( *b, &norm_b[0] );
269 A->apply( Thyra::NONCONJ_ELE, *x, &*y );
272 Thyra::update( -one, *b, &*y );
275 Thyra::norms_2( *y, &norm_res[0] );
279 *out <<
"Final relative residual norms" << std::endl;
280 for (
int i=0; i<numRhs; ++i) {
281 rel_res = norm_res[i]/norm_b[i];
282 if (rel_res > maxResid)
284 *out <<
"RHS " << i+1 <<
" : " 285 << std::setw(16) << std::right << rel_res << std::endl;
288 return ( success ? 0 : 1 );
int main(int argc, char *argv[])
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.