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 //
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 "test_single_amesos_thyra_solver.hpp"
43 
44 #ifndef SUN_CXX
45 
46 #include "Thyra_AmesosLinearOpWithSolveFactory.hpp"
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 
63 bool Thyra::test_single_amesos_thyra_solver(
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 }
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 by   doxygen 1.8.5