Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_LinearOpWithSolveTester_def.hpp
1// @HEADER
2// *****************************************************************************
3// Thyra: Interfaces and Support for Abstract Numerical Algorithms
4//
5// Copyright 2004 NTESS and the Thyra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
11#define THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
12
13
14#include "Thyra_LinearOpWithSolveTester_decl.hpp"
15#include "Thyra_LinearOpWithSolveBase.hpp"
16#include "Thyra_DefaultScaledAdjointLinearOp.hpp"
17#include "Thyra_describeLinearOp.hpp"
18#include "Thyra_VectorStdOps.hpp"
19#include "Thyra_MultiVectorStdOps.hpp"
20#include "Thyra_TestingTools.hpp"
21#include "Teuchos_Time.hpp"
22#include "Teuchos_TestingHelpers.hpp"
23
24#include <functional>
25
26namespace Thyra {
27
28
29// Constructors/initializers
30
31
32template<class Scalar>
33LinearOpWithSolveTester<Scalar>::LinearOpWithSolveTester()
34 :check_forward_default_(check_forward_default_default_),
35 forward_default_residual_warning_tol_(warning_tol_default_),
36 forward_default_residual_error_tol_(error_tol_default_),
37 forward_default_solution_error_warning_tol_(warning_tol_default_),
38 forward_default_solution_error_error_tol_(error_tol_default_),
39 check_forward_residual_(check_forward_residual_default_),
40 forward_residual_solve_tol_(solve_tol_default_),
41 forward_residual_slack_warning_tol_(slack_warning_tol_default_),
42 forward_residual_slack_error_tol_(slack_error_tol_default_),
43 check_adjoint_default_(check_adjoint_default_default_),
44 adjoint_default_residual_warning_tol_(warning_tol_default_),
45 adjoint_default_residual_error_tol_(error_tol_default_),
46 adjoint_default_solution_error_warning_tol_(warning_tol_default_),
47 adjoint_default_solution_error_error_tol_(error_tol_default_),
48 check_adjoint_residual_(check_adjoint_residual_default_),
49 adjoint_residual_solve_tol_(solve_tol_default_),
50 adjoint_residual_slack_warning_tol_(slack_warning_tol_default_),
51 adjoint_residual_slack_error_tol_(slack_error_tol_default_),
52 num_random_vectors_(num_random_vectors_default_),
53 show_all_tests_(show_all_tests_default_),
54 dump_all_(dump_all_default_),
55 num_rhs_(num_rhs_default_)
56{}
57
58
59template<class Scalar>
60void LinearOpWithSolveTester<Scalar>::turn_off_all_tests()
61{
62 check_forward_default_ = false;
63 check_forward_residual_ = false;
64 check_adjoint_default_ = false;
65 check_adjoint_residual_ = false;
66}
67
68
69template<class Scalar>
70void
71LinearOpWithSolveTester<Scalar>::set_all_solve_tol(
72 const ScalarMag solve_tol )
73{
74 forward_residual_solve_tol_ = solve_tol;
75 forward_residual_solve_tol_ = solve_tol;
76 adjoint_residual_solve_tol_ = solve_tol;
77}
78
79
80template<class Scalar>
81void
82LinearOpWithSolveTester<Scalar>::set_all_slack_warning_tol(
83 const ScalarMag slack_warning_tol )
84{
85 forward_default_residual_warning_tol_ = slack_warning_tol;
86 forward_default_solution_error_warning_tol_ = slack_warning_tol;
87 forward_residual_slack_warning_tol_ = slack_warning_tol;
88 adjoint_default_residual_warning_tol_ = slack_warning_tol;
89 adjoint_default_solution_error_warning_tol_ = slack_warning_tol;
90 adjoint_residual_slack_warning_tol_ = slack_warning_tol;
91}
92
93
94template<class Scalar>
95void
96LinearOpWithSolveTester<Scalar>::set_all_slack_error_tol(
97 const ScalarMag slack_error_tol )
98{
99 forward_default_residual_error_tol_ = slack_error_tol;
100 forward_default_solution_error_error_tol_ = slack_error_tol;
101 forward_residual_slack_error_tol_ = slack_error_tol;
102 adjoint_default_residual_error_tol_ = slack_error_tol;
103 adjoint_default_solution_error_error_tol_ = slack_error_tol;
104 adjoint_residual_slack_error_tol_ = slack_error_tol;
105}
106
107
108// Overridden from ParameterListAcceptor
109
110
111template<class Scalar>
112void LinearOpWithSolveTester<Scalar>::setParameterList(
113 const RCP<ParameterList>& paramList )
114{
115 using Teuchos::getParameter;
116 ParameterList &pl = *paramList;
117 this->setMyParamList(paramList);
118 paramList->validateParametersAndSetDefaults(*getValidParameters());
119 set_all_solve_tol(getParameter<ScalarMag>(pl, AllSolveTol_name_));
120 set_all_slack_warning_tol(getParameter<ScalarMag>(pl, AllSlackWarningTol_name_));
121 set_all_slack_error_tol(getParameter<ScalarMag>(pl, AllSlackErrorTol_name_));
122 show_all_tests(getParameter<bool>(pl, ShowAllTests_name_));
123 dump_all(getParameter<bool>(pl, DumpAll_name_));
124 // ToDo: Add more parameters as you need them
125}
126
127
128template<class Scalar>
129RCP<const ParameterList>
130LinearOpWithSolveTester<Scalar>::getValidParameters() const
131{
132 static RCP<const ParameterList> validPL;
133 if (is_null(validPL) ) {
134 RCP<ParameterList> pl = Teuchos::parameterList();
135 pl->set(AllSolveTol_name_, solve_tol_default_,
136 "Sets all of the solve tolerances to the same value. Note: This option\n"
137 "is applied before any specific tolerance which may override it.");
138 pl->set(AllSlackWarningTol_name_, slack_warning_tol_default_,
139 "Sets all of the slack warning values to the same value. Note: This option\n"
140 "is applied before any specific tolerance which may override it.");
141 pl->set(AllSlackErrorTol_name_, slack_error_tol_default_,
142 "Sets all of the slack error values to the same value. Note: This option\n"
143 "is applied before any specific tolerance which may override it.");
144 pl->set(ShowAllTests_name_, false,
145 "If true, then all tests be traced to the output stream.");
146 pl->set(DumpAll_name_, false,
147 "If true, then all qualtities will be dumped to the output stream. Warning!\n"
148 "only do this to debug smaller problems as this can create a lot of output");
149 // ToDo: Add more parameters as you need them (don't forget to test them)
150 validPL = pl;
151 }
152 return validPL;
153}
154
155
156// LOWS testing
157
158
159template<class Scalar>
160bool LinearOpWithSolveTester<Scalar>::check(
161 const LinearOpWithSolveBase<Scalar> &op,
162 Teuchos::FancyOStream *out_arg ) const
163{
164
165 using std::endl;
166 using Teuchos::as;
167 using Teuchos::optInArg;
168 using Teuchos::inoutArg;
169 using Teuchos::FancyOStream;
170 using Teuchos::OSTab;
171 typedef Teuchos::VerboseObjectTempState<LinearOpWithSolveBase<Scalar> > VOTS;
172 typedef Teuchos::ScalarTraits<Scalar> ST;
173
174 bool success = true, result;
175 const int l_num_rhs = this->num_rhs();
176 Teuchos::RCP<FancyOStream> out = Teuchos::rcp(out_arg,false);
177 const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
178 Teuchos::Time timer("");
179
180 OSTab tab(out,1,"THYRA");
181
182 if(out.get()) {
183 *out <<endl<< "*** Entering LinearOpWithSolveTester<"<<ST::name()<<">::check(op,...) ...\n";
184 if(show_all_tests()) {
185 *out <<endl<< "describe forward op:\n" << Teuchos::describe(op,verbLevel);
186 if(opSupported(op, CONJTRANS) && verbLevel==Teuchos::VERB_EXTREME) {
187 *out <<endl<< "describe adjoint op:\n";
188 describeLinearOp<Scalar>(
189 *adjoint(RCP<const LinearOpBase<Scalar> >(Teuchos::rcp(&op,false))),
190 *out, verbLevel
191 );
192 }
193 }
194 else {
195 *out <<endl<< "describe op: " << op.description() << endl;
196 }
197 }
198
199 Teuchos::RCP<const VectorSpaceBase<Scalar> > range = op.range();
200 Teuchos::RCP<const VectorSpaceBase<Scalar> > domain = op.domain();
201
202 if( check_forward_default() ) {
203
204 if(out.get()) *out <<endl<< "this->check_forward_default()==true: Checking the default forward solve ... ";
205
206 TestResultsPrinter testResultsPrinter(out, show_all_tests());
207 const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
208
209 bool these_results = true;
210
211 result = true;
212 TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(NOTRANS), true, *testOut, result);
213 if (!result) these_results = false;
214
215 if(result) {
216
217 *testOut
218 <<endl<< "Checking that the forward default solve matches the forward operator:\n"
219 <<endl<< " inv(Op)*Op*v1 == v1"
220 <<endl<< " \\___/"
221 <<endl<< " v2"
222 <<endl<< " \\___________/"
223 <<endl<< " v3"
224 <<endl<< ""
225 <<endl<< " v4 = v3-v1"
226 <<endl<< " v5 = Op*v3-v2"
227 <<endl<< ""
228 <<endl<< " norm(v4)/norm(v1) <= forward_default_solution_error_error_tol()"
229 <<endl<< " norm(v5)/norm(v2) <= forward_default_residual_error_tol()"
230 <<endl;
231
232 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
233
234 OSTab tab2(testOut);
235
236 *testOut <<endl<< "Random vector tests = " << rand_vec_i << endl;
237
238 tab.incrTab();
239
240 *testOut <<endl<< "v1 = randomize(-1,+1); ...\n" ;
241 Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,l_num_rhs);
242 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.ptr() );
243 if(dump_all()) *testOut <<endl<< "v1 =\n" << describe(*v1,verbLevel);
244
245 *testOut <<endl<< "v2 = Op*v1 ...\n" ;
246 Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,l_num_rhs);
247 timer.start(true);
248 Thyra::apply(op, NOTRANS, *v1, v2.ptr());
249 timer.stop();
250 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
251 if(dump_all()) *testOut <<endl<< "v2 =\n" << describe(*v2,verbLevel);
252
253 *testOut <<endl<< "v3 = inv(Op)*v2 ...\n" ;
254 Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,l_num_rhs);
255 assign(v3.ptr(), ST::zero());
256 SolveStatus<Scalar> solveStatus;
257 {
258 VOTS lowsTempState(Teuchos::rcp(&op,false),testOut,verbLevel);
259 timer.start(true);
260 solveStatus = solve<Scalar>(op, NOTRANS, *v2, v3.ptr());
261 timer.stop();
262 OSTab(testOut).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
263 }
264 if(dump_all()) *testOut <<endl<< "v3 =\n" << describe(*v3,verbLevel);
265 *testOut
266 <<endl<< "solve status:\n";
267 OSTab(testOut).o() << solveStatus;
268
269 *testOut <<endl<< "v4 = v3 - v1 ...\n" ;
270 Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,l_num_rhs);
271 V_VmV( v4.ptr(), *v3, *v1 );
272 if(dump_all()) *testOut <<endl<< "v4 =\n" << describe(*v4,verbLevel);
273
274 *testOut <<endl<< "v5 = Op*v3 - v2 ...\n" ;
275 Teuchos::RCP<MultiVectorBase<Scalar> > v5 = createMembers(range, l_num_rhs);
276 assign( v5.ptr(), *v2 );
277 timer.start(true);
278 Thyra::apply(op, NOTRANS, *v3, v5.ptr(), as<Scalar>(1.0) ,as<Scalar>(-1.0));
279 timer.stop();
280 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
281 if(dump_all()) *testOut <<endl<< "v5 =\n" << describe(*v5,verbLevel);
282
283 Array<ScalarMag> norms_v1(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v1(l_num_rhs);
284 norms(*v1, norms_v1());
285 norms(*v4, norms_v4());
286 std::transform(
287 norms_v4.begin(), norms_v4.end(), norms_v1.begin(), norms_v4_norms_v1.begin()
288 ,std::divides<ScalarMag>()
289 );
290
291 result = testMaxErrors<Scalar>(
292 "norm(v4)/norm(v1)", norms_v4_norms_v1(),
293 "forward_default_solution_error_error_tol()",
294 forward_default_solution_error_error_tol(),
295 "forward_default_solution_error_warning_tol()",
296 forward_default_solution_error_warning_tol(),
297 testOut.ptr()
298 );
299 if(!result) these_results = false;
300
301 Array<ScalarMag> norms_v2(l_num_rhs), norms_v5(l_num_rhs), norms_v5_norms_v2(l_num_rhs);
302 norms(*v2, norms_v2());
303 norms(*v5, norms_v5());
304 std::transform(
305 norms_v5.begin(), norms_v5.end(), norms_v2.begin(), norms_v5_norms_v2.begin()
306 ,std::divides<ScalarMag>()
307 );
308
309 result = testMaxErrors<Scalar>(
310 "norm(v5)/norm(v2)", norms_v5_norms_v2(),
311 "forward_default_residual_error_tol()",
312 forward_default_residual_error_tol(),
313 "forward_default_residual_warning_tol()",
314 forward_default_residual_warning_tol(),
315 testOut.ptr()
316 );
317 if(!result) these_results = false;
318
319 }
320 }
321 else {
322 *testOut <<endl<< "Forward operator not supported, skipping check!\n";
323 }
324
325 testResultsPrinter.printTestResults(these_results, inoutArg(success));
326
327 }
328 else {
329 if(out.get()) *out <<endl<< "this->check_forward_default()==false: Skipping the check of the default forward solve!\n";
330 }
331
332 if( check_forward_residual() ) {
333
334 if(out.get()) *out <<endl<< "this->check_forward_residual()==true: Checking the forward solve with a tolerance on the residual ... ";
335
336 TestResultsPrinter testResultsPrinter(out, show_all_tests());
337 const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
338
339 bool these_results = true;
340
341 result = true;
342 TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(NOTRANS), true, *testOut, result);
343 if (!result) these_results = false;
344
345 if(result) {
346
347 *testOut
348 <<endl<< "Checking that the forward solve matches the forward operator to a residual tolerance:\n"
349 <<endl<< " v3 = inv(Op)*Op*v1"
350 <<endl<< " \\___/"
351 <<endl<< " v2"
352 <<endl<< ""
353 <<endl<< " v4 = Op*v3-v2"
354 <<endl<< ""
355 <<endl<< " norm(v4)/norm(v2) <= forward_residual_solve_tol() + forward_residual_slack_error_tol()"
356 <<endl;
357
358 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
359
360 OSTab tab2(testOut);
361
362 *testOut <<endl<< "Random vector tests = " << rand_vec_i << endl;
363
364 tab.incrTab();
365
366 *testOut <<endl<< "v1 = randomize(-1,+1); ...\n" ;
367 Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,l_num_rhs);
368 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.ptr() );
369 if(dump_all()) *testOut <<endl<< "v1 =\n" << describe(*v1,verbLevel);
370
371 *testOut <<endl<< "v2 = Op*v1 ...\n" ;
372 Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,l_num_rhs);
373 timer.start(true);
374 Thyra::apply(op, NOTRANS, *v1, v2.ptr());
375 timer.stop();
376 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
377 if(dump_all()) *testOut <<endl<< "v2 =\n" << describe(*v2,verbLevel);
378
379 *testOut <<endl<< "v3 = inv(Op)*v2 ...\n" ;
380 Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,l_num_rhs);
381 SolveCriteria<Scalar> solveCriteria(
382 SolveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
383 ,forward_residual_solve_tol()
384 );
385 assign(v3.ptr(), ST::zero());
386 SolveStatus<Scalar> solveStatus;
387 {
388 VOTS lowsTempState(Teuchos::rcp(&op,false),testOut,verbLevel);
389 timer.start(true);
390 solveStatus = solve<Scalar>(op, NOTRANS, *v2, v3.ptr(),
391 optInArg(solveCriteria));
392 timer.stop();
393 OSTab(testOut).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
394 }
395 if(dump_all()) *testOut <<endl<< "v3 =\n" << describe(*v3,verbLevel);
396 *testOut
397 <<endl<< "solve status:\n";
398 OSTab(testOut).o() << solveStatus;
399 result = solveStatus.solveStatus==SOLVE_STATUS_CONVERGED;
400 if(!result) these_results = false;
401 *testOut
402 <<endl<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
403 << passfail(result)<<endl;
404
405 *testOut <<endl<< "v4 = Op*v3 - v2 ...\n" ;
406 Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(range,l_num_rhs);
407 assign( v4.ptr(), *v2 );
408 timer.start(true);
409 Thyra::apply(op, NOTRANS, *v3, v4.ptr(), as<Scalar>(1.0), as<Scalar>(-1.0));
410 timer.stop();
411 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
412 if(dump_all()) *testOut <<endl<< "v4 =\n" << describe(*v4,verbLevel);
413
414 Array<ScalarMag> norms_v2(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v2(l_num_rhs);
415 norms(*v2, norms_v2());
416 norms(*v4, norms_v4());
417 std::transform(
418 norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.begin()
419 ,std::divides<ScalarMag>()
420 );
421
422 result = testMaxErrors<Scalar>(
423 "norm(v4)/norm(v2)", norms_v4_norms_v2(),
424 "forward_residual_solve_tol()+forward_residual_slack_error_tol()",
425 as<ScalarMag>(forward_residual_solve_tol()+forward_residual_slack_error_tol()),
426 "forward_residual_solve_tol()_slack_warning_tol()",
427 as<ScalarMag>(forward_residual_solve_tol()+forward_residual_slack_warning_tol()),
428 testOut.ptr()
429 );
430 if(!result) these_results = false;
431
432 }
433 }
434 else {
435 *testOut <<endl<< "Forward operator not supported, skipping check!\n";
436 }
437
438 testResultsPrinter.printTestResults(these_results, inoutArg(success));
439
440 }
441 else {
442 if(out.get()) *out <<endl<< "this->check_forward_residual()==false: Skipping the check of the forward solve with a tolerance on the residual!\n";
443 }
444
445 if( check_adjoint_default() ) {
446
447 if(out.get()) *out <<endl<< "this->check_adjoint_default()==true: Checking the default adjoint solve ... ";
448
449 TestResultsPrinter testResultsPrinter(out, show_all_tests());
450 const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
451
452 bool these_results = true;
453
454 result = true;
455 TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(CONJTRANS), true, *testOut, result);
456 if (!result) these_results = false;
457
458 if(result) {
459
460 *testOut
461 <<endl<< "Checking that the adjoint default solve matches the adjoint operator:\n"
462 <<endl<< " inv(Op')*Op'*v1 == v1"
463 <<endl<< " \\____/"
464 <<endl<< " v2"
465 <<endl<< " \\_____________/"
466 <<endl<< " v3"
467 <<endl<< ""
468 <<endl<< " v4 = v3-v1"
469 <<endl<< " v5 = Op'*v3-v2"
470 <<endl<< ""
471 <<endl<< " norm(v4)/norm(v1) <= adjoint_default_solution_error_error_tol()"
472 <<endl<< " norm(v5)/norm(v2) <= adjoint_default_residual_error_tol()"
473 <<endl;
474
475 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
476
477 OSTab tab2(testOut);
478
479 *testOut <<endl<< "Random vector tests = " << rand_vec_i << endl;
480
481 tab.incrTab();
482
483 *testOut <<endl<< "v1 = randomize(-1,+1); ...\n" ;
484 Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(range,l_num_rhs);
485 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.ptr() );
486 if(dump_all()) *testOut <<endl<< "v1 =\n" << describe(*v1,verbLevel);
487
488 *testOut <<endl<< "v2 = Op'*v1 ...\n" ;
489 Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,l_num_rhs);
490 timer.start(true);
491 Thyra::apply(op, CONJTRANS, *v1, v2.ptr());
492 timer.stop();
493 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
494 if(dump_all()) *testOut <<endl<< "v2 =\n" << describe(*v2,verbLevel);
495
496 *testOut <<endl<< "v3 = inv(Op')*v2 ...\n" ;
497 Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,l_num_rhs);
498 assign(v3.ptr(), ST::zero());
499 SolveStatus<Scalar> solveStatus;
500 {
501 VOTS lowsTempState(Teuchos::rcp(&op,false),testOut,verbLevel);
502 timer.start(true);
503 solveStatus = solve<Scalar>(op, CONJTRANS, *v2, v3.ptr());
504 timer.stop();
505 OSTab(testOut).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
506 }
507 if(dump_all()) *testOut <<endl<< "v3 =\n" << describe(*v3,verbLevel);
508 *testOut
509 <<endl<< "solve status:\n";
510 OSTab(testOut).o() << solveStatus;
511
512 *testOut <<endl<< "v4 = v3 - v1 ...\n" ;
513 Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(range,l_num_rhs);
514 V_VmV( v4.ptr(), *v3, *v1 );
515 if(dump_all()) *testOut <<endl<< "v4 =\n" << describe(*v4,verbLevel);
516
517 *testOut <<endl<< "v5 = Op'*v3 - v2 ...\n" ;
518 Teuchos::RCP<MultiVectorBase<Scalar> > v5 = createMembers(domain,l_num_rhs);
519 assign( v5.ptr(), *v2 );
520 timer.start(true);
521 Thyra::apply(op, CONJTRANS, *v3, v5.ptr(), as<Scalar>(1.0), as<Scalar>(-1.0));
522 timer.stop();
523 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
524 if(dump_all()) *testOut <<endl<< "v5 =\n" << describe(*v5,verbLevel);
525
526 Array<ScalarMag> norms_v1(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v1(l_num_rhs);
527 norms(*v1, norms_v1());
528 norms(*v4, norms_v4());
529 std::transform(
530 norms_v4.begin(),norms_v4.end(),norms_v1.begin(),norms_v4_norms_v1.begin()
531 ,std::divides<ScalarMag>()
532 );
533
534 result = testMaxErrors<Scalar>(
535 "norm(v4)/norm(v1)", norms_v4_norms_v1(),
536 "adjoint_default_solution_error_error_tol()",
537 adjoint_default_solution_error_error_tol(),
538 "adjoint_default_solution_error_warning_tol()",
539 adjoint_default_solution_error_warning_tol(),
540 testOut.ptr()
541 );
542 if(!result) these_results = false;
543
544 Array<ScalarMag> norms_v2(l_num_rhs), norms_v5(l_num_rhs), norms_v5_norms_v2(l_num_rhs);
545 norms(*v2, norms_v2());
546 norms(*v5, norms_v5());
547 std::transform(
548 norms_v5.begin(), norms_v5.end(), norms_v2.begin(), norms_v5_norms_v2.begin()
549 ,std::divides<ScalarMag>()
550 );
551
552 result = testMaxErrors<Scalar>(
553 "norm(v5)/norm(v2)", norms_v5_norms_v2(),
554 "adjoint_default_residual_error_tol()",
555 adjoint_default_residual_error_tol(),
556 "adjoint_default_residual_warning_tol()",
557 adjoint_default_residual_warning_tol(),
558 testOut.ptr()
559 );
560 if(!result) these_results = false;
561
562 }
563 }
564 else {
565 *testOut <<endl<< "Adjoint operator not supported, skipping check!\n";
566 }
567
568 testResultsPrinter.printTestResults(these_results, inoutArg(success));
569
570 }
571 else {
572 if(out.get()) *out <<endl<< "this->check_adjoint_default()==false: Skipping the check of the adjoint solve with a default tolerance!\n";
573 }
574
575 if( check_adjoint_residual() ) {
576
577 if(out.get()) *out <<endl<< "this->check_adjoint_residual()==true: Checking the adjoint solve with a tolerance on the residual ... ";
578
579 TestResultsPrinter testResultsPrinter(out, show_all_tests());
580 const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
581
582 bool these_results = true;
583
584 result = true;
585 TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(CONJTRANS), true, *testOut, result);
586 if (!result) these_results = false;
587
588 if(result) {
589
590 *testOut
591 <<endl<< "Checking that the adjoint solve matches the adjoint operator to a residual tolerance:\n"
592 <<endl<< " v3 = inv(Op')*Op'*v1"
593 <<endl<< " \\____/"
594 <<endl<< " v2"
595 <<endl<< ""
596 <<endl<< " v4 = Op'*v3-v2"
597 <<endl<< ""
598 <<endl<< " norm(v4)/norm(v2) <= adjoint_residual_solve_tol() + adjoint_residual_slack_error_tol()"
599 <<endl;
600
601 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
602
603 OSTab tab2(testOut);
604
605 *testOut <<endl<< "Random vector tests = " << rand_vec_i << endl;
606
607 tab.incrTab();
608
609 *testOut <<endl<< "v1 = randomize(-1,+1); ...\n" ;
610 Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(range,l_num_rhs);
611 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.ptr() );
612 if(dump_all()) *testOut <<endl<< "v1 =\n" << describe(*v1,verbLevel);
613
614 *testOut <<endl<< "v2 = Op'*v1 ...\n" ;
615 Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,l_num_rhs);
616 timer.start(true);
617 Thyra::apply(op, CONJTRANS, *v1, v2.ptr());
618 timer.stop();
619 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
620 if(dump_all()) *testOut <<endl<< "v2 =\n" << describe(*v2,verbLevel);
621
622 *testOut <<endl<< "v3 = inv(Op')*v2 ...\n" ;
623 Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,l_num_rhs);
624 SolveCriteria<Scalar> solveCriteria(
625 SolveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
626 ,adjoint_residual_solve_tol()
627 );
628 assign(v3.ptr(), ST::zero());
629 SolveStatus<Scalar> solveStatus;
630 {
631 VOTS lowsTempState(Teuchos::rcp(&op,false),testOut,verbLevel);
632 timer.start(true);
633 solveStatus = solve<Scalar>(op, CONJTRANS, *v2, v3.ptr(), optInArg(solveCriteria));
634 timer.stop();
635 OSTab(testOut).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
636 }
637 if(dump_all()) *testOut <<endl<< "v3 =\n" << describe(*v3,verbLevel);
638 *testOut
639 <<endl<< "solve status:\n";
640 OSTab(testOut).o() << solveStatus;
641 result = solveStatus.solveStatus==SOLVE_STATUS_CONVERGED;
642 if(!result) these_results = false;
643 *testOut
644 <<endl<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
645 << passfail(result)<<endl;
646 if(solveStatus.achievedTol==SolveStatus<Scalar>::unknownTolerance())
647 *testOut <<endl<<"achievedTol==unknownTolerance(): Setting achievedTol = adjoint_residual_solve_tol() = "<<adjoint_residual_solve_tol()<<endl;
648
649 *testOut <<endl<< "v4 = Op'*v3 - v2 ...\n" ;
650 Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,l_num_rhs);
651 assign( v4.ptr(), *v2 );
652 timer.start(true);
653 Thyra::apply(op, CONJTRANS, *v3, v4.ptr(), as<Scalar>(1.0), as<Scalar>(-1.0));
654 timer.stop();
655 OSTab(testOut).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
656 if(dump_all()) *testOut <<endl<< "v4 =\n" << describe(*v4,verbLevel);
657
658 Array<ScalarMag> norms_v2(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v2(l_num_rhs);
659 norms(*v2, norms_v2());
660 norms(*v4, norms_v4());
661 std::transform(
662 norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.begin()
663 ,std::divides<ScalarMag>()
664 );
665
666 result = testMaxErrors<Scalar>(
667 "norm(v4)/norm(v2)", norms_v4_norms_v2(),
668 "adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol()",
669 as<ScalarMag>(adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol()),
670 "adjoint_residual_solve_tol()_slack_warning_tol()",
671 as<ScalarMag>(adjoint_residual_solve_tol()+adjoint_residual_slack_warning_tol()),
672 testOut.ptr()
673 );
674 if(!result) these_results = false;
675
676 }
677 }
678 else {
679 *testOut <<endl<< "Adjoint operator not supported, skipping check!\n";
680 }
681
682 testResultsPrinter.printTestResults(these_results, inoutArg(success));
683
684 }
685 else {
686 if(out.get())
687 *out <<endl<< "this->check_adjoint_residual()==false: Skipping the check of the adjoint solve with a tolerance on the residual!\n";
688 }
689
690 if(out.get()) {
691 if(success)
692 *out <<endl<<"Congratulations, this LinearOpWithSolveBase object seems to check out!\n";
693 else
694 *out <<endl<<"Oh no, at least one of the tests performed with this LinearOpWithSolveBase object failed (see above failures)!\n";
695 *out <<endl<< "*** Leaving LinearOpWithSolveTester<"<<ST::name()<<">::check(...)\n";
696 }
697
698 return success;
699
700}
701
702
703// private static members
704
705
706// Define local macros used in the next few lines and then undefined
707
708#define LOWST_DEFINE_RAW_STATIC_MEMBER( DATA_TYPE, MEMBER_NAME, DEFAULT_VALUE ) \
709template<class Scalar> \
710const DATA_TYPE \
711LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
712
713#define LOWST_DEFINE_MTD_STATIC_MEMBER( DATA_TYPE, MEMBER_NAME, DEFAULT_VALUE ) \
714template<class Scalar> \
715const typename LinearOpWithSolveTester<Scalar>::DATA_TYPE \
716LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
717
718LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_default_default_, true);
719LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_residual_default_, true);
720LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_default_default_, true);
721LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_residual_default_, true);
722
723LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, warning_tol_default_, 1e-6);
724LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, error_tol_default_, 1e-5);
725LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, solve_tol_default_, 1e-5);
726LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_warning_tol_default_, 1e-6);
727LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_error_tol_default_, 1e-5);
728
729LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_random_vectors_default_, 1);
730LOWST_DEFINE_RAW_STATIC_MEMBER(bool, show_all_tests_default_, false);
731LOWST_DEFINE_RAW_STATIC_MEMBER(bool, dump_all_default_, false);
732LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_rhs_default_, 1);
733
734LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSolveTol_name_, "All Solve Tol");
735LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackWarningTol_name_, "All Slack Warning Tol");
736LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackErrorTol_name_, "All Slack Error Tol");
737LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, ShowAllTests_name_, "Show All Tests");
738LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, DumpAll_name_, "Dump All");
739
740#undef LOWST_DEFINE_MTD_STATIC_MEMBER
741
742#undef LOWST_DEFINE_RAW_STATIC_MEMBER
743
744
745} // namespace Thyra
746
747
748#endif // THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
Teuchos::Array::begin
iterator begin()
Teuchos::Describable::description
virtual std::string description() const
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos::RCP::get
T * get() const
Teuchos::Time::totalElapsedTime
double totalElapsedTime(bool readCurrentTime=false) const
Teuchos::Time::start
void start(bool reset=false)
Teuchos::Time::stop
double stop()
Teuchos::basic_FancyOStream
Teuchos::basic_OSTab
Thyra::LinearOpBase
Base class for all linear operators.
Definition Thyra_LinearOpBase_decl.hpp:159
Thyra::LinearOpBase::range
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
Thyra::LinearOpBase::domain
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Thyra::LinearOpWithSolveBase
Base class for all linear operators that can support a high-level solve operation.
Definition Thyra_LinearOpWithSolveBase_decl.hpp:279
Thyra::LinearOpWithSolveBase::solveSupports
bool solveSupports(EOpTransp transp) const
Return if solve() supports the argument transp.
Definition Thyra_LinearOpWithSolveBase_decl.hpp:350
Thyra::LinearOpWithSolveTester::setParameterList
void setParameterList(const RCP< ParameterList > &paramList)
Definition Thyra_LinearOpWithSolveTester_def.hpp:112
Thyra::LinearOpWithSolveTester::set_all_slack_warning_tol
void set_all_slack_warning_tol(const ScalarMag slack_warning_tol)
Set all the warning tolerances to the same value.
Definition Thyra_LinearOpWithSolveTester_def.hpp:82
Thyra::LinearOpWithSolveTester::set_all_solve_tol
void set_all_solve_tol(const ScalarMag solve_tol)
Set all the solve tolerances to the same value.
Definition Thyra_LinearOpWithSolveTester_def.hpp:71
Thyra::LinearOpWithSolveTester::ScalarMag
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
Definition Thyra_LinearOpWithSolveTester_decl.hpp:69
Thyra::LinearOpWithSolveTester::check
bool check(const LinearOpWithSolveBase< Scalar > &op, Teuchos::FancyOStream *out) const
Check a LinearOpWithSolveBase object.
Definition Thyra_LinearOpWithSolveTester_def.hpp:160
Thyra::LinearOpWithSolveTester::turn_off_all_tests
void turn_off_all_tests()
Turn off all tests so that individual tests can be set.
Definition Thyra_LinearOpWithSolveTester_def.hpp:60
Thyra::LinearOpWithSolveTester::set_all_slack_error_tol
void set_all_slack_error_tol(const ScalarMag slack_error_tol)
Set all the error tolerances to the same value.
Definition Thyra_LinearOpWithSolveTester_def.hpp:96
Thyra::LinearOpWithSolveTester::getValidParameters
RCP< const ParameterList > getValidParameters() const
Definition Thyra_LinearOpWithSolveTester_def.hpp:130
Thyra::TestResultsPrinter
Control printing of test results.
Definition Thyra_TestingToolsDecl.hpp:207
Thyra::TestResultsPrinter::getTestOStream
RCP< FancyOStream > getTestOStream()
Return the stream used for testing.
Definition Thyra_TestingTools.cpp:102
Thyra::TestResultsPrinter::printTestResults
void printTestResults(const bool this_result, const Ptr< bool > &success)
Print the test result.
Definition Thyra_TestingTools.cpp:108
Thyra::SOLVE_STATUS_CONVERGED
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
Definition Thyra_SolveSupportTypes.hpp:361
Thyra::SOLVE_MEASURE_NORM_RESIDUAL
@ SOLVE_MEASURE_NORM_RESIDUAL
Norm of the current residual vector (i.e. ||A*x-b||)
Definition Thyra_SolveSupportTypes.hpp:32
Thyra::SOLVE_MEASURE_NORM_RHS
@ SOLVE_MEASURE_NORM_RHS
Norm of the right-hand side (i.e. ||b||)
Definition Thyra_SolveSupportTypes.hpp:38
Thyra::passfail
const std::string passfail(const bool result)
Definition Thyra_TestingToolsDecl.hpp:24
Thyra::toString
const char * toString(EConj conj)
Return a string name for a EOpTransp value. `*.
Definition Thyra_OperatorVectorTypes.hpp:113
Thyra::NOTRANS
@ NOTRANS
Use the non-transposed operator.
Definition Thyra_OperatorVectorTypes.hpp:130
Thyra::CONJTRANS
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
Definition Thyra_OperatorVectorTypes.hpp:140
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::VERB_MEDIUM
VERB_MEDIUM
Teuchos::VERB_EXTREME
VERB_EXTREME
TEUCHOS_TEST_EQUALITY_CONST
#define TEUCHOS_TEST_EQUALITY_CONST(v1, v2, out, success)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Thyra::SolveCriteria
Simple struct that defines the requested solution criteria for a solve.
Definition Thyra_SolveSupportTypes.hpp:280
Thyra::SolveStatus
Simple struct for the return status from a solve.
Definition Thyra_SolveSupportTypes.hpp:391
Thyra::SolveStatus::solveStatus
ESolveStatus solveStatus
The return status of the solve.
Definition Thyra_SolveSupportTypes.hpp:397
Thyra::SolveStatus::achievedTol
ScalarMag achievedTol
The maximum final tolerance actually achieved by the (block) linear solve. A value of unknownToleranc...
Definition Thyra_SolveSupportTypes.hpp:401
Generated by doxygen 1.9.8