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 //
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 
59 namespace Thyra {
60 
61 
62 // Constructors/initializers
63 
64 
65 template<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 
92 template<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 
102 template<class Scalar>
103 void
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 
113 template<class Scalar>
114 void
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 
127 template<class Scalar>
128 void
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 
144 template<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 
161 template<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 
192 template<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;
202  using Teuchos::FancyOStream;
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;
679  if(solveStatus.achievedTol==SolveStatus<Scalar>::unknownTolerance())
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 ) \
742 template<class Scalar> \
743 const DATA_TYPE \
744 LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
745 
746 #define LOWST_DEFINE_MTD_STATIC_MEMBER( DATA_TYPE, MEMBER_NAME, DEFAULT_VALUE ) \
747 template<class Scalar> \
748 const typename LinearOpWithSolveTester<Scalar>::DATA_TYPE \
749 LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
750 
751 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_default_default_, true);
752 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_residual_default_, true);
753 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_default_default_, true);
754 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_residual_default_, true);
755 
756 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, warning_tol_default_, 1e-6);
757 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, error_tol_default_, 1e-5);
758 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, solve_tol_default_, 1e-5);
759 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_warning_tol_default_, 1e-6);
760 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_error_tol_default_, 1e-5);
761 
762 LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_random_vectors_default_, 1);
763 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, show_all_tests_default_, false);
764 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, dump_all_default_, false);
765 LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_rhs_default_, 1);
766 
767 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSolveTol_name_, "All Solve Tol");
768 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackWarningTol_name_, "All Slack Warning Tol");
769 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackErrorTol_name_, "All Slack Error Tol");
770 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, ShowAllTests_name_, "Show All Tests");
771 LOWST_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
Control printing of test results.
void turn_off_all_tests()
Turn off all tests so that individual tests can be set.
void set_all_slack_warning_tol(const ScalarMag slack_warning_tol)
Set all the warning tolerances to the same value.
Base class for all linear operators that can support a high-level solve operation.
bool is_null(const boost::shared_ptr< T > &p)
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
basic_OSTab< char > OSTab
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
basic_FancyOStream< char > FancyOStream
const std::string passfail(const bool result)
Use the non-transposed operator.
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
T * get() const
RCP< const ParameterList > getValidParameters() const
bool check(const LinearOpWithSolveBase< Scalar > &op, Teuchos::FancyOStream *out) const
Check a LinearOpWithSolveBase object.
void set_all_solve_tol(const ScalarMag solve_tol)
Set all the solve tolerances to the same value.
void start(bool reset=false)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
double stop()
Ptr< T > ptr() const
void set_all_slack_error_tol(const ScalarMag slack_error_tol)
Set all the error tolerances to the same value.
void printTestResults(const bool this_result, const Ptr< bool > &success)
Print the test result.
virtual std::string description() const
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
Simple struct for the return status from a solve.
Base class for all linear operators.
Norm of the right-hand side (i.e. ||b||)
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
bool solveSupports(EOpTransp transp) const
Return if solve() supports the argument transp.
TypeTo as(const TypeFrom &t)
double totalElapsedTime(bool readCurrentTime=false) const
The requested solution criteria has likely been achieved.
void setParameterList(const RCP< ParameterList > &paramList)
Simple struct that defines the requested solution criteria for a solve.
RCP< FancyOStream > getTestOStream()
Return the stream used for testing.
iterator begin()
#define TEUCHOS_TEST_EQUALITY_CONST(v1, v2, out, success)
Norm of the current residual vector (i.e. ||A*x-b||)