Thyra Version of the Day
Loading...
Searching...
No Matches
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
23namespace Thyra {
24
25
26template<class Scalar>
27class SymmetricLinearOpTester {
28public:
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
121template<class Scalar>
122LinearOpTester<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
141template<class Scalar>
142void 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
150template<class Scalar>
151void 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
159template<class Scalar>
160void 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
168template<class Scalar>
169bool 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
545template<class Scalar>
546bool 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
556template<class Scalar>
557bool 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
715template<class Scalar>
716bool 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
730template<class ScalarMag>
731inline
732void 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
742template<class Scalar>
743void 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
Teuchos::basic_FancyOStream
Teuchos::basic_OSTab
Teuchos::basic_oblackholestream
Thyra::LinearOpBase
Base class for all linear operators.
Definition Thyra_LinearOpBase_decl.hpp:159
Thyra::LinearOpBase::range
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
Thyra::LinearOpBase::domain
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Thyra::LinearOpBase::opSupported
bool opSupported(EOpTransp M_trans) const
Return if the M_trans operation of apply() is supported or not.
Definition Thyra_LinearOpBase_decl.hpp:210
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.
Thyra::LinearOpTester::LinearOpTester
LinearOpTester()
Default constructor which sets default parameter values.
Definition Thyra_LinearOpTester_def.hpp:122
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
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
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
Thyra::LinearOpTester::ScalarMag
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
Local typedef for promoted scalar magnitude.
Definition Thyra_LinearOpTester_decl.hpp:46
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.
Thyra::MultiVectorRandomizerBase
Base interface for a strategy object for randomizing a multi-vector.
Definition Thyra_MultiVectorRandomizerBase.hpp:32
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Teuchos::nonnull
bool nonnull(const std::shared_ptr< T > &p)
Thyra::passfail
const std::string passfail(const bool result)
Definition Thyra_TestingToolsDecl.hpp:24
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::NOTRANS
@ NOTRANS
Use the non-transposed operator.
Definition Thyra_OperatorVectorTypes.hpp:130
Thyra::CONJTRANS
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
Definition Thyra_OperatorVectorTypes.hpp:140
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::dyn_cast
T_To & dyn_cast(T_From &from)
Teuchos::OSTab
basic_OSTab< char > OSTab
Teuchos::EVerbosityLevel
EVerbosityLevel
Teuchos::VERB_MEDIUM
VERB_MEDIUM
Teuchos::VERB_EXTREME
VERB_EXTREME
Teuchos::VERB_LOW
VERB_LOW
TEUCHOS_TEST_EQUALITY_CONST
#define TEUCHOS_TEST_EQUALITY_CONST(v1, v2, out, success)
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Generated by doxygen 1.9.8