Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
26 namespace Thyra {
27 
28 
29 // Constructors/initializers
30 
31 
32 template<class Scalar>
33 LinearOpWithSolveTester<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 
59 template<class Scalar>
60 void 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 
69 template<class Scalar>
70 void
71 LinearOpWithSolveTester<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 
80 template<class Scalar>
81 void
82 LinearOpWithSolveTester<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 
94 template<class Scalar>
95 void
96 LinearOpWithSolveTester<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 
111 template<class Scalar>
112 void 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 
128 template<class Scalar>
129 RCP<const ParameterList>
130 LinearOpWithSolveTester<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 
159 template<class Scalar>
160 bool 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 ) \
709 template<class Scalar> \
710 const DATA_TYPE \
711 LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
712 
713 #define LOWST_DEFINE_MTD_STATIC_MEMBER( DATA_TYPE, MEMBER_NAME, DEFAULT_VALUE ) \
714 template<class Scalar> \
715 const typename LinearOpWithSolveTester<Scalar>::DATA_TYPE \
716 LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
717 
718 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_default_default_, true);
719 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_residual_default_, true);
720 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_default_default_, true);
721 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_residual_default_, true);
722 
723 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, warning_tol_default_, 1e-6);
724 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, error_tol_default_, 1e-5);
725 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, solve_tol_default_, 1e-5);
726 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_warning_tol_default_, 1e-6);
727 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_error_tol_default_, 1e-5);
728 
729 LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_random_vectors_default_, 1);
730 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, show_all_tests_default_, false);
731 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, dump_all_default_, false);
732 LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_rhs_default_, 1);
733 
734 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSolveTol_name_, "All Solve Tol");
735 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackWarningTol_name_, "All Slack Warning Tol");
736 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackErrorTol_name_, "All Slack Error Tol");
737 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, ShowAllTests_name_, "Show All Tests");
738 LOWST_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
Thyra::TestResultsPrinter
Control printing of test results.
Definition: Thyra_TestingToolsDecl.hpp:207
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_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::LinearOpWithSolveBase
Base class for all linear operators that can support a high-level solve operation.
Definition: Thyra_LinearOpWithSolveBase_decl.hpp:276
is_null
bool is_null(const boost::shared_ptr< T > &p)
Thyra::LinearOpBase::range
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
Teuchos::OSTab
basic_OSTab< char > OSTab
Thyra::LinearOpWithSolveTester::ScalarMag
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
Definition: Thyra_LinearOpWithSolveTester_decl.hpp:69
Teuchos::FancyOStream
basic_FancyOStream< char > FancyOStream
passfail
const std::string passfail(const bool result)
Thyra::NOTRANS
Use the non-transposed operator.
Definition: Thyra_OperatorVectorTypes.hpp:130
Teuchos::VERB_EXTREME
Thyra::CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
Definition: Thyra_OperatorVectorTypes.hpp:140
Teuchos::RCP::get
T * get() const
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Thyra::LinearOpWithSolveTester::getValidParameters
RCP< const ParameterList > getValidParameters() const
Definition: Thyra_LinearOpWithSolveTester_def.hpp:130
Teuchos::VERB_MEDIUM
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::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
Teuchos::Time::start
void start(bool reset=false)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::Time::stop
double stop()
Teuchos::basic_FancyOStream
Teuchos::RCP::ptr
Ptr< T > ptr() const
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::TestResultsPrinter::printTestResults
void printTestResults(const bool this_result, const Ptr< bool > &success)
Print the test result.
Definition: Thyra_TestingTools.cpp:108
Teuchos::Describable::description
virtual std::string description() const
toString
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
Teuchos::ParameterList::validateParametersAndSetDefaults
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Thyra::SolveStatus
Simple struct for the return status from a solve.
Definition: Thyra_SolveSupportTypes.hpp:391
Thyra::LinearOpBase
Base class for all linear operators.
Definition: Thyra_LinearOpBase_decl.hpp:159
Thyra::SOLVE_MEASURE_NORM_RHS
Norm of the right-hand side (i.e. ||b||)
Definition: Thyra_SolveSupportTypes.hpp:38
Thyra::LinearOpBase::domain
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Thyra::LinearOpWithSolveBase::solveSupports
bool solveSupports(EOpTransp transp) const
Return if solve() supports the argument transp.
Definition: Thyra_LinearOpWithSolveBase_decl.hpp:348
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::Time::totalElapsedTime
double totalElapsedTime(bool readCurrentTime=false) const
Thyra::SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
Definition: Thyra_SolveSupportTypes.hpp:361
Thyra::LinearOpWithSolveTester::setParameterList
void setParameterList(const RCP< ParameterList > &paramList)
Definition: Thyra_LinearOpWithSolveTester_def.hpp:112
Thyra::SolveCriteria
Simple struct that defines the requested solution criteria for a solve.
Definition: Thyra_SolveSupportTypes.hpp:280
Thyra::TestResultsPrinter::getTestOStream
RCP< FancyOStream > getTestOStream()
Return the stream used for testing.
Definition: Thyra_TestingTools.cpp:102
Teuchos::Array::begin
iterator begin()
TEUCHOS_TEST_EQUALITY_CONST
#define TEUCHOS_TEST_EQUALITY_CONST(v1, v2, out, success)
Thyra::SOLVE_MEASURE_NORM_RESIDUAL
Norm of the current residual vector (i.e. ||A*x-b||)
Definition: Thyra_SolveSupportTypes.hpp:32
Generated on Thu Nov 21 2024 09:20:52 for Thyra by   doxygen 1.8.5