Thyra Version of the Day
Loading...
Searching...
No Matches
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
55namespace Thyra {
56
57
58template<class Scalar>
59class SymmetricLinearOpTester {
60public:
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
153template<class Scalar>
154LinearOpTester<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
173template<class Scalar>
174void 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
182template<class Scalar>
183void 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
191template<class Scalar>
192void 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
200template<class Scalar>
201bool 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
577template<class Scalar>
578bool 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
588template<class Scalar>
589bool 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
747template<class Scalar>
748bool 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
762template<class ScalarMag>
763inline
764void 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
774template<class Scalar>
775void 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
Teuchos_TestingHelpers.hpp
Teuchos::Describable::description
virtual std::string description() const
Teuchos::RCP::ptr
Ptr< T > ptr() const
Teuchos::RCP::get
T * get() const
Thyra::LinearOpBase
Base class for all linear operators.
Definition Thyra_LinearOpBase_decl.hpp:191
Thyra::LinearOpBase::apply
void apply(const LinearOpBase< Scalar > &M, const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha=static_cast< Scalar >(1.0), const Scalar beta=static_cast< Scalar >(0.0))
Non-member function call for M.apply(...).
Definition Thyra_LinearOpBase_def.hpp:73
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:242
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.
Definition Thyra_LinearOpTester_def.hpp:589
Thyra::LinearOpTester::LinearOpTester
LinearOpTester()
Default constructor which sets default parameter values.
Definition Thyra_LinearOpTester_def.hpp:154
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
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
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
Thyra::LinearOpTester::ScalarMag
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
Local typedef for promoted scalar magnitude.
Definition Thyra_LinearOpTester_decl.hpp:79
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.
Definition Thyra_LinearOpTester_def.hpp:201
Thyra::MultiVectorBase::sums
void sums(const MultiVectorBase< Scalar > &V, const ArrayView< Scalar > &sums)
Multi-vector column sum.
Definition Thyra_MultiVectorStdOps_def.hpp:81
Thyra::MultiVectorBase::V_VpV
void V_VpV(const Ptr< MultiVectorBase< Scalar > > &Z, const MultiVectorBase< Scalar > &X, const MultiVectorBase< Scalar > &Y)
Z(i,j) = X(i,j) + Y(i,j), i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
Definition Thyra_MultiVectorStdOps_def.hpp:273
Thyra::MultiVectorRandomizerBase
Base interface for a strategy object for randomizing a multi-vector.
Definition Thyra_MultiVectorRandomizerBase.hpp:64
Thyra::TestResultsPrinter
Control printing of test results.
Definition Thyra_TestingToolsDecl.hpp:239
Thyra::TestResultsPrinter::getTestOStream
RCP< FancyOStream > getTestOStream()
Return the stream used for testing.
Definition Thyra_TestingTools.cpp:134
Thyra::TestResultsPrinter::printTestResults
void printTestResults(const bool this_result, const Ptr< bool > &success)
Print the test result.
Definition Thyra_TestingTools.cpp:140
Thyra::UniversalMultiVectorRandomizer::universalMultiVectorRandomizer
RCP< UniversalMultiVectorRandomizer< Scalar > > universalMultiVectorRandomizer()
Nonmember constructor.
Definition Thyra_UniversalMultiVectorRandomizer.hpp:92
Thyra::VectorSpaceBase::createMembers
RCP< MultiVectorBase< Scalar > > createMembers(const RCP< const VectorSpaceBase< Scalar > > &vs, int numMembers, const std::string &label="")
Create a set of vector members (a MultiVectorBase) from the vector space.
Thyra::testRelErrors
bool testRelErrors(const std::string &v1_name, const ArrayView< const Scalar1 > &v1, const std::string &v2_name, const ArrayView< const Scalar2 > &v2, const std::string &maxRelErr_error_name, const ScalarMag &maxRelErr_error, const std::string &maxRelErr_warning_name, const ScalarMag &maxRelErr_warning, const Ptr< std::ostream > &out, const std::string &leadingIndent=std::string(""))
Compute, check and optionally print the relative errors in two scalar arays.
Definition Thyra_TestingTools.hpp:82
Thyra::passfail
const std::string passfail(const bool result)
Definition Thyra_TestingToolsDecl.hpp:56
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::NOTRANS
@ NOTRANS
Use the non-transposed operator.
Definition Thyra_OperatorVectorTypes.hpp:162
Thyra::CONJTRANS
@ CONJTRANS
Use the transposed operator with complex-conjugate clements (same as TRANS for real scalar types).
Definition Thyra_OperatorVectorTypes.hpp:172
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::FancyOStream
basic_FancyOStream< char > FancyOStream
Teuchos::OSTab
basic_OSTab< char > OSTab
Teuchos::oblackholestream
basic_oblackholestream< char, std::char_traits< char > > oblackholestream
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.13.2