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