Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test_single_amesos_thyra_solver.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stratimikos: Thyra-based strategies for linear solvers
4 //
5 // Copyright 2006 NTESS and the Stratimikos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "test_single_amesos_thyra_solver.hpp"
11 
12 #ifndef SUN_CXX
13 
14 #include "Thyra_AmesosLinearOpWithSolveFactory.hpp"
15 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
16 #include "Thyra_LinearOpWithSolveFactoryExamples.hpp"
17 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
18 #include "Thyra_DefaultInverseLinearOp.hpp"
19 #include "Thyra_EpetraLinearOp.hpp"
20 #include "Thyra_LinearOpTester.hpp"
21 #include "Thyra_LinearOpWithSolveTester.hpp"
22 #include "EpetraExt_readEpetraLinearSystem.h"
23 #ifdef EPETRA_MPI
24 #include "Epetra_MpiComm.h"
25 #else
26 #include "Epetra_SerialComm.h"
27 #endif
28 
29 #endif // SUN_CXX
30 
31 bool Thyra::test_single_amesos_thyra_solver(
32  const std::string matrixFile
33  ,Teuchos::ParameterList *amesosLOWSFPL
34  ,const bool testTranspose_in
35  ,const int numRandomVectors
36  ,const double maxFwdError
37  ,const double maxError
38  ,const double maxResid
39  ,const bool showAllTests
40  ,const bool dumpAll
41  ,Teuchos::FancyOStream *out_arg
42  )
43 {
44 
45  using Teuchos::RCP;
46  using Teuchos::OSTab;
47 
48  bool result, success = true;
49 
50  RCP<Teuchos::FancyOStream>
51  out = Teuchos::rcp(out_arg,false);
52 
53 #ifndef SUN_CXX
54 
55  if(out.get()) {
56  *out
57  << "\n***"
58  << "\n*** Testing Thyra::AmesosLinearOpWithSolveFactory (and Thyra::AmesosLinearOpWithSolve)"
59  << "\n***\n"
60  << "\nEchoing input options:"
61  << "\n matrixFile = " << matrixFile
62  << "\n testTranspose = " << testTranspose_in
63  << "\n numRandomVectors = " << numRandomVectors
64  << "\n maxFwdError = " << maxFwdError
65  << "\n maxError = " << maxError
66  << "\n maxResid = " << maxResid
67  << "\n showAllTests = " << showAllTests
68  << "\n dumpAll = " << dumpAll
69  << std::endl;
70  if(amesosLOWSFPL) {
71  OSTab tab(out);
72  *out
73  << "amesosLOWSFPL:\n";
74  amesosLOWSFPL->print(OSTab(out).o(),0,true);
75  }
76  }
77 
78  if(out.get()) *out << "\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<"\' ...\n";
79 
80 #ifdef EPETRA_MPI
81  Epetra_MpiComm comm(MPI_COMM_WORLD);
82 #else
83  Epetra_SerialComm comm;
84 #endif
85  RCP<Epetra_CrsMatrix> epetra_A;
86  EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
87 
88  RCP<const LinearOpBase<double> >
89  A = Thyra::epetraLinearOp(epetra_A);
90 
91  if(out.get() && dumpAll) *out << "\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
92 
93  if(out.get()) *out << "\nB) Creating a AmesosLinearOpWithSolveFactory object lowsFactory ...\n";
94 
95  RCP<LinearOpWithSolveFactoryBase<double> >
96  lowsFactory = Teuchos::rcp(new AmesosLinearOpWithSolveFactory());
97 
98  lowsFactory->setParameterList(Teuchos::rcp(amesosLOWSFPL,false));
99 
100  if(out.get()) {
101  *out << "\nlowsFactory.description() = " << lowsFactory->description() << std::endl;
102  *out << "\nlowsFactory.getValidParameters() =\n";
103  lowsFactory->getValidParameters()->print(OSTab(out).o(),0,true,false);
104  }
105 
106  if(out.get()) *out << "\nC) Creating a AmesosLinearOpWithSolve object nsA ...\n";
107 
108  RCP<LinearOpWithSolveBase<double> >
109  nsA = lowsFactory->createOp();
110 
111  initializeOp<double>(*lowsFactory, A, nsA.ptr());
112 
113  bool testTranspose = testTranspose_in;
114  if (testTranspose_in && !Thyra::solveSupports(*nsA, Thyra::CONJTRANS)) {
115  if(out.get()) *out
116  << "\nChanging testTranspose=false since "
117  << nsA->description() << " does not support adjoint solves!\n";
118  testTranspose = false;
119  }
120 
121  if(out.get()) *out << "\nD) Testing the LinearOpBase interface of nsA ...\n";
122 
123  Thyra::seed_randomize<double>(0);
124 
125  LinearOpTester<double> linearOpTester;
126  linearOpTester.check_adjoint(testTranspose);
127  linearOpTester.num_random_vectors(numRandomVectors);
128  linearOpTester.set_all_error_tol(maxFwdError);
129  linearOpTester.set_all_warning_tol(1e-2*maxFwdError);
130  linearOpTester.show_all_tests(showAllTests);
131  linearOpTester.dump_all(dumpAll);
132  result = linearOpTester.check(*nsA,out());
133  if(!result) success = false;
134 
135  if(out.get()) *out << "\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
136 
137  LinearOpWithSolveTester<double> linearOpWithSolveTester;
138  linearOpWithSolveTester.turn_off_all_tests();
139  linearOpWithSolveTester.check_forward_default(true);
140  linearOpWithSolveTester.forward_default_residual_error_tol(1.1*maxResid);
141  linearOpWithSolveTester.forward_default_residual_warning_tol(2.0*maxResid);
142  linearOpWithSolveTester.check_forward_residual(true);
143  linearOpWithSolveTester.forward_residual_solve_tol(maxResid);
144  linearOpWithSolveTester.forward_residual_slack_error_tol(1e-1*maxResid);
145  linearOpWithSolveTester.forward_residual_slack_warning_tol(maxResid);
146  linearOpWithSolveTester.check_adjoint_default(testTranspose);
147  linearOpWithSolveTester.adjoint_default_residual_error_tol(1.1*maxResid);
148  linearOpWithSolveTester.adjoint_default_residual_warning_tol(2.0*maxResid);
149  linearOpWithSolveTester.check_adjoint_residual(testTranspose);
150  linearOpWithSolveTester.adjoint_residual_solve_tol(maxResid);
151  linearOpWithSolveTester.adjoint_residual_slack_error_tol(1e-1*maxResid);
152  linearOpWithSolveTester.adjoint_residual_slack_warning_tol(maxResid);
153  linearOpWithSolveTester.num_random_vectors(numRandomVectors);
154  linearOpWithSolveTester.show_all_tests(showAllTests);
155  linearOpWithSolveTester.dump_all(dumpAll);
156 
157  LinearOpWithSolveTester<double> adjLinearOpWithSolveTester(linearOpWithSolveTester);
158  adjLinearOpWithSolveTester.check_forward_default(testTranspose);
159  adjLinearOpWithSolveTester.check_forward_residual(testTranspose);
160 
161  result = linearOpWithSolveTester.check(*nsA,out.get());
162  if(!result) success = false;
163 
164  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";
165 
166  uninitializeOp<double>(*lowsFactory, nsA.ptr()); // Optional call but a good idea if changing the operator
167  epetra_A->Scale(2.5);
168  initializeOp<double>(*lowsFactory, A, nsA.ptr());
169 
170  if(out.get()) *out << "\nG) Testing the LinearOpBase interface of nsA ...\n";
171 
172  Thyra::seed_randomize<double>(0);
173 
174  result = linearOpTester.check(*nsA, out());
175  if(!result) success = false;
176 
177  if(out.get()) *out << "\nH) Testing the LinearOpWithSolveBase interface of nsA ...\n";
178 
179  result = linearOpWithSolveTester.check(*nsA,out.get());
180  if(!result) success = false;
181 
182  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";
183 
184  RCP<Epetra_CrsMatrix>
185  epetra_A2 = Teuchos::rcp(new Epetra_CrsMatrix(*epetra_A));
186  epetra_A2->Scale(2.5);
187  RCP<const LinearOpBase<double> >
188  A2 = Thyra::epetraLinearOp(epetra_A2);
189  initializeOp<double>(*lowsFactory, A2, nsA.ptr());
190 
191  if(out.get()) *out << "\nJ) Testing the LinearOpBase interface of nsA ...\n";
192 
193  Thyra::seed_randomize<double>(0);
194 
195  result = linearOpTester.check(*nsA, out());
196  if(!result) success = false;
197 
198  if(out.get()) *out << "\nK) Testing the LinearOpWithSolveBase interface of nsA ...\n";
199 
200  result = linearOpWithSolveTester.check(*nsA, out.get());
201  if(!result) success = false;
202 
203  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";
204 
205  RCP<const LinearOpBase<double> >
206  A3 = scale<double>(2.5,Thyra::transpose<double>(A));
207  RCP<LinearOpWithSolveBase<double> >
208  nsA2 = linearOpWithSolve(*lowsFactory,A3);
209 
210  if(out.get()) *out << "\nM) Testing the LinearOpBase interface of nsA2 ...\n";
211 
212  Thyra::seed_randomize<double>(0);
213 
214  result = linearOpTester.check(*nsA2,out());
215  if(!result) success = false;
216 
217  if(out.get()) *out << "\nN) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
218 
219  result = adjLinearOpWithSolveTester.check(*nsA2,out.get());
220  if(!result) success = false;
221 
222  if(out.get()) *out << "\nO) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
223 
224  result = linearOpTester.compare(
225  *transpose(Teuchos::rcp_implicit_cast<const LinearOpBase<double> >(nsA)),*nsA2
226  ,out()
227  );
228  if(!result) success = false;
229 
230  if(out.get()) *out << "\nP) Running example use cases ...\n";
231 
232  nonExternallyPreconditionedLinearSolveUseCases(
233  *A,*lowsFactory,true,*out
234  );
235 
236  if(out.get()) *out << "\nQ) Creating a DefaultInverseLinearOp object from nsA and testing the LinearOpBase interface ...\n";
237 
238  RCP<const LinearOpBase<double> >
239  invA = inverse<double>(nsA.getConst());
240 
241  result = linearOpTester.check(*invA,out());
242  if(!result) success = false;
243 
244 #else // SUN_CXX
245 
246  if(out.get()) *out << "\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
247  success = false;
248 
249 #endif // SUN_CXX
250 
251  return success;
252 
253 }
Thyra_EpetraLinearOp.hpp
Teuchos::OSTab
basic_OSTab< char > OSTab
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::basic_FancyOStream
test_single_amesos_thyra_solver.hpp
Thyra_AmesosLinearOpWithSolveFactory.hpp
Generated on Fri Dec 20 2024 09:21:48 for Stratimikos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5