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_)
200 using Teuchos::optInArg;
201 using Teuchos::inoutArg;
207 bool success =
true, result;
208 const int l_num_rhs = this->num_rhs();
213 OSTab tab(out,1,
"THYRA");
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);
220 *out <<endl<<
"describe adjoint op:\n";
221 describeLinearOp<Scalar>(
228 *out <<endl<<
"describe op: " << op.
description() << endl;
235 if( check_forward_default() ) {
237 if(out.
get()) *out <<endl<<
"this->check_forward_default()==true: Checking the default forward solve ... ";
242 bool these_results =
true;
246 if (!result) these_results =
false;
251 <<endl<<
"Checking that the forward default solve matches the forward operator:\n"
252 <<endl<<
" inv(Op)*Op*v1 == v1"
255 <<endl<<
" \\___________/"
258 <<endl<<
" v4 = v3-v1"
259 <<endl<<
" v5 = Op*v3-v2"
261 <<endl<<
" norm(v4)/norm(v1) <= forward_default_solution_error_error_tol()"
262 <<endl<<
" norm(v5)/norm(v2) <= forward_default_residual_error_tol()"
265 for(
int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
269 *testOut <<endl<<
"Random vector tests = " << rand_vec_i << endl;
273 *testOut <<endl<<
"v1 = randomize(-1,+1); ...\n" ;
275 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.
ptr() );
276 if(dump_all()) *testOut <<endl<<
"v1 =\n" << describe(*v1,verbLevel);
278 *testOut <<endl<<
"v2 = Op*v1 ...\n" ;
283 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
284 if(dump_all()) *testOut <<endl<<
"v2 =\n" << describe(*v2,verbLevel);
286 *testOut <<endl<<
"v3 = inv(Op)*v2 ...\n" ;
288 assign(v3.
ptr(), ST::zero());
291 VOTS lowsTempState(
Teuchos::rcp(&op,
false),testOut,verbLevel);
293 solveStatus = solve<Scalar>(op,
NOTRANS, *v2, v3.
ptr());
295 OSTab(testOut).o() <<
"\n=> Solve time = "<<timer.
totalElapsedTime()<<
" sec\n";
297 if(dump_all()) *testOut <<endl<<
"v3 =\n" << describe(*v3,verbLevel);
299 <<endl<<
"solve status:\n";
300 OSTab(testOut).o() << solveStatus;
302 *testOut <<endl<<
"v4 = v3 - v1 ...\n" ;
304 V_VmV( v4.
ptr(), *v3, *v1 );
305 if(dump_all()) *testOut <<endl<<
"v4 =\n" << describe(*v4,verbLevel);
307 *testOut <<endl<<
"v5 = Op*v3 - v2 ...\n" ;
309 assign( v5.
ptr(), *v2 );
311 Thyra::apply(op,
NOTRANS, *v3, v5.
ptr(), as<Scalar>(1.0) ,as<Scalar>(-1.0));
313 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
314 if(dump_all()) *testOut <<endl<<
"v5 =\n" << describe(*v5,verbLevel);
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());
320 norms_v4.begin(), norms_v4.end(), norms_v1.begin(), norms_v4_norms_v1.
begin()
321 ,std::divides<ScalarMag>()
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(),
332 if(!result) these_results =
false;
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());
338 norms_v5.begin(), norms_v5.end(), norms_v2.begin(), norms_v5_norms_v2.
begin()
339 ,std::divides<ScalarMag>()
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(),
350 if(!result) these_results =
false;
355 *testOut <<endl<<
"Forward operator not supported, skipping check!\n";
362 if(out.
get()) *out <<endl<<
"this->check_forward_default()==false: Skipping the check of the default forward solve!\n";
365 if( check_forward_residual() ) {
367 if(out.
get()) *out <<endl<<
"this->check_forward_residual()==true: Checking the forward solve with a tolerance on the residual ... ";
372 bool these_results =
true;
376 if (!result) these_results =
false;
381 <<endl<<
"Checking that the forward solve matches the forward operator to a residual tolerance:\n"
382 <<endl<<
" v3 = inv(Op)*Op*v1"
386 <<endl<<
" v4 = Op*v3-v2"
388 <<endl<<
" norm(v4)/norm(v2) <= forward_residual_solve_tol() + forward_residual_slack_error_tol()"
391 for(
int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
395 *testOut <<endl<<
"Random vector tests = " << rand_vec_i << endl;
399 *testOut <<endl<<
"v1 = randomize(-1,+1); ...\n" ;
401 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.
ptr() );
402 if(dump_all()) *testOut <<endl<<
"v1 =\n" << describe(*v1,verbLevel);
404 *testOut <<endl<<
"v2 = Op*v1 ...\n" ;
409 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
410 if(dump_all()) *testOut <<endl<<
"v2 =\n" << describe(*v2,verbLevel);
412 *testOut <<endl<<
"v3 = inv(Op)*v2 ...\n" ;
416 ,forward_residual_solve_tol()
418 assign(v3.
ptr(), ST::zero());
421 VOTS lowsTempState(
Teuchos::rcp(&op,
false),testOut,verbLevel);
423 solveStatus = solve<Scalar>(op,
NOTRANS, *v2, v3.
ptr(),
424 optInArg(solveCriteria));
426 OSTab(testOut).o() <<
"\n=> Solve time = "<<timer.
totalElapsedTime()<<
" sec\n";
428 if(dump_all()) *testOut <<endl<<
"v3 =\n" << describe(*v3,verbLevel);
430 <<endl<<
"solve status:\n";
431 OSTab(testOut).o() << solveStatus;
433 if(!result) these_results =
false;
435 <<endl<<
"check: solveStatus = " <<
toString(solveStatus.
solveStatus) <<
" == SOLVE_STATUS_CONVERGED : "
438 *testOut <<endl<<
"v4 = Op*v3 - v2 ...\n" ;
440 assign( v4.
ptr(), *v2 );
442 Thyra::apply(op,
NOTRANS, *v3, v4.
ptr(), as<Scalar>(1.0), as<Scalar>(-1.0));
444 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
445 if(dump_all()) *testOut <<endl<<
"v4 =\n" << describe(*v4,verbLevel);
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());
451 norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.
begin()
452 ,std::divides<ScalarMag>()
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()),
463 if(!result) these_results =
false;
468 *testOut <<endl<<
"Forward operator not supported, skipping check!\n";
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";
478 if( check_adjoint_default() ) {
480 if(out.
get()) *out <<endl<<
"this->check_adjoint_default()==true: Checking the default adjoint solve ... ";
485 bool these_results =
true;
489 if (!result) these_results =
false;
494 <<endl<<
"Checking that the adjoint default solve matches the adjoint operator:\n"
495 <<endl<<
" inv(Op')*Op'*v1 == v1"
498 <<endl<<
" \\_____________/"
501 <<endl<<
" v4 = v3-v1"
502 <<endl<<
" v5 = Op'*v3-v2"
504 <<endl<<
" norm(v4)/norm(v1) <= adjoint_default_solution_error_error_tol()"
505 <<endl<<
" norm(v5)/norm(v2) <= adjoint_default_residual_error_tol()"
508 for(
int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
512 *testOut <<endl<<
"Random vector tests = " << rand_vec_i << endl;
516 *testOut <<endl<<
"v1 = randomize(-1,+1); ...\n" ;
518 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.
ptr() );
519 if(dump_all()) *testOut <<endl<<
"v1 =\n" << describe(*v1,verbLevel);
521 *testOut <<endl<<
"v2 = Op'*v1 ...\n" ;
526 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
527 if(dump_all()) *testOut <<endl<<
"v2 =\n" << describe(*v2,verbLevel);
529 *testOut <<endl<<
"v3 = inv(Op')*v2 ...\n" ;
531 assign(v3.
ptr(), ST::zero());
534 VOTS lowsTempState(
Teuchos::rcp(&op,
false),testOut,verbLevel);
536 solveStatus = solve<Scalar>(op,
CONJTRANS, *v2, v3.
ptr());
538 OSTab(testOut).o() <<
"\n=> Solve time = "<<timer.
totalElapsedTime()<<
" sec\n";
540 if(dump_all()) *testOut <<endl<<
"v3 =\n" << describe(*v3,verbLevel);
542 <<endl<<
"solve status:\n";
543 OSTab(testOut).o() << solveStatus;
545 *testOut <<endl<<
"v4 = v3 - v1 ...\n" ;
547 V_VmV( v4.
ptr(), *v3, *v1 );
548 if(dump_all()) *testOut <<endl<<
"v4 =\n" << describe(*v4,verbLevel);
550 *testOut <<endl<<
"v5 = Op'*v3 - v2 ...\n" ;
552 assign( v5.
ptr(), *v2 );
554 Thyra::apply(op,
CONJTRANS, *v3, v5.
ptr(), as<Scalar>(1.0), as<Scalar>(-1.0));
556 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
557 if(dump_all()) *testOut <<endl<<
"v5 =\n" << describe(*v5,verbLevel);
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());
563 norms_v4.begin(),norms_v4.end(),norms_v1.begin(),norms_v4_norms_v1.
begin()
564 ,std::divides<ScalarMag>()
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(),
575 if(!result) these_results =
false;
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());
581 norms_v5.begin(), norms_v5.end(), norms_v2.begin(), norms_v5_norms_v2.
begin()
582 ,std::divides<ScalarMag>()
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(),
593 if(!result) these_results =
false;
598 *testOut <<endl<<
"Adjoint operator not supported, skipping check!\n";
605 if(out.
get()) *out <<endl<<
"this->check_adjoint_default()==false: Skipping the check of the adjoint solve with a default tolerance!\n";
608 if( check_adjoint_residual() ) {
610 if(out.
get()) *out <<endl<<
"this->check_adjoint_residual()==true: Checking the adjoint solve with a tolerance on the residual ... ";
615 bool these_results =
true;
619 if (!result) these_results =
false;
624 <<endl<<
"Checking that the adjoint solve matches the adjoint operator to a residual tolerance:\n"
625 <<endl<<
" v3 = inv(Op')*Op'*v1"
629 <<endl<<
" v4 = Op'*v3-v2"
631 <<endl<<
" norm(v4)/norm(v2) <= adjoint_residual_solve_tol() + adjoint_residual_slack_error_tol()"
634 for(
int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
638 *testOut <<endl<<
"Random vector tests = " << rand_vec_i << endl;
642 *testOut <<endl<<
"v1 = randomize(-1,+1); ...\n" ;
644 Thyra::randomize( as<Scalar>(-1.0), as<Scalar>(+1.0), v1.
ptr() );
645 if(dump_all()) *testOut <<endl<<
"v1 =\n" << describe(*v1,verbLevel);
647 *testOut <<endl<<
"v2 = Op'*v1 ...\n" ;
652 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
653 if(dump_all()) *testOut <<endl<<
"v2 =\n" << describe(*v2,verbLevel);
655 *testOut <<endl<<
"v3 = inv(Op')*v2 ...\n" ;
659 ,adjoint_residual_solve_tol()
661 assign(v3.
ptr(), ST::zero());
664 VOTS lowsTempState(
Teuchos::rcp(&op,
false),testOut,verbLevel);
666 solveStatus = solve<Scalar>(op,
CONJTRANS, *v2, v3.
ptr(), optInArg(solveCriteria));
668 OSTab(testOut).o() <<
"\n=> Solve time = "<<timer.
totalElapsedTime()<<
" sec\n";
670 if(dump_all()) *testOut <<endl<<
"v3 =\n" << describe(*v3,verbLevel);
672 <<endl<<
"solve status:\n";
673 OSTab(testOut).o() << solveStatus;
675 if(!result) these_results =
false;
677 <<endl<<
"check: solveStatus = " <<
toString(solveStatus.
solveStatus) <<
" == SOLVE_STATUS_CONVERGED : "
680 *testOut <<endl<<
"achievedTol==unknownTolerance(): Setting achievedTol = adjoint_residual_solve_tol() = "<<adjoint_residual_solve_tol()<<endl;
682 *testOut <<endl<<
"v4 = Op'*v3 - v2 ...\n" ;
684 assign( v4.
ptr(), *v2 );
686 Thyra::apply(op,
CONJTRANS, *v3, v4.
ptr(), as<Scalar>(1.0), as<Scalar>(-1.0));
688 OSTab(testOut).o() <<
"\n=> Apply time = "<<timer.
totalElapsedTime()<<
" sec\n";
689 if(dump_all()) *testOut <<endl<<
"v4 =\n" << describe(*v4,verbLevel);
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());
695 norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.
begin()
696 ,std::divides<ScalarMag>()
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()),
707 if(!result) these_results =
false;
712 *testOut <<endl<<
"Adjoint operator not supported, skipping check!\n";
720 *out <<endl<<
"this->check_adjoint_residual()==false: Skipping the check of the adjoint solve with a tolerance on the residual!\n";
725 *out <<endl<<
"Congratulations, this LinearOpWithSolveBase object seems to check out!\n";
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";