Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
test_tpetra_belos.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2006) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Teuchos_ParameterList.hpp"
43#include "Teuchos_VerboseObject.hpp"
44#include "Teuchos_CommandLineProcessor.hpp"
45#include "Teuchos_XMLParameterListHelpers.hpp"
46
47#include "Tpetra_Core.hpp"
48#include "Tpetra_CrsMatrix.hpp"
49#include "MatrixMarket_Tpetra.hpp"
50
52#include "Thyra_TpetraThyraWrappers.hpp"
53#include "Thyra_LinearOpBase.hpp"
54#include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
55#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
56
58
59int main(int argc, char* argv[])
60{
61
62 using Teuchos::rcp;
63 using Teuchos::RCP;
64 using Teuchos::OSTab;
65 using Teuchos::ParameterList;
66 using Teuchos::getParameter;
67 typedef double Scalar;
68 typedef Tpetra::Map<>::local_ordinal_type LO;
69 typedef Tpetra::Map<>::global_ordinal_type GO;
70 typedef Tpetra::Map<>::node_type NO;
71
72 using Teuchos::CommandLineProcessor;
73
74 bool success = true;
75 bool verbose = true;
76
77
78 Teuchos::RCP<Teuchos::FancyOStream>
79 out = Teuchos::VerboseObjectBase::getDefaultOStream();
80
81 try {
82 Tpetra::ScopeGuard tpetraScope(&argc,&argv);
83 {
84
85 //
86 // Read options from command-line
87 //
88
89 std::string inputFile = "";
90 std::string extraParams = "";
91 bool printUnusedParams = false;
92 RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
93 int MyPID = rank(*comm);
94 verbose = ( out && (MyPID==0) );
95
96 CommandLineProcessor clp(false); // Don't throw exceptions
97 clp.addOutputSetupOptions(true);
98 clp.setOption( "input-file", &inputFile, "Input file [Required].", true );
99 clp.setOption( "extra-params", &extraParams, "Extra parameters overriding the parameters read in from --input-file");
100 clp.setOption( "print-unused-params", "no-print-unused-params", &printUnusedParams, "Whether to print all unused parameters at the end.");
101 clp.setDocString(
102 "Testing program for Trilinos Tpetra linear solvers access through Thyra."
103 );
104
105 CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
106 if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
107
108 RCP<ParameterList> paramList = rcp(new Teuchos::ParameterList());
109 if(verbose) *out << "\nReading parameters from XML file \""<<inputFile<<"\" ...\n";
110 Teuchos::updateParametersFromXmlFile(inputFile, Teuchos::inOutArg(*paramList));
111 if(extraParams.length()) {
112 if(verbose) *out << "\nAppending extra parameters from the XML string \""<<extraParams<<"\" ...\n";
113 Teuchos::updateParametersFromXmlString(extraParams, Teuchos::inOutArg(*paramList));
114 }
115
116 if(verbose) {
117 *out << "\nEchoing input parameters ...\n";
118 paramList->print(*out,1,true,false);
119 }
120
121 // Create list of valid parameter sublists
122 Teuchos::ParameterList validParamList("test_tpetra_stratimikos_solver");
123 validParamList.set("Matrix File","fileName");
124 validParamList.sublist("Linear Solver Builder").disableRecursiveValidation();
125
126 if(verbose) *out << "\nValidating top-level input parameters ...\n";
127 paramList->validateParametersAndSetDefaults(validParamList);
128
129 const std::string
130 &matrixFile = getParameter<std::string>(*paramList,"Matrix File");
131 RCP<ParameterList>
132 solverBuilderSL = sublist(paramList,"Linear Solver Builder",true);
133
134 if(verbose) *out << "\nReading in an tpetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
135
136 RCP<const Tpetra::CrsMatrix<Scalar, LO, GO, NO>> tpetra_A = Tpetra::MatrixMarket::Reader<Tpetra::CrsMatrix<Scalar, LO, GO, NO>>::readSparseFile( matrixFile, comm );
137
138 RCP<Tpetra::Vector<Scalar,LO,GO,NO>> tpetra_x0 = rcp(new Tpetra::Vector<Scalar,LO,GO,NO>(tpetra_A->getDomainMap()));
139 RCP<Tpetra::Vector<Scalar,LO,GO,NO>> tpetra_b = rcp(new Tpetra::Vector<Scalar,LO,GO,NO>(tpetra_A->getRangeMap()));
140 tpetra_x0->putScalar(1.0);
141 tpetra_A->apply(*tpetra_x0, *tpetra_b);
142 tpetra_x0->putScalar(0.0);
143
144 RCP<const Thyra::LinearOpBase<double> >
145 A = Thyra::createConstLinearOp(Teuchos::rcp_dynamic_cast<const Tpetra::Operator<Scalar,LO,GO,NO>>(tpetra_A));
146
147 RCP<Thyra::VectorBase<double> >
148 x0 = Thyra::createVector( tpetra_x0);
149 RCP<const Thyra::VectorBase<double> >
150 b = Thyra::createVector( tpetra_b);
151
152 if(verbose) *out << "\nCreating a Stratimikos::DefaultLinearSolverBuilder object ...\n";
153
154 // This is the Stratimikos main class (= factory of solver factory).
155 Stratimikos::LinearSolverBuilder<Scalar> linearSolverBuilder;
156
157 // Register Belos+Tpetra as preconditioner:
158 Stratimikos::enableBelosPrecTpetra<Tpetra::CrsMatrix<Scalar,LO,GO,NO>>(linearSolverBuilder);
159
160 linearSolverBuilder.setParameterList(solverBuilderSL);
161
162 if(verbose) *out << "\nCreating the LinearOpWithSolveFactoryBase object lowsFactory ...\n";
163 RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
164 lowsFactory = createLinearSolveStrategy(linearSolverBuilder);
165
166 if(verbose) *out << "\nChecking the LOWSB interface ...\n";
167 RCP<Thyra::LinearOpWithSolveBase<Scalar> >
168 lowsA = Thyra::linearOpWithSolve<Scalar>(*lowsFactory, A);
169
170 // Solve the linear system
171 Thyra::SolveStatus<double> status;
172 status = Thyra::solve<double>(*lowsA, Thyra::NOTRANS, *b, x0.ptr());
173
174 success = (status.solveStatus == Thyra::SOLVE_STATUS_CONVERGED ? true : false);
175
176 if(verbose && printUnusedParams) {
177 *out << "\nPrinting the parameter list (showing what was used) ...\n";
178 paramList->print(*out,1,true,true);
179 }
180
181 } //End Tpetra scope
182 }
183 catch( const std::exception &excpt ) {
184 std::cerr << "*** Caught standard exception : " << excpt.what() << std::endl;
185 success = false;
186 }
187
188 if (verbose) {
189 if(success) *out << "\nCongratulations! The test passed!\n";
190 else *out << "\nOh no! At least one of the solves failed!\n";
191 }
192
193 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
194}
map_type::global_ordinal_type GO
map_type::local_ordinal_type LO
int main(int argc, char *argv[])