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