Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
galeri_xpetra_driver.cpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42#include <iostream>
43
44/*
45 Call Ifpack2 via the Stratimikos interface.
46
47Usage:
48./Ifpack2_Stratimikos.exe : use xml configuration file stratimikos_ParameterList.xml
49
50Note:
51The source code is not MueLu specific and can be used with any Stratimikos strategy.
52*/
53
54// Teuchos includes
55#include <Teuchos_ConfigDefs.hpp>
56#include <Teuchos_GlobalMPISession.hpp>
57#include <Teuchos_StackedTimer.hpp>
58#include <Teuchos_RCP.hpp>
59#include <Teuchos_XMLParameterListHelpers.hpp>
60#include <Teuchos_YamlParameterListHelpers.hpp>
61#include <Teuchos_StandardCatchMacros.hpp>
62#include "Teuchos_AbstractFactoryStd.hpp"
63
64// Thyra includes
65#include <Thyra_LinearOpWithSolveBase.hpp>
66#include <Thyra_VectorBase.hpp>
67#include <Thyra_SolveSupportTypes.hpp>
68
69// Stratimikos includes
70#include <Stratimikos_LinearSolverBuilder.hpp>
71
72// Xpetra include
73#include <Xpetra_Parameters.hpp>
74#include <Xpetra_ThyraUtils.hpp>
75
76// Galeri includes
77#include <Galeri_XpetraParameters.hpp>
78#include <Galeri_XpetraProblemFactory.hpp>
79#include <Galeri_XpetraUtils.hpp>
80#include <Galeri_XpetraMaps.hpp>
81
82// Ifpack2 includes
83#ifdef HAVE_STRATIMIKOS_IFPACK2
84# include <Thyra_Ifpack2PreconditionerFactory.hpp>
85#endif
86
87template <class Scalar>
88int
89main_(int argc, char *argv[], Teuchos::CommandLineProcessor& clp) {
90 typedef Tpetra::Map<> map_type;
91 typedef map_type::local_ordinal_type LocalOrdinal;
92 typedef map_type::global_ordinal_type GlobalOrdinal;
93 typedef map_type::node_type Node;
94 using Teuchos::RCP;
95 using Teuchos::rcp;
96 using Teuchos::ParameterList;
97 using Teuchos::TimeMonitor;
98 #include <Xpetra_UseShortNames.hpp>
99
100 bool success = false;
101 bool verbose = true;
102 try {
103
104 //
105 // MPI initialization
106 //
107 // Teuchos::CommandLineProcessor clp(false);
108 const auto comm = Teuchos::DefaultComm<int>::getComm ();
109
110 //
111 // Parameters
112 //
113 // manage parameters of the test case
114 Galeri::Xpetra::Parameters<GlobalOrdinal> galeriParameters(clp, 100, 100, 100, "Laplace2D");
115 // manage parameters of Xpetra
116 Xpetra::Parameters xpetraParameters(clp);
117
118 // command line parameters
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");
126
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;
132 }
133
134 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
135 Teuchos::FancyOStream& out = *fancy;
136 out.setOutputToRootOnly(0);
137
138 // Set up timers
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);
143
144 // Read in parameter list
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);
150 else
151 Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, paramList.ptr(), *comm);
152
153 //
154 // Construct the problem
155 //
156
157 RCP<Matrix> A;
158 RCP<const Map> map;
159 RCP<MultiVector> X, B;
160
161 std::ostringstream galeriStream;
162 Teuchos::ParameterList galeriList = galeriParameters.GetParameterList();
163 galeriStream << "========================================================\n" << xpetraParameters;
164 galeriStream << galeriParameters;
165
166 // Galeri will attempt to create a square-as-possible distribution of subdomains di, e.g.,
167 // d1 d2 d3
168 // d4 d5 d6
169 // d7 d8 d9
170 // d10 d11 d12
171 // A perfect distribution is only possible when the #processors is a perfect square.
172 // This *will* result in "strip" distribution if the #processors is a prime number or if the factors are very different in
173 // size. For example, np=14 will give a 7-by-2 distribution.
174 // If you don't want Galeri to do this, specify mx or my on the galeriList.
175 std::string matrixType = galeriParameters.GetMatrixType();
176
177 // Create map
178 if (matrixType == "Laplace1D") {
179 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(), "Cartesian1D", comm, galeriList);
180
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);
184
185 } else if (matrixType == "Laplace3D" || matrixType == "Brick3D" || matrixType == "Elasticity3D") {
186 map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(), "Cartesian3D", comm, galeriList);
187 }
188
189 // Expand map to do multiple DOF per node for block problems
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);
194
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;
199
200 if (matrixType == "Elasticity2D" || matrixType == "Elasticity3D") {
201 // Our default test case for elasticity: all boundaries of a square/cube have Neumann b.c. except left which has Dirichlet
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");
207 }
208
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();
212
213 if (matrixType == "Elasticity2D" ||
214 matrixType == "Elasticity3D") {
215 A->SetFixedBlockSize((galeriParameters.GetMatrixType() == "Elasticity2D") ? 2 : 3);
216 }
217
218 out << galeriStream.str();
219 X = MultiVectorFactory::Build(map, numVectors);
220 B = MultiVectorFactory::Build(map, numVectors);
221 B->putScalar(1);
222 X->putScalar(0);
223
224 //
225 // Build Thyra linear algebra objects
226 //
227
228 RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > xpCrsA = Teuchos::rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(A);
229
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);
233
234 //
235 // Build Stratimikos solver
236 //
237
238 // This is the Stratimikos main class (= factory of solver factory).
239 Stratimikos::LinearSolverBuilder<Scalar> linearSolverBuilder;
240 // Register Ifpack2 as a Stratimikos preconditioner strategy.
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");
244
245 // Setup solver parameters using a Stratimikos parameter list.
246 linearSolverBuilder.setParameterList(paramList);
247
248 // Build a new "solver factory" according to the previously specified parameter list.
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();
255
256 // Build a Thyra operator corresponding to A^{-1} computed using the Stratimikos solver.
257 Thyra::initializePrec<Scalar>(*precFactory, thyraA, prec.ptr());
258 thyraInverseA = solverFactory->createOp();
259 Thyra::initializePreconditionedOp<Scalar>(*solverFactory, thyraA, prec, thyraInverseA.ptr());
260 } else {
261 thyraInverseA = Thyra::linearOpWithSolve(*solverFactory, thyraA);
262 }
263
264 // Solve Ax = b.
265 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr());
266
267 success = (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED);
268
269 for (int solveno = 1; solveno < numSolves; solveno++) {
270 if (!precFactory.is_null())
271 Thyra::initializePrec<Scalar>(*precFactory, thyraA, prec.ptr());
272 thyraX->assign(0.);
273
274 status = Thyra::solve<Scalar>(*thyraInverseA, Thyra::NOTRANS, *thyraB, thyraX.ptr());
275
276 success = success && (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED);
277 }
278
279 // print timings
280 if (printTimings) {
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);
286 } else {
287 RCP<ParameterList> reportParams = rcp(new ParameterList);
288 if (timingsFormat == "yaml") {
289 reportParams->set("Report format", "YAML"); // "Table" or "YAML"
290 reportParams->set("YAML style", "compact"); // "spacious" or "compact"
291 }
292 reportParams->set("How to merge timer sets", "Union");
293 reportParams->set("alwaysWriteLocal", false);
294 reportParams->set("writeGlobalStats", true);
295 reportParams->set("writeZeroTimers", false);
296
297 const std::string filter = "";
298
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);
304 }
305 }
306
307 TimeMonitor::clearCounters();
308 out << std::endl;
309
310 }
311 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
312
313 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
314}
315
316
323
324int
325main(int argc, char *argv[]) {
326 Teuchos::CommandLineProcessor clp(false);
327 scalarType scalar = DOUBLE;
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);
333#endif
334#ifdef HAVE_TPETRA_INST_FLOAT
335 availableScalarTypeStrings.push_back("float");
336 availableScalarTypes.push_back(FLOAT);
337#endif
338#ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
339 availableScalarTypeStrings.push_back("complex<double>");
340 availableScalarTypes.push_back(COMPLEX_DOUBLE);
341#endif
342#ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
343 availableScalarTypeStrings.push_back("complex<float>");
344 availableScalarTypes.push_back(COMPLEX_FLOAT);
345#endif
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;
353 }
354
355 Teuchos::GlobalMPISession session (&argc, &argv, NULL);
356
357#ifdef HAVE_TPETRA_INST_DOUBLE
358 if (scalar == DOUBLE)
359 return main_<double>(argc, argv, clp);
360#endif
361#ifdef HAVE_TPETRA_INST_FLOAT
362 if (scalar == FLOAT)
363 return main_<float>(argc, argv, clp);
364#endif
365#ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
366 if (scalar == COMPLEX_DOUBLE)
367 return main_<std::complex<double> >(argc, argv, clp);
368#endif
369#ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
370 if (scalar == COMPLEX_FLOAT)
371 return main_<std::complex<float> >(argc, argv, clp);
372#endif
373}
int main(int argc, char *argv[])
int main_(int argc, char *argv[], Teuchos::CommandLineProcessor &clp)
Tpetra::Map< int, int > map_type
map_type::global_ordinal_type GO