Thyra_LinearOpTester.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_LINEAR_OP_TESTER_HPP
00030 #define THYRA_LINEAR_OP_TESTER_HPP
00031 
00032 #include "Thyra_LinearOpTesterDecl.hpp"
00033 #include "Thyra_LinearOpBase.hpp"
00034 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00035 #include "Thyra_describeLinearOp.hpp"
00036 #include "Thyra_VectorStdOps.hpp"
00037 #include "Thyra_TestingTools.hpp"
00038 #include "Thyra_UniversalMultiVectorRandomizer.hpp"
00039 
00040 namespace Thyra {
00041 
00042 // SymmetricLinearOpTester (using partial specialization only test symmetry on
00043 // operators where RangeScalar and DomainScalar are the same)
00044 
00045 template<class RangeScalar, class DomainScalar>
00046 class SymmetricLinearOpTester {
00047 public:
00048   typedef typename Teuchos::PromotionTraits<RangeScalar,DomainScalar>::promote Scalar;
00049   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00050   static void checkSymmetry(
00051     const LinearOpBase<RangeScalar,DomainScalar>  &op
00052     ,MultiVectorRandomizerBase<DomainScalar>      *dRand
00053     ,Teuchos::FancyOStream                        &oss
00054     ,const int                                    num_rhs
00055     ,const int                                    num_random_vectors
00056     ,const Teuchos::EVerbosityLevel               verbLevel
00057     ,const bool                                   dump_all
00058     ,const ScalarMag                              &symmetry_error_tol
00059     ,const ScalarMag                              &symmetry_warning_tol
00060     ,bool                                         *these_results
00061     )
00062     {
00063       using std::endl;
00064       typedef Teuchos::ScalarTraits<RangeScalar>  RST;
00065       typedef Teuchos::ScalarTraits<DomainScalar> DST;
00066       oss << endl << "RangeScalar = "<<RST::name()<<" == DomainScalar = "<<DST::name()<<": failed, the opeator can not be symmetric!\n";
00067       *these_results = false;
00068     }
00069 };
00070 
00071 template<class Scalar>
00072 class SymmetricLinearOpTester<Scalar,Scalar> {
00073 public:
00074   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00075   static void checkSymmetry(
00076     const LinearOpBase<Scalar>                    &op
00077     ,MultiVectorRandomizerBase<Scalar>            *dRand
00078     ,Teuchos::FancyOStream                        &oss
00079     ,const int                                    num_rhs
00080     ,const int                                    num_random_vectors
00081     ,const Teuchos::EVerbosityLevel               verbLevel
00082     ,const bool                                   dump_all
00083     ,const ScalarMag                              &symmetry_error_tol
00084     ,const ScalarMag                              &symmetry_warning_tol
00085     ,bool                                         *these_results
00086     )
00087     {
00088 
00089       using std::endl;
00090 
00091       bool result;
00092       using Teuchos::OSTab;
00093       typedef Teuchos::ScalarTraits<Scalar> ST;
00094       const Scalar half = Scalar(0.4)*ST::one();
00095       Teuchos::RCP<const VectorSpaceBase<Scalar> > domain = op.domain();
00096       
00097       oss << endl << "op.domain()->isCompatible(*op.range()) == true : ";
00098       result = op.domain()->isCompatible(*op.range());
00099       if(!result) *these_results = false;
00100       oss << passfail(result) << endl;
00101       
00102       if(result) {
00103         
00104         oss
00105           << endl << "Checking that the operator is symmetric as:\n"
00106           << endl << "  <0.5*op*v2,v1> == <v2,0.5*op*v1>"
00107           << endl << "   \\_______/            \\_______/"
00108           << endl << "      v4                    v3"
00109           << endl << ""
00110           << endl << "         <v4,v1> == <v2,v3>"
00111           << endl;
00112         
00113         for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors; ++rand_vec_i ) {
00114           
00115           oss << endl << "Random vector tests = " << rand_vec_i << endl;
00116 
00117           OSTab tab(Teuchos::rcp(&oss,false));
00118 
00119           if(dump_all) oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00120           Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,num_rhs);
00121           dRand->randomize(&*v1);
00122           if(dump_all) oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00123           
00124           if(dump_all) oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00125           Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,num_rhs);
00126           dRand->randomize(&*v2);
00127           if(dump_all) oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00128           
00129           if(dump_all) oss << endl << "v3 = 0.5*op*v1 ...\n" ;
00130           Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,num_rhs);
00131           apply( op, NONCONJ_ELE, *v1, &*v3, half );
00132          if(dump_all) oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00133           
00134           if(dump_all) oss << endl << "v4 = 0.5*op*v2 ...\n" ;
00135           Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,num_rhs);
00136           apply( op, NONCONJ_ELE, *v2, &*v4, half );
00137           if(dump_all) oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00138 
00139           std::vector<Scalar> prod1(num_rhs), prod2(num_rhs);
00140           domain->scalarProds(*v4,*v1,&prod1[0]);
00141           domain->scalarProds(*v2,*v3,&prod2[0]);
00142           
00143           result = testRelErrors(
00144             num_rhs
00145             ,"<v4,v1>", &prod1[0]
00146             ,"<v2,v3>", &prod2[0]
00147             ,"symmetry_error_tol()", symmetry_error_tol
00148             ,"symmetry_warning_tol()", symmetry_warning_tol
00149             ,&oss
00150             );
00151           if(!result) *these_results = false;
00152         
00153         }
00154       }
00155       else {
00156         oss << endl << "Range and domain spaces are different, skipping check!\n";
00157       }
00158     }
00159 };
00160 
00161 // LinearOpTester
00162 
00163 template<class RangeScalar, class DomainScalar>
00164 LinearOpTester<RangeScalar,DomainScalar>::LinearOpTester(
00165   const bool          check_linear_properties
00166   ,const ScalarMag    linear_properties_warning_tol
00167   ,const ScalarMag    linear_properties_error_tol
00168   ,const bool         check_adjoint
00169   ,const ScalarMag    adjoint_warning_tol
00170   ,const ScalarMag    adjoint_error_tol
00171   ,const bool         check_for_symmetry
00172   ,const ScalarMag    symmetry_warning_tol
00173   ,const ScalarMag    symmetry_error_tol
00174   ,const int          num_random_vectors
00175   ,const bool         show_all_tests
00176   ,const bool         dump_all
00177   ,const int          num_rhs
00178   )
00179   :check_linear_properties_(check_linear_properties)
00180   ,linear_properties_warning_tol_(linear_properties_warning_tol)
00181   ,linear_properties_error_tol_(linear_properties_error_tol)
00182   ,check_adjoint_(check_adjoint)
00183   ,adjoint_warning_tol_(adjoint_warning_tol)
00184   ,adjoint_error_tol_(adjoint_error_tol)
00185   ,check_for_symmetry_(check_for_symmetry)
00186   ,symmetry_warning_tol_(symmetry_warning_tol)
00187   ,symmetry_error_tol_(symmetry_error_tol)
00188   ,num_random_vectors_(num_random_vectors)
00189   ,show_all_tests_(show_all_tests)
00190   ,dump_all_(dump_all)
00191   ,num_rhs_(num_rhs)
00192 {}
00193 
00194 template<class RangeScalar, class DomainScalar>
00195 void LinearOpTester<RangeScalar,DomainScalar>::enable_all_tests( const bool enable_all_tests )
00196 {
00197   check_linear_properties_ = enable_all_tests;
00198   check_adjoint_           = enable_all_tests;
00199   check_for_symmetry_      = enable_all_tests;
00200 }
00201 
00202 template<class RangeScalar, class DomainScalar>
00203 void LinearOpTester<RangeScalar,DomainScalar>::set_all_warning_tol( const ScalarMag warning_tol )
00204 {
00205   linear_properties_warning_tol_  = warning_tol;
00206   adjoint_warning_tol_            = warning_tol;
00207   symmetry_warning_tol_           = warning_tol;
00208 }
00209 
00210 template<class RangeScalar, class DomainScalar>
00211 void LinearOpTester<RangeScalar,DomainScalar>::set_all_error_tol( const ScalarMag error_tol )
00212 {
00213   linear_properties_error_tol_  = error_tol;
00214   adjoint_error_tol_            = error_tol;
00215   symmetry_error_tol_           = error_tol;
00216 }
00217 
00218 template<class RangeScalar, class DomainScalar>
00219 bool LinearOpTester<RangeScalar,DomainScalar>::check(
00220   const LinearOpBase<RangeScalar,DomainScalar>  &op
00221   ,MultiVectorRandomizerBase<RangeScalar>       *rangeRandomizer
00222   ,MultiVectorRandomizerBase<DomainScalar>      *domainRandomizer
00223   ,Teuchos::FancyOStream                        *out_arg
00224   ) const
00225 {
00226 
00227   using std::endl;
00228   using Teuchos::FancyOStream;
00229   using Teuchos::OSTab;
00230   typedef Teuchos::ScalarTraits<RangeScalar>  RST;
00231   typedef Teuchos::ScalarTraits<DomainScalar> DST;
00232   bool success = true, result;
00233   const int          num_rhs = this->num_rhs();
00234   const RangeScalar  r_one  = RST::one();
00235   const DomainScalar d_one  = DST::one();
00236   const RangeScalar  r_half = RangeScalar(0.5)*r_one;
00237   const DomainScalar d_half = DomainScalar(0.5)*d_one;
00238   Teuchos::RCP<FancyOStream> out = Teuchos::rcp(out_arg,false);
00239   const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00240 
00241   OSTab tab(out,1,"THYRA");
00242 
00243   // ToDo 04/28/2005:
00244   // * Test the MultiVectorBase apply() function and output to the VectorBase apply() function!
00245 
00246   if(out.get()) {
00247     *out << endl << "*** Entering LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::check(op,...) ...\n";
00248     if(show_all_tests()) {
00249       *out << endl << "describe op:\n" << Teuchos::describe(op,verbLevel);
00250 /*
00251       if(op.applyTransposeSupports(CONJ_ELE) && verbLevel==Teuchos::VERB_EXTREME) {
00252         *out << endl << "describe adjoint op:\n";
00253         describeLinearOp(*adjoint(Teuchos::rcp(&op,false)),*out,verbLevel);
00254       }
00255 */
00256     }
00257     else {
00258       *out << endl << "describe op:\n" << Teuchos::describe(op,Teuchos::VERB_LOW);
00259     }
00260   }
00261 
00262   Teuchos::RCP< MultiVectorRandomizerBase<RangeScalar> >  rRand;
00263   if(rangeRandomizer)   rRand = Teuchos::rcp(rangeRandomizer,false);
00264   else                  rRand = Teuchos::rcp(new UniversalMultiVectorRandomizer<RangeScalar>());
00265   Teuchos::RCP< MultiVectorRandomizerBase<DomainScalar> > dRand;
00266   if(domainRandomizer)  dRand = Teuchos::rcp(domainRandomizer,false);
00267   else                  dRand = Teuchos::rcp(new UniversalMultiVectorRandomizer<DomainScalar>());
00268   
00269   if(out.get())
00270     *out << endl << "Checking the domain and range spaces ... ";
00271 
00272   Teuchos::RCP<const VectorSpaceBase<RangeScalar> >  range  = op.range();
00273   Teuchos::RCP<const VectorSpaceBase<DomainScalar> > domain = op.domain();
00274   
00275   {
00276 
00277     std::ostringstream ossStore;
00278     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00279     if(out.get()) ossStore.copyfmt(*out);
00280     bool these_results = true;
00281 
00282     *oss << endl << "op.domain().get() != NULL ? ";
00283     result = domain.get() != NULL;
00284     if(!result) these_results = false;
00285     *oss << passfail(result) << endl;
00286     
00287     *oss << endl << "op.range().get() != NULL ? ";
00288     result = range.get() != NULL;
00289     if(!result) these_results = false;
00290     *oss << passfail(result) << endl;
00291 
00292     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00293 
00294   }
00295 
00296   if( check_linear_properties() ) {
00297 
00298     if(out.get()) *out << endl << "this->check_linear_properties()==true: Checking the linear properties of the forward linear operator ... ";
00299 
00300     std::ostringstream ossStore;
00301     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00302     if(out.get()) ossStore.copyfmt(*out);
00303     bool these_results = true;
00304 
00305     *oss << endl << "op.applySupports(NONCONJ_ELE) == true ? ";
00306     result = op.applySupports(NONCONJ_ELE);
00307     if(!result) these_results = false;
00308     *oss << passfail(result) << endl;
00309 
00310     if(result) {
00311     
00312       *oss
00313         << endl << "Checking that the forward operator is truly linear:\n"
00314         << endl << "  0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2"
00315         << endl << "          \\_____/         \\___/"
00316         << endl << "             v3            v5"
00317         << endl << "  \\_____________/     \\___________________/"
00318         << endl << "         v4                    v5"
00319         << endl << ""
00320         << endl << "           sum(v4) == sum(v5)"
00321         << endl;
00322       
00323       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00324         
00325         *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00326 
00327         OSTab tab(oss);
00328         
00329         *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00330         Teuchos::RCP<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00331         dRand->randomize(&*v1);
00332         if(dump_all()) *oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00333         
00334         *oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00335         Teuchos::RCP<MultiVectorBase<DomainScalar> > v2 = createMembers(domain,num_rhs);
00336         dRand->randomize(&*v2);
00337         if(dump_all()) *oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00338         
00339         *oss << endl << "v3 = v1 + v2 ...\n" ;
00340         Teuchos::RCP<MultiVectorBase<DomainScalar> > v3 = createMembers(domain,num_rhs);
00341         V_VpV(&*v3,*v1,*v2);
00342         if(dump_all()) *oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00343         
00344         *oss << endl << "v4 = 0.5*op*v3 ...\n" ;
00345         Teuchos::RCP<MultiVectorBase<RangeScalar> > v4 = createMembers(range,num_rhs);
00346         apply( op, NONCONJ_ELE, *v3, &*v4, r_half );
00347         if(dump_all()) *oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00348         
00349         *oss << endl << "v5 = op*v1 ...\n" ;
00350         Teuchos::RCP<MultiVectorBase<RangeScalar> > v5 = createMembers(range,num_rhs);
00351         apply( op, NONCONJ_ELE, *v1, &*v5 );
00352         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00353         
00354         *oss << endl << "v5 = 0.5*op*v2 + 0.5*v5 ...\n" ;
00355         apply( op, NONCONJ_ELE, *v2, &*v5, r_half, r_half );
00356         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00357 
00358         std::vector<RangeScalar> sum_v4(num_rhs), sum_v5(num_rhs);
00359         sums(*v4,&sum_v4[0]);
00360         sums(*v5,&sum_v5[0]);
00361         
00362         result = testRelErrors(
00363           num_rhs
00364           ,"sum(v4)", &sum_v4[0]
00365           ,"sum(v5)", &sum_v5[0]
00366           ,"linear_properties_error_tol()", linear_properties_error_tol()
00367           ,"linear_properties_warning_tol()", linear_properties_warning_tol()
00368           ,&*oss
00369           );
00370         if(!result) these_results = false;
00371         
00372       }
00373     }
00374     else {
00375       *oss << endl << "Forward operator not supported, skipping check!\n";
00376     }
00377 
00378     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00379 
00380   }
00381   else {
00382     if(out.get()) *out << endl << "this->check_linear_properties()==false: Skipping the check of the linear properties of the forward operator!\n";
00383   }
00384 
00385   if( check_linear_properties() && check_adjoint() ) {
00386 
00387     if(out.get()) *out << endl << "(this->check_linear_properties()&&this->check_adjoint())==true: Checking the linear properties of the adjoint operator ... ";
00388 
00389     std::ostringstream ossStore;
00390     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00391     if(out.get()) ossStore.copyfmt(*out);
00392     bool these_results = true;
00393 
00394     *oss << endl << "op.applyTransposeSupports(CONJ_ELE) == true ? ";
00395     result = op.applyTransposeSupports(CONJ_ELE);
00396     if(!result) these_results = false;
00397     *oss << passfail(result) << endl;
00398 
00399     if(result) {
00400     
00401       *oss
00402         << endl << "Checking that the adjoint operator is truly linear:\n"
00403         << endl << "  0.5*op'*(v1 + v2) == 0.5*op'*v1 + 0.5*op'*v2"
00404         << endl << "           \\_____/         \\____/"
00405         << endl << "              v3             v5"
00406         << endl << "  \\_______________/    \\_____________________/"
00407         << endl << "         v4                      v5"
00408         << endl << ""
00409         << endl << "           sum(v4) == sum(v5)"
00410         << endl;
00411       
00412       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00413         
00414         *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00415 
00416         OSTab tab(oss);
00417         
00418         *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00419         Teuchos::RCP<MultiVectorBase<RangeScalar> > v1 = createMembers(range,num_rhs);
00420         rRand->randomize(&*v1);
00421         if(dump_all()) *oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00422         
00423         *oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00424         Teuchos::RCP<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00425         rRand->randomize(&*v2);
00426         if(dump_all()) *oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00427         
00428         *oss << endl << "v3 = v1 + v2 ...\n" ;
00429         Teuchos::RCP<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00430         V_VpV(&*v3,*v1,*v2);
00431         if(dump_all()) *oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00432         
00433         *oss << endl << "v4 = 0.5*op'*v3 ...\n" ;
00434         Teuchos::RCP<MultiVectorBase<DomainScalar> > v4 = createMembers(domain,num_rhs);
00435         applyTranspose( op, CONJ_ELE, *v3, &*v4, d_half );
00436         if(dump_all()) *oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00437         
00438         *oss << endl << "v5 = op'*v1 ...\n" ;
00439         Teuchos::RCP<MultiVectorBase<DomainScalar> > v5 = createMembers(domain,num_rhs);
00440         applyTranspose( op, CONJ_ELE, *v1, &*v5 );
00441         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00442         
00443         *oss << endl << "v5 = 0.5*op'*v2 + 0.5*v5 ...\n" ;
00444         applyTranspose( op, CONJ_ELE, *v2, &*v5, d_half, d_half );
00445         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00446         
00447 
00448         std::vector<DomainScalar> sum_v4(num_rhs), sum_v5(num_rhs);
00449         sums(*v4,&sum_v4[0]);
00450         sums(*v5,&sum_v5[0]);
00451         
00452         result = testRelErrors(
00453           num_rhs
00454           ,"sum(v4)", &sum_v4[0]
00455           ,"sum(v5)", &sum_v5[0]
00456           ,"linear_properties_error_tol()", linear_properties_error_tol()
00457           ,"linear_properties_warning_tol()", linear_properties_warning_tol()
00458           ,&*oss
00459           );
00460         if(!result) these_results = false;
00461         
00462       }
00463     }
00464     else {
00465       *oss << endl << "Adjoint operator not supported, skipping check!\n";
00466     }
00467 
00468     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00469 
00470   }
00471   else {
00472     if(out.get()) *out << endl << "(this->check_linear_properties()&&this->check_adjoint())==false: Skipping the check of the linear properties of the adjoint operator!\n";
00473   }
00474   
00475   if( check_adjoint() ) {
00476 
00477     if(out.get()) *out << endl << "this->check_adjoint()==true: Checking the agreement of the adjoint and forward operators ... ";
00478 
00479     std::ostringstream ossStore;
00480     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00481     if(out.get()) ossStore.copyfmt(*out);
00482     bool these_results = true;
00483     
00484     *oss << endl << "op.applyTransposeSupports(CONJ_ELE) == true ? ";
00485     result = op.applyTransposeSupports(CONJ_ELE);
00486     if(!result) these_results = false;
00487     *oss << passfail(result) << endl;
00488 
00489     if(result) {
00490     
00491       *oss
00492         << endl << "Checking that the adjoint agrees with the non-adjoint operator as:\n"
00493         << endl << "  <0.5*op'*v2,v1> == <v2,0.5*op*v1>"
00494         << endl << "   \\________/            \\_______/"
00495         << endl << "       v4                   v3"
00496         << endl << ""
00497         << endl << "         <v4,v1>  == <v2,v3>"
00498         << endl;
00499     
00500       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00501       
00502         *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00503 
00504         OSTab tab(oss);
00505       
00506         *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00507         Teuchos::RCP<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00508         dRand->randomize(&*v1);
00509         if(dump_all()) *oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00510       
00511         *oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00512         Teuchos::RCP<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00513         rRand->randomize(&*v2);
00514         if(dump_all()) *oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00515       
00516         *oss << endl << "v3 = 0.5*op*v1 ...\n" ;
00517         Teuchos::RCP<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00518         apply( op, NONCONJ_ELE, *v1, &*v3, r_half );
00519         if(dump_all()) *oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00520       
00521         *oss << endl << "v4 = 0.5*op'*v2 ...\n" ;
00522         Teuchos::RCP<MultiVectorBase<DomainScalar> > v4 = createMembers(domain,num_rhs);
00523         applyTranspose( op, CONJ_ELE, *v2, &*v4, d_half );
00524         if(dump_all()) *oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00525 
00526         std::vector<DomainScalar> prod_v4_v1(num_rhs);
00527         domain->scalarProds(*v4,*v1,&prod_v4_v1[0]);
00528         std::vector<RangeScalar> prod_v2_v3(num_rhs);
00529         range->scalarProds(*v2,*v3,&prod_v2_v3[0]);
00530         
00531         result = testRelErrors(
00532           num_rhs
00533           ,"<v4,v1>", &prod_v4_v1[0]
00534           ,"<v2,v3>", &prod_v2_v3[0]
00535           ,"adjoint_error_tol()", adjoint_error_tol()
00536           ,"adjoint_warning_tol()", adjoint_warning_tol()
00537           ,&*oss
00538           );
00539         if(!result) these_results = false;
00540         
00541       }
00542     }
00543     else {
00544       *oss << endl << "Adjoint operator not supported, skipping check!\n";
00545     }
00546 
00547     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00548 
00549   }
00550   else {
00551     if(out.get()) *out << endl << "this->check_adjoint()==false: Skipping check for the agreement of the adjoint and forward operators!\n";
00552   }
00553 
00554   if( check_for_symmetry() ) {
00555 
00556     if(out.get()) *out << endl << "this->check_for_symmetry()==true: Performing check of symmetry ... ";
00557 
00558 
00559     std::ostringstream ossStore;
00560     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00561     if(out.get()) ossStore.copyfmt(*out);
00562     bool these_results = true;
00563 
00564     SymmetricLinearOpTester<RangeScalar,DomainScalar>::checkSymmetry(
00565       op,&*dRand,*oss,num_rhs,num_random_vectors(),verbLevel,dump_all(),symmetry_error_tol(),symmetry_warning_tol(),&these_results
00566       );
00567     
00568     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00569     
00570   }
00571   else {
00572     if(out.get()) *out << endl << "this->check_for_symmetry()==false: Skipping check of symmetry ...\n";
00573   }
00574   
00575   if(out.get()) {
00576     if(success)
00577       *out << endl <<"Congratulations, this LinearOpBase object seems to check out!\n";
00578     else
00579       *out << endl <<"Oh no, at least one of the tests performed with this LinearOpBase object failed (see above failures)!\n";
00580     *out << endl << "*** Leaving LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::check(...)\n";
00581   }
00582 
00583   return success;
00584 }
00585 
00586 
00587 template<class RangeScalar, class DomainScalar>
00588 bool LinearOpTester<RangeScalar,DomainScalar>::check(
00589   const LinearOpBase<RangeScalar,DomainScalar>  &op
00590   ,Teuchos::FancyOStream                        *out
00591   ) const
00592 {
00593   return check(op,NULL,NULL,out);
00594 }
00595 
00596 template<class RangeScalar, class DomainScalar>
00597 bool LinearOpTester<RangeScalar,DomainScalar>::compare(
00598   const LinearOpBase<RangeScalar,DomainScalar>  &op1
00599   ,const LinearOpBase<RangeScalar,DomainScalar> &op2
00600   ,MultiVectorRandomizerBase<DomainScalar>      *domainRandomizer
00601   ,Teuchos::FancyOStream                        *out_arg
00602   ) const
00603 {
00604 
00605   using std::endl;
00606   using Teuchos::arrayArg;
00607   using Teuchos::FancyOStream;
00608   using Teuchos::OSTab;
00609   typedef Teuchos::ScalarTraits<RangeScalar>  RST;
00610   typedef Teuchos::ScalarTraits<DomainScalar> DST;
00611   bool success = true, result;
00612   const int num_rhs = this->num_rhs();
00613   const RangeScalar  r_half = RangeScalar(0.5)*RST::one();
00614   Teuchos::RCP<FancyOStream> out = Teuchos::rcp(out_arg,false);
00615   const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00616 
00617   OSTab tab(out,1,"THYRA");
00618 
00619   if(out.get()) {
00620     *out
00621       << endl << "*** Entering LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::compare(op1,op2,...) ...\n";
00622     if(show_all_tests())
00623       *out << endl << "describe op1:\n" << Teuchos::describe(op1,verbLevel);
00624     else
00625       *out << endl << "describe op1: " << op1.description() << endl;
00626     if(show_all_tests())
00627       *out << endl << "describe op2:\n" << Teuchos::describe(op2,verbLevel);
00628     else
00629       *out << endl << "describe op2: " << op2.description() << endl;
00630   }
00631 
00632   Teuchos::RCP< MultiVectorRandomizerBase<DomainScalar> > dRand;
00633   if(domainRandomizer)  dRand = Teuchos::rcp(domainRandomizer,false);
00634   else                  dRand = Teuchos::rcp(new UniversalMultiVectorRandomizer<DomainScalar>());
00635 
00636   Teuchos::RCP<const VectorSpaceBase<RangeScalar> >  range  = op1.range();
00637   Teuchos::RCP<const VectorSpaceBase<DomainScalar> > domain = op1.domain();
00638 
00639   if(out.get()) *out << endl << "Checking that range and domain spaces are compatible ... ";
00640 
00641   {
00642 
00643     std::ostringstream ossStore;
00644     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00645     if(out.get()) ossStore.copyfmt(*out);
00646     bool these_results = true;
00647 
00648     *oss << endl << "op1.domain()->isCompatible(*op2.domain()) ? ";
00649     result = op1.domain()->isCompatible(*op2.domain());
00650     if(!result) these_results = false;
00651     *oss << passfail(result) << endl;
00652     
00653     *oss << endl << "op1.range()->isCompatible(*op2.range()) ? ";
00654     result = op1.range()->isCompatible(*op2.range());
00655     if(!result) these_results = false;
00656     *oss << passfail(result) << endl;
00657 
00658     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00659 
00660   }
00661 
00662   if(!success) {
00663     if(out.get()) *out << endl << "Skipping further checks since operators are not compatible!\n";
00664     return success;
00665   }
00666 
00667   if(out.get()) *out << endl << "Checking that op1 == op2 ... ";
00668 
00669   {
00670 
00671     std::ostringstream ossStore;
00672     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00673     if(out.get()) ossStore.copyfmt(*out);
00674     bool these_results = true;
00675 
00676     *oss
00677       << endl << "Checking that op1 and op2 produce the same results:\n"
00678       << endl << "  0.5*op1*v1 == 0.5*op2*v1"
00679       << endl << "  \\________/    \\________/"
00680       << endl << "      v2            v3"
00681       << endl << ""
00682       << endl << "   |sum(v2)| == |sum(v3)|"
00683       << endl;
00684 
00685     for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00686       
00687       *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00688 
00689       OSTab tab(oss);
00690       
00691       if(dump_all()) *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00692       Teuchos::RCP<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00693       dRand->randomize(&*v1);
00694       if(dump_all()) *oss << endl << "v1 =\n" << *v1;
00695       
00696       if(dump_all()) *oss << endl << "v2 = 0.5*op1*v1 ...\n" ;
00697       Teuchos::RCP<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00698       apply( op1, NONCONJ_ELE, *v1, &*v2, r_half );
00699       if(dump_all()) *oss << endl << "v2 =\n" << *v2;
00700       
00701       if(dump_all()) *oss << endl << "v3 = 0.5*op2*v1 ...\n" ;
00702       Teuchos::RCP<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00703       apply( op2, NONCONJ_ELE, *v1, &*v3, r_half );
00704       if(dump_all()) *oss << endl << "v3 =\n" << *v3;
00705       
00706       std::vector<RangeScalar> sum_v2(num_rhs), sum_v3(num_rhs);
00707       sums(*v2,&sum_v2[0]);
00708       sums(*v3,&sum_v3[0]);
00709       
00710       result = testRelErrors(
00711         num_rhs
00712         ,"sum(v2)", &sum_v2[0]
00713         ,"sum(v3)", &sum_v3[0]
00714         ,"linear_properties_error_tol()", linear_properties_error_tol()
00715         ,"linear_properties_warning_tol()", linear_properties_warning_tol()
00716         ,&*oss
00717         );
00718       if(!result) these_results = false;
00719       
00720     }
00721 
00722     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00723 
00724   }
00725   
00726   if(out.get()) {
00727     if(success)
00728       *out << endl <<"Congratulations, these two LinearOpBase objects seem to be the same!\n";
00729     else
00730       *out << endl <<"Oh no, these two LinearOpBase objects seem to be different (see above failures)!\n";
00731     *out << endl << "*** Leaving LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::compare(...)\n";
00732   }
00733 
00734   return success;
00735 
00736 }
00737 
00738 template<class RangeScalar, class DomainScalar>
00739 bool LinearOpTester<RangeScalar,DomainScalar>::compare(
00740   const LinearOpBase<RangeScalar,DomainScalar>  &op1
00741   ,const LinearOpBase<RangeScalar,DomainScalar> &op2
00742   ,Teuchos::FancyOStream                        *out
00743   ) const
00744 {
00745   return compare(op1,op2,NULL,out);
00746 }
00747 
00748 } // namespace Thyra
00749 
00750 #endif // THYRA_LINEAR_OP_TESTER_HPP

Generated on Tue Oct 20 12:47:26 2009 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.4.7