Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
test_single_aztecoo_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
42
43#ifndef SUN_CXX
44
46#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
47#include "Thyra_PreconditionerFactoryHelpers.hpp"
48#include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
49#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
50#include "Thyra_EpetraLinearOp.hpp"
51#include "Thyra_LinearOpTester.hpp"
52#include "Thyra_LinearOpWithSolveBase.hpp"
53#include "Thyra_LinearOpWithSolveTester.hpp"
54#include "EpetraExt_readEpetraLinearSystem.h"
55#include "Epetra_SerialComm.h"
56#include "Teuchos_ParameterList.hpp"
57#ifdef HAVE_AZTECOO_IFPACK
59#endif
60#ifdef HAVE_MPI
61# include "Epetra_MpiComm.h"
62#endif
63
64#endif // SUN_CXX
65
67 const std::string matrixFile
68 ,const bool testTranspose
69 ,const int numRandomVectors
70 ,const double maxFwdError
71 ,const double maxResid
72 ,const double maxSolutionError
73 ,const bool showAllTests
74 ,const bool dumpAll
75 ,Teuchos::ParameterList *aztecooLOWSFPL
76 ,Teuchos::FancyOStream *out_arg
77 )
78{
79 using Teuchos::rcp;
80 using Teuchos::RCP;
81 using Teuchos::OSTab;
82 typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
83 bool result, success = true;
84
85 RCP<Teuchos::FancyOStream>
86 out = Teuchos::rcp(out_arg,false);
87
88 try {
89
90#ifndef SUN_CXX
91
92 if(out.get()) {
93 *out
94 << "\n***"
95 << "\n*** Testing Thyra::AztecOOLinearOpWithSolveFactory (and Thyra::AztecOOLinearOpWithSolve)"
96 << "\n***\n"
97 << "\nEchoing input options:"
98 << "\n matrixFile = " << matrixFile
99 << "\n testTranspose = " << testTranspose
100 << "\n numRandomVectors = " << numRandomVectors
101 << "\n maxFwdError = " << maxFwdError
102 << "\n maxResid = " << maxResid
103 << "\n showAllTests = " << showAllTests
104 << "\n dumpAll = " << dumpAll
105 << std::endl;
106 }
107
108 const bool useAztecPrec = (
109 aztecooLOWSFPL
110 &&
111 aztecooLOWSFPL->sublist("Forward Solve")
112 .sublist("AztecOO Settings")
113 .get("Aztec Preconditioner","none")!="none"
114 );
115
116 if(out.get()) {
117 if(useAztecPrec)
118 *out << "\nUsing aztec preconditioning so we will not test adjoint solves using internal preconditioning ...\n";
119 }
120
121 if(out.get()) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
122
123#ifdef HAVE_MPI
124 Epetra_MpiComm comm(MPI_COMM_WORLD);
125#else
126 Epetra_SerialComm comm;
127#endif
128 RCP<Epetra_CrsMatrix> epetra_A;
129 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
130
131 RCP<const LinearOpBase<double> > A = Thyra::epetraLinearOp(epetra_A);
132
133 if(out.get() && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
134
135 if(out.get()) *out << "\nB) Creating a AztecOOLinearOpWithSolveFactory object opFactory ...\n";
136
137 RCP<LinearOpWithSolveFactoryBase<double> >
138 lowsFactory = Teuchos::rcp(new AztecOOLinearOpWithSolveFactory());
139 if(out.get()) {
140 *out << "\nlowsFactory.getValidParameters() initially:\n";
141 lowsFactory->getValidParameters()->print(OSTab(out).o(),PLPrintOptions().showTypes(true).showDoc(true));
142 }
143 TEUCHOS_ASSERT(aztecooLOWSFPL != NULL);
144 aztecooLOWSFPL->sublist("Forward Solve").set("Tolerance",maxResid);
145 aztecooLOWSFPL->sublist("Adjoint Solve").set("Tolerance",maxResid);
146 if(showAllTests) {
147 aztecooLOWSFPL->set("Output Every RHS",bool(true));
148 }
149 if(out.get()) {
150 *out << "\naztecooLOWSFPL before setting parameters:\n";
151 aztecooLOWSFPL->print(OSTab(out).o(),0,true);
152 }
153 if(aztecooLOWSFPL) lowsFactory->setParameterList(Teuchos::rcp(aztecooLOWSFPL,false));
154
155 if(out.get()) *out << "\nC) Creating a AztecOOLinearOpWithSolve object nsA from A ...\n";
156
157 RCP<LinearOpWithSolveBase<double> >
158 nsA = lowsFactory->createOp();
159
160 Thyra::initializeOp<double>(*lowsFactory, A, nsA.ptr());
161
162 if(out.get()) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
163
164 LinearOpTester<double> linearOpTester;
165 linearOpTester.check_adjoint(testTranspose);
166 linearOpTester.num_random_vectors(numRandomVectors);
167 linearOpTester.set_all_error_tol(maxFwdError);
168 linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
169 linearOpTester.show_all_tests(showAllTests);
170 linearOpTester.dump_all(dumpAll);
171 Thyra::seed_randomize<double>(0);
172 result = linearOpTester.check(*nsA,out());
173 if(!result) success = false;
174
175 if(out.get()) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
176
177 LinearOpWithSolveTester<double> linearOpWithSolveTester;
178 linearOpWithSolveTester.turn_off_all_tests();
179 linearOpWithSolveTester.check_forward_default(true);
180 linearOpWithSolveTester.check_forward_residual(true);
181 if(testTranspose && useAztecPrec) {
182 linearOpWithSolveTester.check_adjoint_default(true);
183 linearOpWithSolveTester.check_adjoint_residual(true);
184 }
185 else {
186 linearOpWithSolveTester.check_adjoint_default(false);
187 linearOpWithSolveTester.check_adjoint_residual(false);
188 }
189 linearOpWithSolveTester.set_all_solve_tol(maxResid);
190 linearOpWithSolveTester.set_all_slack_error_tol(maxResid);
191 linearOpWithSolveTester.set_all_slack_warning_tol(10.0*maxResid);
192 linearOpWithSolveTester.forward_default_residual_error_tol(2.5*maxResid);
193 linearOpWithSolveTester.forward_default_solution_error_error_tol(maxSolutionError);
194 linearOpWithSolveTester.adjoint_default_residual_error_tol(2.5*maxResid);
195 linearOpWithSolveTester.adjoint_default_solution_error_error_tol(maxSolutionError);
196 linearOpWithSolveTester.show_all_tests(showAllTests);
197 linearOpWithSolveTester.dump_all(dumpAll);
198 Thyra::seed_randomize<double>(0);
199 result = linearOpWithSolveTester.check(*nsA,out.get());
200 if(!result) success = false;
201
202 if(out.get()) *out << "\nF) Uninitialize nsA, create preconditioner for diagonal scaled by 0.99 and then reinitialize nsA reusing the old preconditioner ...\n";
203
204 // Scale the diagonal of the matrix and then create the preconditioner for it
205 Thyra::uninitializeOp<double>(*lowsFactory, nsA.ptr());
206 // Above is not required but a good idea since we are changing the matrix
207 {
208 Epetra_Vector diag(epetra_A->RowMap());
209 epetra_A->ExtractDiagonalCopy(diag);
210 diag.Scale(0.5);
211 epetra_A->ReplaceDiagonalValues(diag);
212 }
213 Thyra::initializeOp<double>(*lowsFactory, A, nsA.ptr());
214
215 // Scale the matrix back again and then reuse the preconditioner
216 Thyra::uninitializeOp<double>(*lowsFactory, nsA.ptr());
217 // Above is not required but a good idea since we are changing the matrix
218 {
219 Epetra_Vector diag(epetra_A->RowMap());
220 epetra_A->ExtractDiagonalCopy(diag);
221 diag.Scale(1.0/0.5);
222 epetra_A->ReplaceDiagonalValues(diag);
223 }
224 initializeAndReuseOp<double>(*lowsFactory, A, nsA.ptr());
225
226 if(out.get()) *out << "\nG) Testing the LinearOpWithSolveBase interface of nsA ...\n";
227
228 Thyra::seed_randomize<double>(0);
229 result = linearOpWithSolveTester.check(*nsA,out.get());
230 if(!result) success = false;
231
232 if(useAztecPrec) {
233
234 if(out.get()) *out << "\nH) Reinitialize (A,A,PRECONDITIONER_INPUT_TYPE_AS_MATRIX) => nsA ...\n";
235
236 initializeApproxPreconditionedOp<double>(*lowsFactory, A, A, nsA.ptr());
237
238 if(out.get()) *out << "\nI) Testing the LinearOpWithSolveBase interface of nsA ...\n";
239
240 Thyra::seed_randomize<double>(0);
241 result = linearOpWithSolveTester.check(*nsA,out.get());
242 if(!result) success = false;
243
244 if(testTranspose && useAztecPrec) {
245 linearOpWithSolveTester.check_adjoint_default(true);
246 linearOpWithSolveTester.check_adjoint_residual(true);
247 }
248 else {
249 linearOpWithSolveTester.check_adjoint_default(false);
250 linearOpWithSolveTester.check_adjoint_residual(false);
251 }
252
253 }
254 else {
255
256 if(out.get()) *out << "\nSkipping testing steps H and I since we are not using aztec preconditioning and therefore will not test with an external preconditioner matrix!\n";
257
258 }
259
260
261 RCP<PreconditionerFactoryBase<double> >
262 precFactory;
263
264#ifdef HAVE_AZTECOO_IFPACK
265
266 if(useAztecPrec) {
267
268 if(testTranspose) {
269 linearOpWithSolveTester.check_adjoint_default(true);
270 linearOpWithSolveTester.check_adjoint_residual(true);
271 }
272
273 if(out.get()) *out << "\nJ) Create an ifpack preconditioner precA for A ...\n";
274
275 precFactory = Teuchos::rcp(new IfpackPreconditionerFactory());
276
277 if(out.get()) {
278 *out << "\nprecFactory.description() = " << precFactory->description() << std::endl;
279 *out << "\nprecFactory.getValidParameters() =\n";
280 precFactory->getValidParameters()->print(OSTab(out).o(),0,true,false);
281 }
282
283 RCP<Teuchos::ParameterList>
284 ifpackPFPL = Teuchos::rcp(new Teuchos::ParameterList("IfpackPreconditionerFactory"));
285 ifpackPFPL->set("Prec Type","ILUT");
286 ifpackPFPL->set("Overlap",int(1));
287 if(out.get()) {
288 *out << "\nifpackPFPL before setting parameters =\n";
289 ifpackPFPL->print(OSTab(out).o(),0,true);
290 }
291
292 precFactory->setParameterList(ifpackPFPL);
293
294 RCP<PreconditionerBase<double> >
295 precA = precFactory->createPrec();
296 Thyra::initializePrec<double>(*precFactory,A,&*precA);
297
298 if(out.get()) {
299 *out << "\nifpackPFPL after setting parameters =\n";
300 ifpackPFPL->print(OSTab(out).o(),0,true);
301 *out << "\nprecFactory.description() = " << precFactory->description() << std::endl;
302 }
303
304 if(out.get()) *out << "\nprecA.description() = " << precA->description() << std::endl;
305 if(out.get() && dumpAll) *out << "\ndescribe(precA) =\n" << describe(*precA,Teuchos::VERB_EXTREME);
306
307 if(out.get()) *out << "\nK) Reinitialize (A,precA->getUnspecifiedPrecOp(),PRECONDITIONER_INPUT_TYPE_AS_OPERATOR) => nsA ...\n";
308
309 Thyra::initializePreconditionedOp<double>(*lowsFactory,A,precA,&*nsA);
310
311 if(out.get()) *out << "\nL) Testing the LinearOpWithSolveBase interface of nsA ...\n";
312
313 Thyra::seed_randomize<double>(0);
314 result = linearOpWithSolveTester.check(*nsA,out.get());
315 if(!result) success = false;
316
317 if(testTranspose && useAztecPrec) {
318 linearOpWithSolveTester.check_adjoint_default(true);
319 linearOpWithSolveTester.check_adjoint_residual(true);
320 }
321 else {
322 linearOpWithSolveTester.check_adjoint_default(false);
323 linearOpWithSolveTester.check_adjoint_residual(false);
324 }
325
326 }
327 else {
328
329 if(out.get()) *out << "\nSkipping testing steps J, K, and L since we are not using aztec preconditioning and therefore will not test with an ifpack preconditioner!\n";
330
331 }
332
333#else // HAVE_AZTECOO_IFPACK
334
335 if(out.get()) *out << "\nSkipping testing steps J, K, and L since they require ifpack and ifpack has not been enabled!\n";
336
337#endif // HAVE_AZTECOO_IFPACK
338
339
340 if(out.get()) *out << "\nM) Scale the epetra_A object by 2.5, and then reinitialize nsA with epetra_A ...\n";
341
342 Thyra::uninitializeOp<double>(*lowsFactory, nsA.ptr());
343 // Not required but a good idea since we are changing the matrix
344 epetra_A->Scale(2.5);
345 initializeOp<double>(*lowsFactory, A, nsA.ptr());
346
347 if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA ...\n";
348
349 Thyra::seed_randomize<double>(0);
350 result = linearOpWithSolveTester.check(*nsA,out.get());
351 if(!result) success = false;
352
353 if(out.get()) *out << "\nO) Create a scaled (by 2.5) copy epetra_A2 of epetra_A, and then reinitialize nsA with epetra_A2 ...\n";
354
355 RCP<Epetra_CrsMatrix>
356 epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
357 epetra_A2->Scale(2.5);
358 RCP<const LinearOpBase<double> >
359 A2 = Thyra::epetraLinearOp(epetra_A2);
360 initializeOp<double>(*lowsFactory, A2, nsA.ptr());
361 // Note that it was okay not to uninitialize nsA first here since A, which
362 // was used to initialize nsA last, was not changed and therefore the
363 // state of nsA was fine throughout
364
365 if(out.get()) *out << "\nP) Testing the LinearOpWithSolveBase interface of nsA ...\n";
366
367 Thyra::seed_randomize<double>(0);
368 result = linearOpWithSolveTester.check(*nsA,out.get());
369 if(!result) success = false;
370
371 if(!useAztecPrec) {
372
373 if(out.get()) *out << "\nQ) Create an implicitly scaled (by 2.5) and transposed matrix A3 = scale(2.5,transpose(A)) and initialize nsA2 ...\n";
374
375 RCP<const LinearOpBase<double> >
376 A3 = scale<double>(2.5,transpose<double>(A));
377 RCP<LinearOpWithSolveBase<double> >
378 nsA2 = linearOpWithSolve(*lowsFactory,A3);
379
380 if(out.get()) *out << "\nR) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
381
382 Thyra::seed_randomize<double>(0);
383 result = linearOpWithSolveTester.check(*nsA2,out.get());
384 if(!result) success = false;
385
386 if(out.get()) *out << "\nS) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
387
388 result = linearOpTester.compare(
389 *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
390 ,out()
391 );
392 if(!result) success = false;
393
394 }
395 else {
396
397 if(out.get()) *out << "\nSkipping testing steps Q, R, and S because we are using internal AztecOO preconditioners!\n";
398
399 }
400
401 if(out.get()) *out << "\nT) Running example use cases ...\n";
402
403 nonExternallyPreconditionedLinearSolveUseCases(
404 *A,*lowsFactory,testTranspose,*out
405 );
406
407 if(precFactory.get()) {
408 externallyPreconditionedLinearSolveUseCases(
409 *A,*lowsFactory,*precFactory,false,true,*out
410 );
411 }
412
413#else // SUN_CXX
414
415 if(out.get()) *out << "\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
416 success = false;
417
418#endif // SUN_CXX
419
420 }
421 catch( const std::exception &excpt ) {
422 std::cerr << "\n*** Caught standard exception : " << excpt.what() << std::endl;
423 success = false;
424 }
425
426 return success;
427
428}
LinearOpWithSolveFactoryBase subclass implemented in terms of AztecOO.
Concrete preconditioner factory subclass based on Ifpack.
bool test_single_aztecoo_thyra_solver(const std::string matrixFile, const bool testTranspose, const int numRandomVectors, const double maxFwdError, const double maxResid, const double maxSolutionError, const bool showAllTests, const bool dumpAll, Teuchos::ParameterList *paramList, Teuchos::FancyOStream *out)
Testing function for a single aztecoo solver with a single matrix.