66 #include <Thyra_LinearOpWithSolveBase.hpp>
67 #include <Thyra_VectorBase.hpp>
68 #include <Thyra_SolveSupportTypes.hpp>
71 #include <Stratimikos_LinearSolverBuilder.hpp>
72 #include <Stratimikos_InternalConfig.h>
75 #include <Xpetra_Parameters.hpp>
76 #include <Xpetra_ThyraUtils.hpp>
79 #include <Galeri_XpetraParameters.hpp>
80 #include <Galeri_XpetraProblemFactory.hpp>
81 #include <Galeri_XpetraUtils.hpp>
82 #include <Galeri_XpetraMaps.hpp>
85 template <
class Scalar>
89 typedef map_type::local_ordinal_type LocalOrdinal;
90 typedef map_type::global_ordinal_type GlobalOrdinal;
91 typedef map_type::node_type Node;
96 #include <Xpetra_UseShortNames.hpp>
112 Galeri::Xpetra::Parameters<GlobalOrdinal> galeriParameters(clp, 100, 100, 100,
"Laplace2D");
114 Xpetra::Parameters xpetraParameters(clp);
117 std::string xmlFileName =
"stratimikos_ParameterList.xml"; clp.
setOption(
"xml", &xmlFileName,
"read parameters from an xml file");
118 std::string yamlFileName =
""; clp.
setOption(
"yaml", &yamlFileName,
"read parameters from a yaml file");
119 bool printTimings =
false; clp.
setOption(
"timings",
"notimings", &printTimings,
"print timings to screen");
120 bool use_stacked_timer =
false; clp.
setOption(
"stacked-timer",
"no-stacked-timer", &use_stacked_timer,
"Run with or without stacked timer output");
121 std::string timingsFormat =
"table-fixed"; clp.
setOption(
"time-format", &timingsFormat,
"timings format (table-fixed | table-scientific | yaml)");
122 int numVectors = 1; clp.
setOption(
"multivector", &numVectors,
"number of rhs to solve simultaneously");
123 int numSolves = 1; clp.
setOption(
"numSolves", &numSolves,
"number of times the system should be solved");
125 switch (clp.
parse(argc,argv)) {
132 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
138 if (use_stacked_timer)
140 TimeMonitor::setStackedTimer(stacked_timer);
144 "Need to provide xml or yaml input file");
145 RCP<ParameterList> paramList =
rcp(
new ParameterList(
"params"));
146 if (yamlFileName !=
"")
147 Teuchos::updateParametersFromYamlFileAndBroadcast(yamlFileName, paramList.ptr(), *comm);
149 Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, paramList.ptr(), *comm);
157 RCP<MultiVector> X,
B;
159 std::ostringstream galeriStream;
161 galeriStream <<
"========================================================\n" << xpetraParameters;
162 galeriStream << galeriParameters;
173 std::string matrixType = galeriParameters.GetMatrixType();
176 if (matrixType ==
"Laplace1D") {
177 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(),
"Cartesian1D", comm, galeriList);
179 }
else if (matrixType ==
"Laplace2D" || matrixType ==
"Star2D" ||
180 matrixType ==
"BigStar2D" || matrixType ==
"AnisotropicDiffusion" || matrixType ==
"Elasticity2D") {
181 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(),
"Cartesian2D", comm, galeriList);
183 }
else if (matrixType ==
"Laplace3D" || matrixType ==
"Brick3D" || matrixType ==
"Elasticity3D") {
184 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(),
"Cartesian3D", comm, galeriList);
188 if (matrixType ==
"Elasticity2D")
189 map = Xpetra::MapFactory<LO,GO,Node>::Build(map, 2);
190 if (matrixType ==
"Elasticity3D")
191 map = Xpetra::MapFactory<LO,GO,Node>::Build(map, 3);
193 galeriStream <<
"Processor subdomains in x direction: " << galeriList.
get<
GO>(
"mx") << std::endl
194 <<
"Processor subdomains in y direction: " << galeriList.
get<
GO>(
"my") << std::endl
195 <<
"Processor subdomains in z direction: " << galeriList.
get<
GO>(
"mz") << std::endl
196 <<
"========================================================" << std::endl;
198 if (matrixType ==
"Elasticity2D" || matrixType ==
"Elasticity3D") {
200 galeriList.
set(
"right boundary" ,
"Neumann");
201 galeriList.
set(
"bottom boundary",
"Neumann");
202 galeriList.
set(
"top boundary" ,
"Neumann");
203 galeriList.
set(
"front boundary" ,
"Neumann");
204 galeriList.
set(
"back boundary" ,
"Neumann");
207 RCP<Galeri::Xpetra::Problem<Map,CrsMatrixWrap,MultiVector> > Pr =
208 Galeri::Xpetra::BuildProblem<SC,LO,GO,Map,CrsMatrixWrap,MultiVector>(galeriParameters.GetMatrixType(), map, galeriList);
209 A = Pr->BuildMatrix();
211 if (matrixType ==
"Elasticity2D" ||
212 matrixType ==
"Elasticity3D") {
213 A->SetFixedBlockSize((galeriParameters.GetMatrixType() ==
"Elasticity2D") ? 2 : 3);
216 out << galeriStream.str();
217 X = MultiVectorFactory::Build(map, numVectors);
218 B = MultiVectorFactory::Build(map, numVectors);
226 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xpCrsA = Teuchos::rcp_dynamic_cast<
const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(A);
228 RCP<const Thyra::LinearOpBase<Scalar> > thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpCrsA->getCrsMatrix());
229 RCP< Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(X));
230 RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(B);
243 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
244 auto precFactory = solverFactory->getPreconditionerFactory();
245 RCP<Thyra::PreconditionerBase<Scalar> > prec;
247 if (!precFactory.is_null()) {
248 prec = precFactory->createPrec();
251 Thyra::initializePrec<Scalar>(*precFactory, thyraA, prec.ptr());
252 thyraInverseA = solverFactory->createOp();
253 Thyra::initializePreconditionedOp<Scalar>(*solverFactory, thyraA, prec, thyraInverseA.ptr());
255 thyraInverseA = Thyra::linearOpWithSolve(*solverFactory, thyraA);
259 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.
ptr());
261 success = (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED);
263 for (
int solveno = 1; solveno < numSolves; solveno++) {
264 if (!precFactory.is_null())
265 Thyra::initializePrec<Scalar>(*precFactory, thyraA, prec.ptr());
268 status = Thyra::solve<Scalar>(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr());
270 success = success && (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED);
275 if (use_stacked_timer) {
276 stacked_timer->
stop(
"Main");
279 stacked_timer->
report(out, comm, options);
281 RCP<ParameterList> reportParams =
rcp(
new ParameterList);
282 if (timingsFormat ==
"yaml") {
283 reportParams->set(
"Report format",
"YAML");
284 reportParams->set(
"YAML style",
"compact");
286 reportParams->set(
"How to merge timer sets",
"Union");
287 reportParams->set(
"alwaysWriteLocal",
false);
288 reportParams->set(
"writeGlobalStats",
true);
289 reportParams->set(
"writeZeroTimers",
false);
291 const std::string filter =
"";
293 std::ios_base::fmtflags ff(out.flags());
294 if (timingsFormat ==
"table-fixed") out << std::fixed;
295 else out << std::scientific;
296 TimeMonitor::report(comm.ptr(), out, filter, reportParams);
297 out << std::setiosflags(ff);
301 TimeMonitor::clearCounters();
307 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
324 std::vector<const char*> availableScalarTypeStrings;
325 std::vector<scalarType> availableScalarTypes;
326 #ifdef HAVE_TPETRA_INST_DOUBLE
327 availableScalarTypeStrings.push_back(
"double");
328 availableScalarTypes.push_back(
DOUBLE);
330 #ifdef HAVE_TPETRA_INST_FLOAT
331 availableScalarTypeStrings.push_back(
"float");
332 availableScalarTypes.push_back(
FLOAT);
334 #ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
335 availableScalarTypeStrings.push_back(
"complex<double>");
338 #ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
339 availableScalarTypeStrings.push_back(
"complex<float>");
342 clp.
setOption(
"scalarType", &scalar, availableScalarTypes.size(), availableScalarTypes.data(), availableScalarTypeStrings.data(),
"scalar type");
344 switch (clp.
parse(argc, argv, NULL)) {
351 #ifdef HAVE_TPETRA_INST_DOUBLE
353 return main_<double>(argc, argv, clp);
355 #ifdef HAVE_TPETRA_INST_FLOAT
357 return main_<float>(argc, argv, clp);
359 #ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
361 return main_<std::complex<double> >(argc, argv, clp);
363 #ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
365 return main_<std::complex<float> >(argc, argv, clp);
int main(int argc, char *argv[])
void stop(const std::string &name, const bool pop_kokkos_profiling_region=true)
Tpetra::Map< int, int > map_type
void recogniseAllOptions(const bool &recogniseAllOptions)
T & get(ParameterList &l, const std::string &name)
int main_(int argc, char *argv[], Teuchos::CommandLineProcessor &clp)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setParameterList(RCP< ParameterList > const ¶mList)
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
void report(std::ostream &os)
Concrete subclass of Thyra::LinearSolverBuilderBase for creating Thyra::LinearOpWithSolveFactoryBase ...
map_type::global_ordinal_type GO