Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_LinearOpTester_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_TESTER_DEF_HPP
11 #define THYRA_LINEAR_OP_TESTER_DEF_HPP
12 
13 #include "Thyra_LinearOpTester_decl.hpp"
14 #include "Thyra_LinearOpBase.hpp"
15 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
16 #include "Thyra_describeLinearOp.hpp"
17 #include "Thyra_VectorStdOps.hpp"
18 #include "Thyra_TestingTools.hpp"
19 #include "Thyra_UniversalMultiVectorRandomizer.hpp"
20 #include "Teuchos_TestingHelpers.hpp"
21 
22 
23 namespace Thyra {
24 
25 
26 template<class Scalar>
27 class SymmetricLinearOpTester {
28 public:
29  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
30  static void checkSymmetry(
31  const LinearOpBase<Scalar> &op,
32  const Ptr<MultiVectorRandomizerBase<Scalar> > &dRand,
33  Teuchos::FancyOStream &testOut,
34  const int num_rhs,
35  const int num_random_vectors,
36  const Teuchos::EVerbosityLevel verbLevel,
37  const bool dump_all,
38  const ScalarMag &symmetry_error_tol,
39  const ScalarMag &symmetry_warning_tol,
40  const Ptr<bool> &these_results
41  )
42  {
43 
44  using std::endl;
45 
46  bool result;
47  using Teuchos::OSTab;
48  typedef Teuchos::ScalarTraits<Scalar> ST;
49  const Scalar half = Scalar(0.4)*ST::one();
50  RCP<const VectorSpaceBase<Scalar> > domain = op.domain();
51 
52  testOut << endl << "op.domain()->isCompatible(*op.range()) == true : ";
53  result = op.domain()->isCompatible(*op.range());
54  if(!result) *these_results = false;
55  testOut << passfail(result) << endl;
56 
57  if(result) {
58 
59  testOut
60  << endl << "Checking that the operator is symmetric as:\n"
61  << endl << " <0.5*op*v2,v1> == <v2,0.5*op*v1>"
62  << endl << " \\_______/ \\_______/"
63  << endl << " v4 v3"
64  << endl << ""
65  << endl << " <v4,v1> == <v2,v3>"
66  << endl << std::flush;
67 
68  for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors; ++rand_vec_i ) {
69 
70  testOut << endl << "Random vector tests = " << rand_vec_i << endl;
71 
72  OSTab tab(testOut);
73 
74  if(dump_all) testOut << endl << "v1 = randomize(-1,+1); ...\n" ;
75  RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,num_rhs);
76  dRand->randomize(v1.ptr());
77  if(dump_all) testOut << endl << "v1 =\n" << describe(*v1,verbLevel);
78 
79  if(dump_all) testOut << endl << "v2 = randomize(-1,+1); ...\n" ;
80  RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,num_rhs);
81  dRand->randomize(v2.ptr());
82  if(dump_all) testOut << endl << "v2 =\n" << describe(*v2,verbLevel);
83 
84  if(dump_all) testOut << endl << "v3 = 0.5*op*v1 ...\n" ;
85  RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,num_rhs);
86  apply( op, NOTRANS, *v1, v3.ptr(), half );
87  if(dump_all) testOut << endl << "v3 =\n" << describe(*v3,verbLevel);
88 
89  if(dump_all) testOut << endl << "v4 = 0.5*op*v2 ...\n" ;
90  RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,num_rhs);
91  apply( op, NOTRANS, *v2, v4.ptr(), half );
92  if(dump_all) testOut << endl << "v4 =\n" << describe(*v4,verbLevel);
93 
94  Array<Scalar> prod1(num_rhs), prod2(num_rhs);
95  domain->scalarProds(*v4, *v1, prod1());
96  domain->scalarProds(*v2, *v3, prod2());
97 
98  result = testRelErrors<Scalar, Scalar, ScalarMag>(
99  "<v4,v1>", prod1(),
100  "<v2,v3>", prod2(),
101  "symmetry_error_tol()", symmetry_error_tol,
102  "symmetry_warning_tol()", symmetry_warning_tol,
103  inOutArg(testOut)
104  );
105  if(!result) *these_results = false;
106 
107  }
108  }
109  else {
110  testOut << endl << "Range and domain spaces are different, skipping check!\n";
111  }
112  }
113 };
114 
115 
116 //
117 // LinearOpTester
118 //
119 
120 
121 template<class Scalar>
122 LinearOpTester<Scalar>::LinearOpTester()
123  :check_linear_properties_(true),
124  linear_properties_warning_tol_(-1.0),
125  linear_properties_error_tol_(-1.0),
126  check_adjoint_(true),
127  adjoint_warning_tol_(-1.0),
128  adjoint_error_tol_(-1.0),
129  check_for_symmetry_(false),
130  symmetry_warning_tol_(-1.0),
131  symmetry_error_tol_(-1.0),
132  num_random_vectors_(1),
133  show_all_tests_(false),
134  dump_all_(false),
135  num_rhs_(1)
136 {
137  setDefaultTols();
138 }
139 
140 
141 template<class Scalar>
142 void LinearOpTester<Scalar>::enable_all_tests( const bool enable_all_tests_in )
143 {
144  check_linear_properties_ = enable_all_tests_in;
145  check_adjoint_ = enable_all_tests_in;
146  check_for_symmetry_ = enable_all_tests_in;
147 }
148 
149 
150 template<class Scalar>
151 void LinearOpTester<Scalar>::set_all_warning_tol( const ScalarMag warning_tol_in )
152 {
153  linear_properties_warning_tol_ = warning_tol_in;
154  adjoint_warning_tol_ = warning_tol_in;
155  symmetry_warning_tol_ = warning_tol_in;
156 }
157 
158 
159 template<class Scalar>
160 void LinearOpTester<Scalar>::set_all_error_tol( const ScalarMag error_tol_in )
161 {
162  linear_properties_error_tol_ = error_tol_in;
163  adjoint_error_tol_ = error_tol_in;
164  symmetry_error_tol_ = error_tol_in;
165 }
166 
167 
168 template<class Scalar>
169 bool LinearOpTester<Scalar>::check(
170  const LinearOpBase<Scalar> &op,
171  const Ptr<MultiVectorRandomizerBase<Scalar> > &rangeRandomizer,
172  const Ptr<MultiVectorRandomizerBase<Scalar> > &domainRandomizer,
173  const Ptr<Teuchos::FancyOStream> &out_inout
174  ) const
175 {
176 
177  using std::endl;
178  using Teuchos::as;
179  using Teuchos::rcp;
180  using Teuchos::rcpFromPtr;
181  using Teuchos::rcpFromRef;
182  using Teuchos::outArg;
183  using Teuchos::inoutArg;
184  using Teuchos::fancyOStream;
185  using Teuchos::FancyOStream;
186  using Teuchos::OSTab;
187  typedef Teuchos::ScalarTraits<Scalar> ST;
188 
189  bool success = true, result;
190  const int loc_num_rhs = this->num_rhs();
191  const Scalar r_one = ST::one();
192  const Scalar d_one = ST::one();
193  const Scalar r_half = as<Scalar>(0.5)*r_one;
194  const Scalar d_half = as<Scalar>(0.5)*d_one;
195 
196  RCP<FancyOStream> out;
197  if (!is_null(out_inout))
198  out = Teuchos::rcpFromPtr(out_inout);
199  else
200  out = Teuchos::fancyOStream(rcp(new Teuchos::oblackholestream));
201 
202  const Teuchos::EVerbosityLevel verbLevel =
203  (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
204 
205  OSTab tab2(out,1,"THYRA");
206 
207  // ToDo 04/28/2005:
208  // * Test the MultiVectorBase apply() function and output to the VectorBase apply() function!
209 
210  *out << endl << "*** Entering LinearOpTester<"<<ST::name()<<","<<ST::name()<<">::check(op,...) ...\n";
211  if(show_all_tests()) {
212  *out << endl << "describe op:\n" << Teuchos::describe(op,verbLevel);
213 /*
214  if(op.applyTransposeSupports(CONJ_ELE) && verbLevel==Teuchos::VERB_EXTREME) {
215  *out << endl << "describe adjoint op:\n";
216  describeLinearOp(*adjoint(Teuchos::rcp(&op,false)),*out,verbLevel);
217  }
218 */
219  }
220  else {
221  *out << endl << "describe op: " << Teuchos::describe(op, Teuchos::VERB_LOW);
222  }
223 
224  RCP< MultiVectorRandomizerBase<Scalar> > rRand;
225  if (!is_null(rangeRandomizer))
226  rRand = rcpFromPtr(rangeRandomizer);
227  else
228  rRand = universalMultiVectorRandomizer<Scalar>();
229  RCP< MultiVectorRandomizerBase<Scalar> > dRand;
230  if (!is_null(domainRandomizer))
231  dRand = rcpFromPtr(domainRandomizer);
232  else
233  dRand = universalMultiVectorRandomizer<Scalar>();
234 
235  *out << endl << "Checking the domain and range spaces ... ";
236 
237  RCP<const VectorSpaceBase<Scalar> > range = op.range();
238  RCP<const VectorSpaceBase<Scalar> > domain = op.domain();
239 
240  {
241 
242  TestResultsPrinter testResultsPrinter(out, show_all_tests());
243  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
244 
245  bool these_results = true;
246 
247  *testOut << endl << "op.domain().get() != NULL ? ";
248  result = domain.get() != NULL;
249  if(!result) these_results = false;
250  *testOut << passfail(result) << endl;
251 
252  *testOut << endl << "op.range().get() != NULL ? ";
253  result = range.get() != NULL;
254  if(!result) these_results = false;
255  *testOut << passfail(result) << endl;
256 
257  testResultsPrinter.printTestResults(these_results, inoutArg(success));
258 
259  }
260 
261  if( check_linear_properties() ) {
262 
263  *out << endl << "this->check_linear_properties()==true:"
264  << "Checking the linear properties of the forward linear operator ... ";
265 
266  TestResultsPrinter testResultsPrinter(out, show_all_tests());
267  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
268 
269  bool these_results = true;
270 
271  TEUCHOS_TEST_EQUALITY_CONST( op.opSupported(NOTRANS), true, *testOut, these_results );
272 
273  if(result) {
274 
275  *testOut
276  << endl << "Checking that the forward operator is truly linear:\n"
277  << endl << " 0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2"
278  << endl << " \\_____/ \\___/"
279  << endl << " v3 v5"
280  << endl << " \\_____________/ \\___________________/"
281  << endl << " v4 v5"
282  << endl << ""
283  << endl << " sum(v4) == sum(v5)"
284  << endl << std::flush;
285 
286  for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
287 
288  *testOut << endl << "Random vector tests = " << rand_vec_i << endl;
289 
290  OSTab tab3(testOut);
291 
292  *testOut << endl << "v1 = randomize(-1,+1); ...\n" ;
293  RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,loc_num_rhs);
294  dRand->randomize(v1.ptr());
295  if(dump_all()) *testOut << endl << "v1 =\n" << describe(*v1,verbLevel);
296 
297  *testOut << endl << "v2 = randomize(-1,+1); ...\n" ;
298  RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,loc_num_rhs);
299  dRand->randomize(v2.ptr());
300  if(dump_all()) *testOut << endl << "v2 =\n" << describe(*v2,verbLevel);
301 
302  *testOut << endl << "v3 = v1 + v2 ...\n" ;
303  RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,loc_num_rhs);
304  V_VpV(v3.ptr(),*v1,*v2);
305  if(dump_all()) *testOut << endl << "v3 =\n" << describe(*v3,verbLevel);
306 
307  *testOut << endl << "v4 = 0.5*op*v3 ...\n" ;
308  RCP<MultiVectorBase<Scalar> > v4 = createMembers(range,loc_num_rhs);
309  apply( op, NOTRANS, *v3, v4.ptr(), r_half );
310  if(dump_all()) *testOut << endl << "v4 =\n" << describe(*v4,verbLevel);
311 
312  *testOut << endl << "v5 = op*v1 ...\n" ;
313  RCP<MultiVectorBase<Scalar> > v5 = createMembers(range,loc_num_rhs);
314  apply( op, NOTRANS, *v1, v5.ptr() );
315  if(dump_all()) *testOut << endl << "v5 =\n" << describe(*v5,verbLevel);
316 
317  *testOut << endl << "v5 = 0.5*op*v2 + 0.5*v5 ...\n" ;
318  apply( op, NOTRANS, *v2, v5.ptr(), r_half, r_half );
319  if(dump_all()) *testOut << endl << "v5 =\n" << describe(*v5,verbLevel);
320 
321  Array<Scalar> sum_v4(loc_num_rhs), sum_v5(loc_num_rhs);
322  sums(*v4, sum_v4());
323  sums(*v5, sum_v5());
324 
325  result = testRelErrors<Scalar, Scalar, ScalarMag>(
326  "sum(v4)", sum_v4(),
327  "sum(v5)", sum_v5(),
328  "linear_properties_error_tol()", linear_properties_error_tol(),
329  "linear_properties_warning_tol()", linear_properties_warning_tol(),
330  testOut.ptr()
331  );
332  if(!result) these_results = false;
333 
334  }
335  }
336  else {
337  *testOut << endl << "Forward operator not supported, skipping check!\n";
338  }
339 
340  testResultsPrinter.printTestResults(these_results, inoutArg(success));
341 
342  }
343  else {
344  *out << endl << "this->check_linear_properties()==false: Skipping the check of the linear properties of the forward operator!\n";
345  }
346 
347  if( check_linear_properties() && check_adjoint() ) {
348 
349  *out << endl << "(this->check_linear_properties()&&this->check_adjoint())==true:"
350  << " Checking the linear properties of the adjoint operator ... ";
351 
352 
353  TestResultsPrinter testResultsPrinter(out, show_all_tests());
354  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
355 
356  bool these_results = true;
357 
358  TEUCHOS_TEST_EQUALITY_CONST( op.opSupported(CONJTRANS), true, *testOut, these_results );
359 
360  if(result) {
361 
362  *testOut
363  << endl << "Checking that the adjoint operator is truly linear:\n"
364  << endl << " 0.5*op'*(v1 + v2) == 0.5*op'*v1 + 0.5*op'*v2"
365  << endl << " \\_____/ \\____/"
366  << endl << " v3 v5"
367  << endl << " \\_______________/ \\_____________________/"
368  << endl << " v4 v5"
369  << endl << ""
370  << endl << " sum(v4) == sum(v5)"
371  << endl << std::flush;
372 
373  for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
374 
375  *testOut << endl << "Random vector tests = " << rand_vec_i << endl;
376 
377  OSTab tab(testOut);
378 
379  *testOut << endl << "v1 = randomize(-1,+1); ...\n" ;
380  RCP<MultiVectorBase<Scalar> > v1 = createMembers(range,loc_num_rhs);
381  rRand->randomize(v1.ptr());
382  if(dump_all()) *testOut << endl << "v1 =\n" << describe(*v1,verbLevel);
383 
384  *testOut << endl << "v2 = randomize(-1,+1); ...\n" ;
385  RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,loc_num_rhs);
386  rRand->randomize(v2.ptr());
387  if(dump_all()) *testOut << endl << "v2 =\n" << describe(*v2,verbLevel);
388 
389  *testOut << endl << "v3 = v1 + v2 ...\n" ;
390  RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,loc_num_rhs);
391  V_VpV(v3.ptr(),*v1,*v2);
392  if(dump_all()) *testOut << endl << "v3 =\n" << describe(*v3,verbLevel);
393 
394  *testOut << endl << "v4 = 0.5*op'*v3 ...\n" ;
395  RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,loc_num_rhs);
396  apply( op, CONJTRANS, *v3, v4.ptr(), d_half );
397  if(dump_all()) *testOut << endl << "v4 =\n" << describe(*v4,verbLevel);
398 
399  *testOut << endl << "v5 = op'*v1 ...\n" ;
400  RCP<MultiVectorBase<Scalar> > v5 = createMembers(domain,loc_num_rhs);
401  apply( op, CONJTRANS, *v1, v5.ptr() );
402  if(dump_all()) *testOut << endl << "v5 =\n" << describe(*v5,verbLevel);
403 
404  *testOut << endl << "v5 = 0.5*op'*v2 + 0.5*v5 ...\n" ;
405  apply( op, CONJTRANS, *v2, v5.ptr(), d_half, d_half );
406  if(dump_all()) *testOut << endl << "v5 =\n" << describe(*v5,verbLevel);
407 
408  Array<Scalar> sum_v4(loc_num_rhs), sum_v5(loc_num_rhs);
409  sums(*v4, sum_v4());
410  sums(*v5, sum_v5());
411 
412  result = testRelErrors<Scalar, Scalar, ScalarMag>(
413  "sum(v4)", sum_v4(),
414  "sum(v5)", sum_v5(),
415  "linear_properties_error_tol()", linear_properties_error_tol(),
416  "linear_properties_warning_tol()", linear_properties_warning_tol(),
417  testOut.ptr()
418  );
419  if(!result) these_results = false;
420 
421  }
422  }
423  else {
424  *testOut << endl << "Adjoint operator not supported, skipping check!\n";
425  }
426 
427  testResultsPrinter.printTestResults(these_results, inoutArg(success));
428 
429  }
430  else {
431  *out << endl << "(this->check_linear_properties()&&this->check_adjoint())==false: Skipping the check of the linear properties of the adjoint operator!\n";
432  }
433 
434  if( check_adjoint() ) {
435 
436  *out << endl << "this->check_adjoint()==true:"
437  << " Checking the agreement of the adjoint and forward operators ... ";
438 
439  TestResultsPrinter testResultsPrinter(out, show_all_tests());
440  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
441 
442  bool these_results = true;
443 
444  TEUCHOS_TEST_EQUALITY_CONST( op.opSupported(CONJTRANS), true, *testOut, these_results );
445 
446  if(result) {
447 
448  *testOut
449  << endl << "Checking that the adjoint agrees with the non-adjoint operator as:\n"
450  << endl << " <0.5*op'*v2,v1> == <v2,0.5*op*v1>"
451  << endl << " \\________/ \\_______/"
452  << endl << " v4 v3"
453  << endl << ""
454  << endl << " <v4,v1> == <v2,v3>"
455  << endl << std::flush;
456 
457  for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
458 
459  *testOut << endl << "Random vector tests = " << rand_vec_i << endl;
460 
461  OSTab tab(testOut);
462 
463  *testOut << endl << "v1 = randomize(-1,+1); ...\n" ;
464  RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,loc_num_rhs);
465  dRand->randomize(v1.ptr());
466  if(dump_all()) *testOut << endl << "v1 =\n" << describe(*v1,verbLevel);
467 
468  *testOut << endl << "v2 = randomize(-1,+1); ...\n" ;
469  RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,loc_num_rhs);
470  rRand->randomize(v2.ptr());
471  if(dump_all()) *testOut << endl << "v2 =\n" << describe(*v2,verbLevel);
472 
473  *testOut << endl << "v3 = 0.5*op*v1 ...\n" ;
474  RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,loc_num_rhs);
475  apply( op, NOTRANS, *v1, v3.ptr(), r_half );
476  if(dump_all()) *testOut << endl << "v3 =\n" << describe(*v3,verbLevel);
477 
478  *testOut << endl << "v4 = 0.5*op'*v2 ...\n" ;
479  RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,loc_num_rhs);
480  apply( op, CONJTRANS, *v2, v4.ptr(), d_half );
481  if(dump_all()) *testOut << endl << "v4 =\n" << describe(*v4,verbLevel);
482 
483  Array<Scalar> prod_v4_v1(loc_num_rhs);
484  domain->scalarProds(*v4, *v1, prod_v4_v1());
485  Array<Scalar> prod_v2_v3(loc_num_rhs);
486  range->scalarProds(*v2, *v3, prod_v2_v3());
487 
488  result = testRelErrors<Scalar, Scalar, ScalarMag>(
489  "<v4,v1>", prod_v4_v1(),
490  "<v2,v3>", prod_v2_v3(),
491  "adjoint_error_tol()", adjoint_error_tol(),
492  "adjoint_warning_tol()", adjoint_warning_tol(),
493  testOut.ptr()
494  );
495  if(!result) these_results = false;
496 
497  }
498  }
499  else {
500  *testOut << endl << "Adjoint operator not supported, skipping check!\n";
501  }
502 
503  testResultsPrinter.printTestResults(these_results, inoutArg(success));
504 
505  }
506  else {
507  *out << endl << "this->check_adjoint()==false:"
508  << " Skipping check for the agreement of the adjoint and forward operators!\n";
509  }
510 
511 
512  if( check_for_symmetry() ) {
513 
514  *out << endl << "this->check_for_symmetry()==true: Performing check of symmetry ... ";
515 
516  TestResultsPrinter testResultsPrinter(out, show_all_tests());
517  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
518 
519  bool these_results = true;
520 
521  SymmetricLinearOpTester<Scalar>::checkSymmetry(
522  op, dRand.ptr(), *testOut, loc_num_rhs,num_random_vectors(), verbLevel,dump_all(),
523  symmetry_error_tol(), symmetry_warning_tol(),
524  outArg(these_results)
525  );
526 
527  testResultsPrinter.printTestResults(these_results, inoutArg(success));
528 
529  }
530  else {
531  *out << endl << "this->check_for_symmetry()==false: Skipping check of symmetry ...\n";
532  }
533 
534  if(success)
535  *out << endl <<"Congratulations, this LinearOpBase object seems to check out!\n";
536  else
537  *out << endl <<"Oh no, at least one of the tests performed with this LinearOpBase object failed (see above failures)!\n";
538  *out << endl << "*** Leaving LinearOpTester<"<<ST::name()<<","<<ST::name()<<">::check(...)\n";
539 
540  return success;
541 
542 }
543 
544 
545 template<class Scalar>
546 bool LinearOpTester<Scalar>::check(
547  const LinearOpBase<Scalar> &op,
548  const Ptr<Teuchos::FancyOStream> &out
549  ) const
550 {
551  using Teuchos::null;
552  return check(op, null, null, out);
553 }
554 
555 
556 template<class Scalar>
557 bool LinearOpTester<Scalar>::compare(
558  const LinearOpBase<Scalar> &op1,
559  const LinearOpBase<Scalar> &op2,
560  const Ptr<MultiVectorRandomizerBase<Scalar> > &domainRandomizer,
561  const Ptr<Teuchos::FancyOStream> &out_arg
562  ) const
563 {
564 
565  using std::endl;
566  using Teuchos::rcpFromPtr;
567  using Teuchos::inoutArg;
568  using Teuchos::FancyOStream;
569  using Teuchos::OSTab;
570  typedef Teuchos::ScalarTraits<Scalar> ST;
571  bool success = true, result;
572  const int loc_num_rhs = this->num_rhs();
573  const Scalar r_half = Scalar(0.5)*ST::one();
574  const RCP<FancyOStream> out = rcpFromPtr(out_arg);
575  const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
576 
577  OSTab tab(out,1,"THYRA");
578 
579  if(out.get()) {
580  *out
581  << endl << "*** Entering LinearOpTester<"<<ST::name()<<","<<ST::name()<<">::compare(op1,op2,...) ...\n";
582  if(show_all_tests())
583  *out << endl << "describe op1:\n" << Teuchos::describe(op1,verbLevel);
584  else
585  *out << endl << "describe op1: " << op1.description() << endl;
586  if(show_all_tests())
587  *out << endl << "describe op2:\n" << Teuchos::describe(op2,verbLevel);
588  else
589  *out << endl << "describe op2: " << op2.description() << endl;
590  }
591 
592  RCP<MultiVectorRandomizerBase<Scalar> > dRand;
593  if (nonnull(domainRandomizer)) dRand = rcpFromPtr(domainRandomizer);
594  else dRand = universalMultiVectorRandomizer<Scalar>();
595 
596  RCP<const VectorSpaceBase<Scalar> > range = op1.range();
597  RCP<const VectorSpaceBase<Scalar> > domain = op1.domain();
598 
599  if(out.get()) *out << endl << "Checking that range and domain spaces are compatible ... ";
600 
601  {
602 
603  TestResultsPrinter testResultsPrinter(out, show_all_tests());
604  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
605 
606  bool these_results = true;
607 
608  *testOut << endl << "op1.domain()->isCompatible(*op2.domain()) ? ";
609  result = op1.domain()->isCompatible(*op2.domain());
610  if(!result) these_results = false;
611  *testOut << passfail(result) << endl;
612 
613  *testOut << endl << "op1.range()->isCompatible(*op2.range()) ? ";
614  result = op1.range()->isCompatible(*op2.range());
615  if(!result) these_results = false;
616  *testOut << passfail(result) << endl;
617 
618  testResultsPrinter.printTestResults(these_results, inoutArg(success));
619 
620  }
621 
622  if(!success) {
623  if(out.get()) *out << endl << "Skipping further checks since operators are not compatible!\n";
624  return success;
625  }
626 
627  if(out.get()) *out << endl << "Checking that op1 == op2 ... ";
628 
629  {
630 
631  TestResultsPrinter testResultsPrinter(out, show_all_tests());
632  const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();
633 
634  bool these_results = true;
635 
636  *testOut
637  << endl << "Checking that op1 and op2 produce the same results:\n"
638  << endl << " 0.5*op1*v1 == 0.5*op2*v1"
639  << endl << " \\________/ \\________/"
640  << endl << " v2 v3"
641  << endl << ""
642  << endl << " norm(v2-v3) ~= 0"
643  // << endl << " |sum(v2)| == |sum(v3)|"
644  << endl << std::flush;
645 
646  for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
647 
648  *testOut << endl << "Random vector tests = " << rand_vec_i << endl;
649 
650  OSTab tab2(testOut);
651 
652  if(dump_all()) *testOut << endl << "v1 = randomize(-1,+1); ...\n" ;
653  RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,loc_num_rhs);
654  dRand->randomize(v1.ptr());
655  if(dump_all()) *testOut << endl << "v1 =\n" << *v1;
656 
657  if(dump_all()) *testOut << endl << "v2 = 0.5*op1*v1 ...\n" ;
658  RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,loc_num_rhs);
659  apply( op1, NOTRANS, *v1, v2.ptr(), r_half );
660  if(dump_all()) *testOut << endl << "v2 =\n" << *v2;
661 
662  if(dump_all()) *testOut << endl << "v3 = 0.5*op2*v1 ...\n" ;
663  RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,loc_num_rhs);
664  apply( op2, NOTRANS, *v1, v3.ptr(), r_half );
665  if(dump_all()) *testOut << endl << "v3 =\n" << *v3;
666 
667  // check error in each column
668  for(int col_id=0;col_id < v1->domain()->dim();col_id++) {
669  std::stringstream ss;
670  ss << ".col[" << col_id << "]";
671 
672  result = Thyra::testRelNormDiffErr(
673  "v2"+ss.str(),*v2->col(col_id),
674  "v3"+ss.str(),*v3->col(col_id),
675  "linear_properties_error_tol()", linear_properties_error_tol(),
676  "linear_properties_warning_tol()", linear_properties_warning_tol(),
677  &*testOut);
678  if(!result) these_results = false;
679  }
680  /*
681  Array<Scalar> sum_v2(loc_num_rhs), sum_v3(loc_num_rhs);
682  sums(*v2,&sum_v2[0]);
683  sums(*v3,&sum_v3[0]);
684 
685  result = testRelErrors(
686  loc_num_rhs
687  ,"sum(v2)", &sum_v2[0]
688  ,"sum(v3)", &sum_v3[0]
689  ,"linear_properties_error_tol()", linear_properties_error_tol()
690  ,"linear_properties_warning_tol()", linear_properties_warning_tol()
691  ,inOutArg(testOut)
692  );
693  */
694  if(!result) these_results = false;
695 
696  }
697 
698  testResultsPrinter.printTestResults(these_results, inoutArg(success));
699 
700  }
701 
702  if(out.get()) {
703  if(success)
704  *out << endl <<"Congratulations, these two LinearOpBase objects seem to be the same!\n";
705  else
706  *out << endl <<"Oh no, these two LinearOpBase objects seem to be different (see above failures)!\n";
707  *out << endl << "*** Leaving LinearOpTester<"<<ST::name()<<","<<ST::name()<<">::compare(...)\n";
708  }
709 
710  return success;
711 
712 }
713 
714 
715 template<class Scalar>
716 bool LinearOpTester<Scalar>::compare(
717  const LinearOpBase<Scalar> &op1,
718  const LinearOpBase<Scalar> &op2,
719  const Ptr<Teuchos::FancyOStream> &out_arg
720  ) const
721 {
722  return compare(op1, op2, Teuchos::null, out_arg);
723 }
724 
725 
726 // private
727 
728 
729 // Nonmember helper
730 template<class ScalarMag>
731 inline
732 void setDefaultTol(const ScalarMag &defaultDefaultTol,
733  ScalarMag &defaultTol)
734 {
735  typedef ScalarTraits<ScalarMag> SMT;
736  if (defaultTol < SMT::zero()) {
737  defaultTol = defaultDefaultTol;
738  }
739 }
740 
741 
742 template<class Scalar>
743 void LinearOpTester<Scalar>::setDefaultTols()
744 {
745  typedef ScalarTraits<ScalarMag> SMT;
746  const ScalarMag defaultWarningTol = 1e+2 * SMT::eps();
747  const ScalarMag defaultErrorTol = 1e+4 * SMT::eps();
748  setDefaultTol(defaultWarningTol, linear_properties_warning_tol_);
749  setDefaultTol(defaultErrorTol, linear_properties_error_tol_);
750  setDefaultTol(defaultWarningTol, adjoint_warning_tol_);
751  setDefaultTol(defaultErrorTol, adjoint_error_tol_);
752  setDefaultTol(defaultWarningTol, symmetry_warning_tol_);
753  setDefaultTol(defaultErrorTol, symmetry_error_tol_);
754 }
755 
756 
757 } // namespace Thyra
758 
759 
760 #endif // THYRA_LINEAR_OP_TESTER_DEF_HPP
Thyra::LinearOpTester::enable_all_tests
void enable_all_tests(const bool enable_all_tests)
Enable or disable all tests.
Definition: Thyra_LinearOpTester_def.hpp:142
is_null
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::OSTab
basic_OSTab< char > OSTab
Teuchos::VERB_LOW
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::basic_oblackholestream
Teuchos::VERB_EXTREME
Thyra::testRelNormDiffErr
bool testRelNormDiffErr(const std::string &v1_name, const VectorBase< Scalar > &v1, const std::string &v2_name, const VectorBase< Scalar > &v2, const std::string &maxRelErr_error_name, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType &maxRelErr_error, const std::string &maxRelErr_warning_name, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType &maxRelErr_warning, std::ostream *out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_LOW, const std::string &leadingIndent=std::string(""))
Compute, check and optionally print the relative errors in two vectors.
Definition: Thyra_TestingTools.hpp:106
Thyra::CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types)...
Definition: Thyra_OperatorVectorTypes.hpp:140
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::VERB_MEDIUM
Thyra::LinearOpTester::compare
bool compare(const LinearOpBase< Scalar > &op1, const LinearOpBase< Scalar > &op2, const Ptr< MultiVectorRandomizerBase< Scalar > > &domainRandomizer, const Ptr< FancyOStream > &out_arg) const
Check if two linear operators are the same or not.
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::basic_FancyOStream
Thyra::LinearOpBase
Base class for all linear operators.
Definition: Thyra_LinearOpBase_decl.hpp:159
Thyra::LinearOpTester::set_all_warning_tol
void set_all_warning_tol(const ScalarMag warning_tol)
Set all the warning tolerances to the same value.
Definition: Thyra_LinearOpTester_def.hpp:151
nonnull
bool nonnull(const boost::shared_ptr< T > &p)
Thyra::LinearOpTester::check
bool check(const LinearOpBase< Scalar > &op, const Ptr< MultiVectorRandomizerBase< Scalar > > &rangeRandomizer, const Ptr< MultiVectorRandomizerBase< Scalar > > &domainRandomizer, const Ptr< FancyOStream > &out) const
Check a linear operator.
Teuchos::as
TypeTo as(const TypeFrom &t)
Thyra::LinearOpTester::ScalarMag
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
Local typedef for promoted scalar magnitude.
Definition: Thyra_LinearOpTester_decl.hpp:46
Thyra::LinearOpTester::LinearOpTester
LinearOpTester()
Default constructor which sets default parameter values.
Definition: Thyra_LinearOpTester_def.hpp:122
Thyra::MultiVectorRandomizerBase
Base interface for a strategy object for randomizing a multi-vector.
Definition: Thyra_MultiVectorRandomizerBase.hpp:32
TEUCHOS_TEST_EQUALITY_CONST
#define TEUCHOS_TEST_EQUALITY_CONST(v1, v2, out, success)
Thyra::LinearOpTester::set_all_error_tol
void set_all_error_tol(const ScalarMag error_tol)
Set all the error tolerances to the same value.
Definition: Thyra_LinearOpTester_def.hpp:160
Generated on Thu Oct 24 2024 09:21:23 for Thyra by   doxygen 1.8.5