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" 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
69 ,Teuchos::FancyOStream *out_arg
76 bool result, success =
true;
78 RCP<Teuchos::FancyOStream>
79 out = Teuchos::rcp(out_arg,
false);
86 <<
"\n*** Testing Thyra::AmesosLinearOpWithSolveFactory (and Thyra::AmesosLinearOpWithSolve)" 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
101 <<
"amesosLOWSFPL:\n";
102 amesosLOWSFPL->print(OSTab(out).o(),0,
true);
106 if(out.get()) *out <<
"\nA) Reading in an epetra matrix A from the file \'"<<matrixFile<<
"\' ...\n";
108 Epetra_SerialComm comm;
109 RCP<Epetra_CrsMatrix> epetra_A;
110 EpetraExt::readEpetraLinearSystem( matrixFile, comm, &epetra_A );
112 RCP<const LinearOpBase<double> >
113 A = Thyra::epetraLinearOp(epetra_A);
115 if(out.get() && dumpAll) *out <<
"\ndescribe(A) =\n" << describe(*A,Teuchos::VERB_EXTREME);
117 if(out.get()) *out <<
"\nB) Creating a AmesosLinearOpWithSolveFactory object lowsFactory ...\n";
119 RCP<LinearOpWithSolveFactoryBase<double> >
122 lowsFactory->setParameterList(Teuchos::rcp(amesosLOWSFPL,
false));
125 *out <<
"\nlowsFactory.description() = " << lowsFactory->description() << std::endl;
126 *out <<
"\nlowsFactory.getValidParameters() =\n";
127 lowsFactory->getValidParameters()->print(OSTab(out).o(),0,
true,
false);
130 if(out.get()) *out <<
"\nC) Creating a AmesosLinearOpWithSolve object nsA ...\n";
132 RCP<LinearOpWithSolveBase<double> >
133 nsA = lowsFactory->createOp();
135 initializeOp<double>(*lowsFactory, A, nsA.ptr());
137 bool testTranspose = testTranspose_in;
138 if (testTranspose_in && !Thyra::solveSupports(*nsA, Thyra::CONJTRANS)) {
140 <<
"\nChanging testTranspose=false since " 141 << nsA->description() <<
" does not support adjoint solves!\n";
142 testTranspose =
false;
145 if(out.get()) *out <<
"\nD) Testing the LinearOpBase interface of nsA ...\n";
147 Thyra::seed_randomize<double>(0);
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;
159 if(out.get()) *out <<
"\nE) Testing the LinearOpWithSolveBase interface of nsA ...\n";
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);
181 LinearOpWithSolveTester<double> adjLinearOpWithSolveTester(linearOpWithSolveTester);
182 adjLinearOpWithSolveTester.check_forward_default(testTranspose);
183 adjLinearOpWithSolveTester.check_forward_residual(testTranspose);
185 result = linearOpWithSolveTester.check(*nsA,out.get());
186 if(!result) success =
false;
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";
190 uninitializeOp<double>(*lowsFactory, nsA.ptr());
191 epetra_A->Scale(2.5);
192 initializeOp<double>(*lowsFactory, A, nsA.ptr());
194 if(out.get()) *out <<
"\nG) Testing the LinearOpBase interface of nsA ...\n";
196 Thyra::seed_randomize<double>(0);
198 result = linearOpTester.check(*nsA, out());
199 if(!result) success =
false;
201 if(out.get()) *out <<
"\nH) Testing the LinearOpWithSolveBase interface of nsA ...\n";
203 result = linearOpWithSolveTester.check(*nsA,out.get());
204 if(!result) success =
false;
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";
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());
215 if(out.get()) *out <<
"\nJ) Testing the LinearOpBase interface of nsA ...\n";
217 Thyra::seed_randomize<double>(0);
219 result = linearOpTester.check(*nsA, out());
220 if(!result) success =
false;
222 if(out.get()) *out <<
"\nK) Testing the LinearOpWithSolveBase interface of nsA ...\n";
224 result = linearOpWithSolveTester.check(*nsA, out.get());
225 if(!result) success =
false;
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";
229 RCP<const LinearOpBase<double> >
230 A3 = scale<double>(2.5,Thyra::transpose<double>(A));
231 RCP<LinearOpWithSolveBase<double> >
232 nsA2 = linearOpWithSolve(*lowsFactory,A3);
234 if(out.get()) *out <<
"\nM) Testing the LinearOpBase interface of nsA2 ...\n";
236 Thyra::seed_randomize<double>(0);
238 result = linearOpTester.check(*nsA2,out());
239 if(!result) success =
false;
241 if(out.get()) *out <<
"\nN) Testing the LinearOpWithSolveBase interface of nsA2 ...\n";
243 result = adjLinearOpWithSolveTester.check(*nsA2,out.get());
244 if(!result) success =
false;
246 if(out.get()) *out <<
"\nO) Testing that LinearOpBase interfaces of transpose(nsA) == nsA2 ...\n";
248 result = linearOpTester.compare(
249 *transpose(Teuchos::rcp_implicit_cast<
const LinearOpBase<double> >(nsA)),*nsA2
252 if(!result) success =
false;
254 if(out.get()) *out <<
"\nP) Running example use cases ...\n";
256 nonExternallyPreconditionedLinearSolveUseCases(
257 *A,*lowsFactory,
true,*out
260 if(out.get()) *out <<
"\nQ) Creating a DefaultInverseLinearOp object from nsA and testing the LinearOpBase interface ...\n";
262 RCP<const LinearOpBase<double> >
263 invA = inverse<double>(nsA.getConst());
265 result = linearOpTester.check(*invA,out());
266 if(!result) success =
false;
270 if(out.get()) *out <<
"\nTest failed since is was not even compiled since SUN_CXX was defined!\n";
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.