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