Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
test_single_amesos_thyra_solver.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
43
44#ifndef SUN_CXX
45
47#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
48#include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
49#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
50#include "Thyra_DefaultInverseLinearOp.hpp"
51#include "Thyra_EpetraLinearOp.hpp"
52#include "Thyra_LinearOpTester.hpp"
53#include "Thyra_LinearOpWithSolveTester.hpp"
54#include "EpetraExt_readEpetraLinearSystem.h"
55#ifdef EPETRA_MPI
56#include "Epetra_MpiComm.h"
57#else
58#include "Epetra_SerialComm.h"
59#endif
60
61#endif // SUN_CXX
62
64 const std::string matrixFile
65 ,Teuchos::ParameterList *amesosLOWSFPL
66 ,const bool testTranspose_in
67 ,const int numRandomVectors
68 ,const double maxFwdError
69 ,const double maxError
70 ,const double maxResid
71 ,const bool showAllTests
72 ,const bool dumpAll
73 ,Teuchos::FancyOStream *out_arg
74 )
75{
76
77 using Teuchos::RCP;
78 using Teuchos::OSTab;
79
80 bool result, success = true;
81
82 RCP<Teuchos::FancyOStream>
83 out = Teuchos::rcp(out_arg,false);
84
85#ifndef SUN_CXX
86
87 if(out.get()) {
88 *out
89 << "\n***"
90 << "\n*** Testing Thyra::AmesosLinearOpWithSolveFactory (and Thyra::AmesosLinearOpWithSolve)"
91 << "\n***\n"
92 << "\nEchoing input options:"
93 << "\n matrixFile = " << matrixFile
94 << "\n testTranspose = " << testTranspose_in
95 << "\n numRandomVectors = " << numRandomVectors
96 << "\n maxFwdError = " << maxFwdError
97 << "\n maxError = " << maxError
98 << "\n maxResid = " << maxResid
99 << "\n showAllTests = " << showAllTests
100 << "\n dumpAll = " << dumpAll
101 << std::endl;
102 if(amesosLOWSFPL) {
103 OSTab tab(out);
104 *out
105 << "amesosLOWSFPL:\n";
106 amesosLOWSFPL->print(OSTab(out).o(),0,true);
107 }
108 }
109
110 if(out.get()) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
111
112#ifdef EPETRA_MPI
113 Epetra_MpiComm comm(MPI_COMM_WORLD);
114#else
115 Epetra_SerialComm comm;
116#endif
117 RCP<Epetra_CrsMatrix> epetra_A;
118 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
119
120 RCP<const LinearOpBase<double> >
121 A = Thyra::epetraLinearOp(epetra_A);
122
123 if(out.get() && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
124
125 if(out.get()) *out << "\nB) Creating a AmesosLinearOpWithSolveFactory object lowsFactory ...\n";
126
127 RCP<LinearOpWithSolveFactoryBase<double> >
128 lowsFactory = Teuchos::rcp(new AmesosLinearOpWithSolveFactory());
129
130 lowsFactory->setParameterList(Teuchos::rcp(amesosLOWSFPL,false));
131
132 if(out.get()) {
133 *out << "\nlowsFactory.description() = " << lowsFactory->description() << std::endl;
134 *out << "\nlowsFactory.getValidParameters() =\n";
135 lowsFactory->getValidParameters()->print(OSTab(out).o(),0,true,false);
136 }
137
138 if(out.get()) *out << "\nC) Creating a AmesosLinearOpWithSolve object nsA ...\n";
139
140 RCP<LinearOpWithSolveBase<double> >
141 nsA = lowsFactory->createOp();
142
143 initializeOp<double>(*lowsFactory, A, nsA.ptr());
144
145 bool testTranspose = testTranspose_in;
146 if (testTranspose_in && !Thyra::solveSupports(*nsA, Thyra::CONJTRANS)) {
147 if(out.get()) *out
148 << "\nChanging testTranspose=false since "
149 << nsA->description() << " does not support adjoint solves!\n";
150 testTranspose = false;
151 }
152
153 if(out.get()) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
154
155 Thyra::seed_randomize<double>(0);
156
157 LinearOpTester<double> linearOpTester;
158 linearOpTester.check_adjoint(testTranspose);
159 linearOpTester.num_random_vectors(numRandomVectors);
160 linearOpTester.set_all_error_tol(maxFwdError);
161 linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
162 linearOpTester.show_all_tests(showAllTests);
163 linearOpTester.dump_all(dumpAll);
164 result = linearOpTester.check(*nsA,out());
165 if(!result) success = false;
166
167 if(out.get()) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
168
169 LinearOpWithSolveTester<double> linearOpWithSolveTester;
170 linearOpWithSolveTester.turn_off_all_tests();
171 linearOpWithSolveTester.check_forward_default(true);
172 linearOpWithSolveTester.forward_default_residual_error_tol(1.1*maxResid);
173 linearOpWithSolveTester.forward_default_residual_warning_tol(2.0*maxResid);
174 linearOpWithSolveTester.check_forward_residual(true);
175 linearOpWithSolveTester.forward_residual_solve_tol(maxResid);
176 linearOpWithSolveTester.forward_residual_slack_error_tol(1e-1*maxResid);
177 linearOpWithSolveTester.forward_residual_slack_warning_tol(maxResid);
178 linearOpWithSolveTester.check_adjoint_default(testTranspose);
179 linearOpWithSolveTester.adjoint_default_residual_error_tol(1.1*maxResid);
180 linearOpWithSolveTester.adjoint_default_residual_warning_tol(2.0*maxResid);
181 linearOpWithSolveTester.check_adjoint_residual(testTranspose);
182 linearOpWithSolveTester.adjoint_residual_solve_tol(maxResid);
183 linearOpWithSolveTester.adjoint_residual_slack_error_tol(1e-1*maxResid);
184 linearOpWithSolveTester.adjoint_residual_slack_warning_tol(maxResid);
185 linearOpWithSolveTester.num_random_vectors(numRandomVectors);
186 linearOpWithSolveTester.show_all_tests(showAllTests);
187 linearOpWithSolveTester.dump_all(dumpAll);
188
189 LinearOpWithSolveTester<double> adjLinearOpWithSolveTester(linearOpWithSolveTester);
190 adjLinearOpWithSolveTester.check_forward_default(testTranspose);
191 adjLinearOpWithSolveTester.check_forward_residual(testTranspose);
192
193 result = linearOpWithSolveTester.check(*nsA,out.get());
194 if(!result) success = false;
195
196 if(out.get()) *out << "\nF) Uninitialize the matrix object nsA, scale the epetra_A object by 2.5, and then refactor nsA with epetra_A ...\n";
197
198 uninitializeOp<double>(*lowsFactory, nsA.ptr()); // Optional call but a good idea if changing the operator
199 epetra_A->Scale(2.5);
200 initializeOp<double>(*lowsFactory, A, nsA.ptr());
201
202 if(out.get()) *out << "\nG) Testing the LinearOpBase interface of nsA ...\n";
203
204 Thyra::seed_randomize<double>(0);
205
206 result = linearOpTester.check(*nsA, out());
207 if(!result) success = false;
208
209 if(out.get()) *out << "\nH) Testing the LinearOpWithSolveBase interface of nsA ...\n";
210
211 result = linearOpWithSolveTester.check(*nsA,out.get());
212 if(!result) success = false;
213
214 if(out.get()) *out << "\nI) Uninitialize the matrix object nsA, create a scaled (by 2.5) copy epetra_A2 of epetra_A, and then refactor nsA with epetra_A2 ...\n";
215
216 RCP<Epetra_CrsMatrix>
217 epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
218 epetra_A2->Scale(2.5);
219 RCP<const LinearOpBase<double> >
220 A2 = Thyra::epetraLinearOp(epetra_A2);
221 initializeOp<double>(*lowsFactory, A2, nsA.ptr());
222
223 if(out.get()) *out << "\nJ) Testing the LinearOpBase interface of nsA ...\n";
224
225 Thyra::seed_randomize<double>(0);
226
227 result = linearOpTester.check(*nsA, out());
228 if(!result) success = false;
229
230 if(out.get()) *out << "\nK) Testing the LinearOpWithSolveBase interface of nsA ...\n";
231
232 result = linearOpWithSolveTester.check(*nsA, out.get());
233 if(!result) success = false;
234
235 if(out.get()) *out << "\nL) Create an implicitly scaled (by 2.5) and transposed matrix A3 = scale(2.5,transpose(A)) and initialize nsA2 ...\n";
236
237 RCP<const LinearOpBase<double> >
238 A3 = scale<double>(2.5,Thyra::transpose<double>(A));
239 RCP<LinearOpWithSolveBase<double> >
240 nsA2 = linearOpWithSolve(*lowsFactory,A3);
241
242 if(out.get()) *out << "\nM) Testing the LinearOpBase interface of nsA2 ...\n";
243
244 Thyra::seed_randomize<double>(0);
245
246 result = linearOpTester.check(*nsA2,out());
247 if(!result) success = false;
248
249 if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
250
251 result = adjLinearOpWithSolveTester.check(*nsA2,out.get());
252 if(!result) success = false;
253
254 if(out.get()) *out << "\nO) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
255
256 result = linearOpTester.compare(
257 *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
258 ,out()
259 );
260 if(!result) success = false;
261
262 if(out.get()) *out << "\nP) Running example use cases ...\n";
263
264 nonExternallyPreconditionedLinearSolveUseCases(
265 *A,*lowsFactory,true,*out
266 );
267
268 if(out.get()) *out << "\nQ) Creating a DefaultInverseLinearOp object from nsA and testing the LinearOpBase interface ...\n";
269
270 RCP<const LinearOpBase<double> >
271 invA = inverse<double>(nsA.getConst());
272
273 result = linearOpTester.check(*invA,out());
274 if(!result) success = false;
275
276#else // SUN_CXX
277
278 if(out.get()) *out << "\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
279 success = false;
280
281#endif // SUN_CXX
282
283 return success;
284
285}
Concrete LinearOpWithSolveFactoryBase adapter subclass that uses Amesos direct solvers.
bool test_single_amesos_thyra_solver(const std::string matrixFile, Teuchos::ParameterList *amesosLOWSFPL, const bool testTranspose, const int numRandomVectors, const double maxFwdError, const double maxError, const double maxResid, const bool showAllTests, const bool dumpAll, Teuchos::FancyOStream *out)
Testing function for a single amesos solver with a single matrix.