89main_(
int argc,
char *argv[], Teuchos::CommandLineProcessor& clp) {
91 typedef map_type::local_ordinal_type LocalOrdinal;
92 typedef map_type::global_ordinal_type GlobalOrdinal;
93 typedef map_type::node_type Node;
96 using Teuchos::ParameterList;
97 using Teuchos::TimeMonitor;
98 #include <Xpetra_UseShortNames.hpp>
100 bool success =
false;
108 const auto comm = Teuchos::DefaultComm<int>::getComm ();
114 Galeri::Xpetra::Parameters<GlobalOrdinal> galeriParameters(clp, 100, 100, 100,
"Laplace2D");
116 Xpetra::Parameters xpetraParameters(clp);
119 std::string xmlFileName =
"stratimikos_ParameterList.xml"; clp.setOption(
"xml", &xmlFileName,
"read parameters from an xml file");
120 std::string yamlFileName =
""; clp.setOption(
"yaml", &yamlFileName,
"read parameters from a yaml file");
121 bool printTimings =
false; clp.setOption(
"timings",
"notimings", &printTimings,
"print timings to screen");
122 bool use_stacked_timer =
false; clp.setOption(
"stacked-timer",
"no-stacked-timer", &use_stacked_timer,
"Run with or without stacked timer output");
123 std::string timingsFormat =
"table-fixed"; clp.setOption(
"time-format", &timingsFormat,
"timings format (table-fixed | table-scientific | yaml)");
124 int numVectors = 1; clp.setOption(
"multivector", &numVectors,
"number of rhs to solve simultaneously");
125 int numSolves = 1; clp.setOption(
"numSolves", &numSolves,
"number of times the system should be solved");
127 switch (clp.parse(argc,argv)) {
128 case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:
return EXIT_SUCCESS;
129 case Teuchos::CommandLineProcessor::PARSE_ERROR:
130 case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION:
return EXIT_FAILURE;
131 case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:
break;
134 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
135 Teuchos::FancyOStream& out = *fancy;
136 out.setOutputToRootOnly(0);
139 Teuchos::RCP<Teuchos::StackedTimer> stacked_timer;
140 if (use_stacked_timer)
141 stacked_timer = rcp(
new Teuchos::StackedTimer(
"Main"));
142 TimeMonitor::setStackedTimer(stacked_timer);
145 TEUCHOS_TEST_FOR_EXCEPTION(xmlFileName ==
"" && yamlFileName ==
"", std::runtime_error,
146 "Need to provide xml or yaml input file");
147 RCP<ParameterList> paramList = rcp(
new ParameterList(
"params"));
148 if (yamlFileName !=
"")
149 Teuchos::updateParametersFromYamlFileAndBroadcast(yamlFileName, paramList.ptr(), *comm);
151 Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, paramList.ptr(), *comm);
159 RCP<MultiVector> X, B;
161 std::ostringstream galeriStream;
162 Teuchos::ParameterList galeriList = galeriParameters.GetParameterList();
163 galeriStream <<
"========================================================\n" << xpetraParameters;
164 galeriStream << galeriParameters;
175 std::string matrixType = galeriParameters.GetMatrixType();
178 if (matrixType ==
"Laplace1D") {
179 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(),
"Cartesian1D", comm, galeriList);
181 }
else if (matrixType ==
"Laplace2D" || matrixType ==
"Star2D" ||
182 matrixType ==
"BigStar2D" || matrixType ==
"AnisotropicDiffusion" || matrixType ==
"Elasticity2D") {
183 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(),
"Cartesian2D", comm, galeriList);
185 }
else if (matrixType ==
"Laplace3D" || matrixType ==
"Brick3D" || matrixType ==
"Elasticity3D") {
186 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(),
"Cartesian3D", comm, galeriList);
190 if (matrixType ==
"Elasticity2D")
191 map = Xpetra::MapFactory<LO,GO,Node>::Build(map, 2);
192 if (matrixType ==
"Elasticity3D")
193 map = Xpetra::MapFactory<LO,GO,Node>::Build(map, 3);
195 galeriStream <<
"Processor subdomains in x direction: " << galeriList.get<
GO>(
"mx") << std::endl
196 <<
"Processor subdomains in y direction: " << galeriList.get<
GO>(
"my") << std::endl
197 <<
"Processor subdomains in z direction: " << galeriList.get<
GO>(
"mz") << std::endl
198 <<
"========================================================" << std::endl;
200 if (matrixType ==
"Elasticity2D" || matrixType ==
"Elasticity3D") {
202 galeriList.set(
"right boundary" ,
"Neumann");
203 galeriList.set(
"bottom boundary",
"Neumann");
204 galeriList.set(
"top boundary" ,
"Neumann");
205 galeriList.set(
"front boundary" ,
"Neumann");
206 galeriList.set(
"back boundary" ,
"Neumann");
209 RCP<Galeri::Xpetra::Problem<Map,CrsMatrixWrap,MultiVector> > Pr =
210 Galeri::Xpetra::BuildProblem<SC,LO,GO,Map,CrsMatrixWrap,MultiVector>(galeriParameters.GetMatrixType(), map, galeriList);
211 A = Pr->BuildMatrix();
213 if (matrixType ==
"Elasticity2D" ||
214 matrixType ==
"Elasticity3D") {
215 A->SetFixedBlockSize((galeriParameters.GetMatrixType() ==
"Elasticity2D") ? 2 : 3);
218 out << galeriStream.str();
219 X = MultiVectorFactory::Build(map, numVectors);
220 B = MultiVectorFactory::Build(map, numVectors);
228 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xpCrsA = Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(A);
230 RCP<const Thyra::LinearOpBase<Scalar> > thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpCrsA->getCrsMatrix());
231 RCP< Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(X));
232 RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(B);
239 Stratimikos::LinearSolverBuilder<Scalar> linearSolverBuilder;
241 typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
242 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Impl;
243 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(),
"Ifpack2");
246 linearSolverBuilder.setParameterList(paramList);
249 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
250 auto precFactory = solverFactory->getPreconditionerFactory();
251 RCP<Thyra::PreconditionerBase<Scalar> > prec;
252 Teuchos::RCP<Thyra::LinearOpWithSolveBase<Scalar> > thyraInverseA;
253 if (!precFactory.is_null()) {
254 prec = precFactory->createPrec();
257 Thyra::initializePrec<Scalar>(*precFactory, thyraA, prec.ptr());
258 thyraInverseA = solverFactory->createOp();
259 Thyra::initializePreconditionedOp<Scalar>(*solverFactory, thyraA, prec, thyraInverseA.ptr());
261 thyraInverseA = Thyra::linearOpWithSolve(*solverFactory, thyraA);
265 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr());
267 success = (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED);
269 for (
int solveno = 1; solveno < numSolves; solveno++) {
270 if (!precFactory.is_null())
271 Thyra::initializePrec<Scalar>(*precFactory, thyraA, prec.ptr());
274 status = Thyra::solve<Scalar>(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr());
276 success = success && (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED);
281 if (use_stacked_timer) {
282 stacked_timer->stop(
"Main");
283 Teuchos::StackedTimer::OutputOptions options;
284 options.output_fraction = options.output_histogram = options.output_minmax =
true;
285 stacked_timer->report(out, comm, options);
287 RCP<ParameterList> reportParams = rcp(
new ParameterList);
288 if (timingsFormat ==
"yaml") {
289 reportParams->set(
"Report format",
"YAML");
290 reportParams->set(
"YAML style",
"compact");
292 reportParams->set(
"How to merge timer sets",
"Union");
293 reportParams->set(
"alwaysWriteLocal",
false);
294 reportParams->set(
"writeGlobalStats",
true);
295 reportParams->set(
"writeZeroTimers",
false);
297 const std::string filter =
"";
299 std::ios_base::fmtflags ff(out.flags());
300 if (timingsFormat ==
"table-fixed") out << std::fixed;
301 else out << std::scientific;
302 TimeMonitor::report(comm.ptr(), out, filter, reportParams);
303 out << std::setiosflags(ff);
307 TimeMonitor::clearCounters();
311 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
313 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
326 Teuchos::CommandLineProcessor clp(
false);
328 std::vector<const char*> availableScalarTypeStrings;
329 std::vector<scalarType> availableScalarTypes;
330#ifdef HAVE_TPETRA_INST_DOUBLE
331 availableScalarTypeStrings.push_back(
"double");
332 availableScalarTypes.push_back(
DOUBLE);
334#ifdef HAVE_TPETRA_INST_FLOAT
335 availableScalarTypeStrings.push_back(
"float");
336 availableScalarTypes.push_back(
FLOAT);
338#ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
339 availableScalarTypeStrings.push_back(
"complex<double>");
342#ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
343 availableScalarTypeStrings.push_back(
"complex<float>");
346 clp.setOption(
"scalarType", &scalar, availableScalarTypes.size(), availableScalarTypes.data(), availableScalarTypeStrings.data(),
"scalar type");
347 clp.recogniseAllOptions(
false);
348 switch (clp.parse(argc, argv, NULL)) {
349 case Teuchos::CommandLineProcessor::PARSE_ERROR:
return EXIT_FAILURE;
350 case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION:
351 case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:
352 case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:
break;
355 Teuchos::GlobalMPISession session (&argc, &argv, NULL);
357#ifdef HAVE_TPETRA_INST_DOUBLE
359 return main_<double>(argc, argv, clp);
361#ifdef HAVE_TPETRA_INST_FLOAT
363 return main_<float>(argc, argv, clp);
365#ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
367 return main_<std::complex<double> >(argc, argv, clp);
369#ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
371 return main_<std::complex<float> >(argc, argv, clp);