Thyra_LinearOpTester.hpp

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       bool result;
00090       using Teuchos::OSTab;
00091       typedef Teuchos::ScalarTraits<Scalar> ST;
00092       const Scalar half = Scalar(0.4)*ST::one();
00093       Teuchos::RefCountPtr<const VectorSpaceBase<Scalar> > domain = op.domain();
00094       
00095       oss << endl << "op.domain()->isCompatible(*op.range()) == true : ";
00096       result = op.domain()->isCompatible(*op.range());
00097       if(!result) *these_results = false;
00098       oss << passfail(result) << endl;
00099       
00100       if(result) {
00101         
00102         oss
00103           << endl << "Checking that the operator is symmetric as:\n"
00104           << endl << "  <0.5*op*v2,v1> == <v2,0.5*op*v1>"
00105           << endl << "   \\_______/            \\_______/"
00106           << endl << "      v4                    v3"
00107           << endl << ""
00108           << endl << "         <v4,v1> == <v2,v3>"
00109           << endl;
00110         
00111         for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors; ++rand_vec_i ) {
00112           
00113           oss << endl << "Random vector tests = " << rand_vec_i << endl;
00114 
00115           OSTab tab(Teuchos::rcp(&oss,false));
00116 
00117           if(dump_all) oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00118           Teuchos::RefCountPtr<MultiVectorBase<Scalar> > v1 = createMembers(domain,num_rhs);
00119           dRand->randomize(&*v1);
00120           if(dump_all) oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00121           
00122           if(dump_all) oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00123           Teuchos::RefCountPtr<MultiVectorBase<Scalar> > v2 = createMembers(domain,num_rhs);
00124           dRand->randomize(&*v2);
00125           if(dump_all) oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00126           
00127           if(dump_all) oss << endl << "v3 = 0.5*op*v1 ...\n" ;
00128           Teuchos::RefCountPtr<MultiVectorBase<Scalar> > v3 = createMembers(domain,num_rhs);
00129           apply( op, NONCONJ_ELE, *v1, &*v3, half );
00130          if(dump_all) oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00131           
00132           if(dump_all) oss << endl << "v4 = 0.5*op*v2 ...\n" ;
00133           Teuchos::RefCountPtr<MultiVectorBase<Scalar> > v4 = createMembers(domain,num_rhs);
00134           apply( op, NONCONJ_ELE, *v2, &*v4, half );
00135           if(dump_all) oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00136 
00137           std::vector<Scalar> prod1(num_rhs), prod2(num_rhs);
00138           domain->scalarProds(*v4,*v1,&prod1[0]);
00139           domain->scalarProds(*v2,*v3,&prod2[0]);
00140           
00141           result = testRelErrors(
00142             num_rhs
00143             ,"<v4,v1>", &prod1[0]
00144             ,"<v2,v3>", &prod2[0]
00145             ,"symmetry_error_tol()", symmetry_error_tol
00146             ,"symmetry_warning_tol()", symmetry_warning_tol
00147             ,&oss
00148             );
00149           if(!result) *these_results = false;
00150         
00151         }
00152       }
00153       else {
00154         oss << endl << "Range and domain spaces are different, skipping check!\n";
00155       }
00156     }
00157 };
00158 
00159 // LinearOpTester
00160 
00161 template<class RangeScalar, class DomainScalar>
00162 LinearOpTester<RangeScalar,DomainScalar>::LinearOpTester(
00163   const bool          check_linear_properties
00164   ,const ScalarMag    linear_properties_warning_tol
00165   ,const ScalarMag    linear_properties_error_tol
00166   ,const bool         check_adjoint
00167   ,const ScalarMag    adjoint_warning_tol
00168   ,const ScalarMag    adjoint_error_tol
00169   ,const bool         check_for_symmetry
00170   ,const ScalarMag    symmetry_warning_tol
00171   ,const ScalarMag    symmetry_error_tol
00172   ,const int          num_random_vectors
00173   ,const bool         show_all_tests
00174   ,const bool         dump_all
00175   ,const int          num_rhs
00176   )
00177   :check_linear_properties_(check_linear_properties)
00178   ,linear_properties_warning_tol_(linear_properties_warning_tol)
00179   ,linear_properties_error_tol_(linear_properties_error_tol)
00180   ,check_adjoint_(check_adjoint)
00181   ,adjoint_warning_tol_(adjoint_warning_tol)
00182   ,adjoint_error_tol_(adjoint_error_tol)
00183   ,check_for_symmetry_(check_for_symmetry)
00184   ,symmetry_warning_tol_(symmetry_warning_tol)
00185   ,symmetry_error_tol_(symmetry_error_tol)
00186   ,num_random_vectors_(num_random_vectors)
00187   ,show_all_tests_(show_all_tests)
00188   ,dump_all_(dump_all)
00189   ,num_rhs_(num_rhs)
00190 {}
00191 
00192 template<class RangeScalar, class DomainScalar>
00193 void LinearOpTester<RangeScalar,DomainScalar>::enable_all_tests( const bool enable_all_tests )
00194 {
00195   check_linear_properties_ = enable_all_tests;
00196   check_adjoint_           = enable_all_tests;
00197   check_for_symmetry_      = enable_all_tests;
00198 }
00199 
00200 template<class RangeScalar, class DomainScalar>
00201 void LinearOpTester<RangeScalar,DomainScalar>::set_all_warning_tol( const ScalarMag warning_tol )
00202 {
00203   linear_properties_warning_tol_  = warning_tol;
00204   adjoint_warning_tol_            = warning_tol;
00205   symmetry_warning_tol_           = warning_tol;
00206 }
00207 
00208 template<class RangeScalar, class DomainScalar>
00209 void LinearOpTester<RangeScalar,DomainScalar>::set_all_error_tol( const ScalarMag error_tol )
00210 {
00211   linear_properties_error_tol_  = error_tol;
00212   adjoint_error_tol_            = error_tol;
00213   symmetry_error_tol_           = error_tol;
00214 }
00215 
00216 template<class RangeScalar, class DomainScalar>
00217 bool LinearOpTester<RangeScalar,DomainScalar>::check(
00218   const LinearOpBase<RangeScalar,DomainScalar>  &op
00219   ,MultiVectorRandomizerBase<RangeScalar>       *rangeRandomizer
00220   ,MultiVectorRandomizerBase<DomainScalar>      *domainRandomizer
00221   ,Teuchos::FancyOStream                        *out_arg
00222   ) const
00223 {
00224 
00225   using std::endl;
00226   using Teuchos::FancyOStream;
00227   using Teuchos::OSTab;
00228   typedef Teuchos::ScalarTraits<RangeScalar>  RST;
00229   typedef Teuchos::ScalarTraits<DomainScalar> DST;
00230   bool success = true, result;
00231   const int          num_rhs = this->num_rhs();
00232   const RangeScalar  r_one  = RST::one();
00233   const DomainScalar d_one  = DST::one();
00234   const RangeScalar  r_half = RangeScalar(0.5)*r_one;
00235   const DomainScalar d_half = DomainScalar(0.5)*d_one;
00236   Teuchos::RefCountPtr<FancyOStream> out = Teuchos::rcp(out_arg,false);
00237   const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00238 
00239   OSTab tab(out,1,"THYRA");
00240 
00241   // ToDo 04/28/2005:
00242   // * Test the MultiVectorBase apply() function and output to the VectorBase apply() function!
00243 
00244   if(out.get()) {
00245     *out << endl << "*** Entering LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::check(op,...) ...\n";
00246     if(show_all_tests()) {
00247       *out << endl << "describe op:\n" << Teuchos::describe(op,verbLevel);
00248 /*
00249       if(op.applyTransposeSupports(CONJ_ELE) && verbLevel==Teuchos::VERB_EXTREME) {
00250         *out << endl << "describe adjoint op:\n";
00251         describeLinearOp(*adjoint(Teuchos::rcp(&op,false)),*out,verbLevel);
00252       }
00253 */
00254     }
00255     else {
00256       *out << endl << "describe op:\n" << Teuchos::describe(op,Teuchos::VERB_LOW);
00257     }
00258   }
00259 
00260   Teuchos::RefCountPtr< MultiVectorRandomizerBase<RangeScalar> >  rRand;
00261   if(rangeRandomizer)   rRand = Teuchos::rcp(rangeRandomizer,false);
00262   else                  rRand = Teuchos::rcp(new UniversalMultiVectorRandomizer<RangeScalar>());
00263   Teuchos::RefCountPtr< MultiVectorRandomizerBase<DomainScalar> > dRand;
00264   if(domainRandomizer)  dRand = Teuchos::rcp(domainRandomizer,false);
00265   else                  dRand = Teuchos::rcp(new UniversalMultiVectorRandomizer<DomainScalar>());
00266   
00267   if(out.get())
00268     *out << endl << "Checking the domain and range spaces ... ";
00269 
00270   Teuchos::RefCountPtr<const VectorSpaceBase<RangeScalar> >  range  = op.range();
00271   Teuchos::RefCountPtr<const VectorSpaceBase<DomainScalar> > domain = op.domain();
00272   
00273   if(1) {
00274 
00275     std::ostringstream ossStore;
00276     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00277     if(out.get()) ossStore.copyfmt(*out);
00278     bool these_results = true;
00279 
00280     *oss << endl << "op.domain().get() != NULL ? ";
00281     result = domain.get() != NULL;
00282     if(!result) these_results = false;
00283     *oss << passfail(result) << endl;
00284     
00285     *oss << endl << "op.range().get() != NULL ? ";
00286     result = range.get() != NULL;
00287     if(!result) these_results = false;
00288     *oss << passfail(result) << endl;
00289 
00290     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00291 
00292   }
00293 
00294   if( check_linear_properties() ) {
00295 
00296     if(out.get()) *out << endl << "this->check_linear_properties()==true: Checking the linear properties of the forward linear operator ... ";
00297 
00298     std::ostringstream ossStore;
00299     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00300     if(out.get()) ossStore.copyfmt(*out);
00301     bool these_results = true;
00302 
00303     *oss << endl << "op.applySupports(NONCONJ_ELE) == true ? ";
00304     result = op.applySupports(NONCONJ_ELE);
00305     if(!result) these_results = false;
00306     *oss << passfail(result) << endl;
00307 
00308     if(result) {
00309     
00310       *oss
00311         << endl << "Checking that the forward operator is truly linear:\n"
00312         << endl << "  0.5*op*(v1 + v2) == 0.5*op*v1 + 0.5*op*v2"
00313         << endl << "          \\_____/         \\___/"
00314         << endl << "             v3            v5"
00315         << endl << "  \\_____________/     \\___________________/"
00316         << endl << "         v4                    v5"
00317         << endl << ""
00318         << endl << "           sum(v4) == sum(v5)"
00319         << endl;
00320       
00321       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00322         
00323         *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00324 
00325         OSTab tab(oss);
00326         
00327         *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00328         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00329         dRand->randomize(&*v1);
00330         if(dump_all()) *oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00331         
00332         *oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00333         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v2 = createMembers(domain,num_rhs);
00334         dRand->randomize(&*v2);
00335         if(dump_all()) *oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00336         
00337         *oss << endl << "v3 = v1 + v2 ...\n" ;
00338         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v3 = createMembers(domain,num_rhs);
00339         V_VpV(&*v3,*v1,*v2);
00340         if(dump_all()) *oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00341         
00342         *oss << endl << "v4 = 0.5*op*v3 ...\n" ;
00343         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v4 = createMembers(range,num_rhs);
00344         apply( op, NONCONJ_ELE, *v3, &*v4, r_half );
00345         if(dump_all()) *oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00346         
00347         *oss << endl << "v5 = op*v1 ...\n" ;
00348         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v5 = createMembers(range,num_rhs);
00349         apply( op, NONCONJ_ELE, *v1, &*v5 );
00350         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00351         
00352         *oss << endl << "v5 = 0.5*op*v2 + 0.5*v5 ...\n" ;
00353         apply( op, NONCONJ_ELE, *v2, &*v5, r_half, r_half );
00354         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00355 
00356         std::vector<RangeScalar> sum_v4(num_rhs), sum_v5(num_rhs);
00357         sums(*v4,&sum_v4[0]);
00358         sums(*v5,&sum_v5[0]);
00359         
00360         result = testRelErrors(
00361           num_rhs
00362           ,"sum(v4)", &sum_v4[0]
00363           ,"sum(v5)", &sum_v5[0]
00364           ,"linear_properties_error_tol()", linear_properties_error_tol()
00365           ,"linear_properties_warning_tol()", linear_properties_warning_tol()
00366           ,&*oss
00367           );
00368         if(!result) these_results = false;
00369         
00370       }
00371     }
00372     else {
00373       *oss << endl << "Forward operator not supported, skipping check!\n";
00374     }
00375 
00376     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00377 
00378   }
00379   else {
00380     if(out.get()) *out << endl << "this->check_linear_properties()==false: Skipping the check of the linear properties of the forward operator!\n";
00381   }
00382 
00383   if( check_linear_properties() && check_adjoint() ) {
00384 
00385     if(out.get()) *out << endl << "(this->check_linear_properties()&&this->check_adjoint())==true: Checking the linear properties of the adjoint operator ... ";
00386 
00387     std::ostringstream ossStore;
00388     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00389     if(out.get()) ossStore.copyfmt(*out);
00390     bool these_results = true;
00391 
00392     *oss << endl << "op.applyTransposeSupports(CONJ_ELE) == true ? ";
00393     result = op.applyTransposeSupports(CONJ_ELE);
00394     if(!result) these_results = false;
00395     *oss << passfail(result) << endl;
00396 
00397     if(result) {
00398     
00399       *oss
00400         << endl << "Checking that the adjoint operator is truly linear:\n"
00401         << endl << "  0.5*op'*(v1 + v2) == 0.5*op'*v1 + 0.5*op'*v2"
00402         << endl << "           \\_____/         \\____/"
00403         << endl << "              v3             v5"
00404         << endl << "  \\_______________/    \\_____________________/"
00405         << endl << "         v4                      v5"
00406         << endl << ""
00407         << endl << "           sum(v4) == sum(v5)"
00408         << endl;
00409       
00410       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00411         
00412         *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00413 
00414         OSTab tab(oss);
00415         
00416         *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00417         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v1 = createMembers(range,num_rhs);
00418         rRand->randomize(&*v1);
00419         if(dump_all()) *oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00420         
00421         *oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00422         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00423         rRand->randomize(&*v2);
00424         if(dump_all()) *oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00425         
00426         *oss << endl << "v3 = v1 + v2 ...\n" ;
00427         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00428         V_VpV(&*v3,*v1,*v2);
00429         if(dump_all()) *oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00430         
00431         *oss << endl << "v4 = 0.5*op'*v3 ...\n" ;
00432         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v4 = createMembers(domain,num_rhs);
00433         applyTranspose( op, CONJ_ELE, *v3, &*v4, d_half );
00434         if(dump_all()) *oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00435         
00436         *oss << endl << "v5 = op'*v1 ...\n" ;
00437         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v5 = createMembers(domain,num_rhs);
00438         applyTranspose( op, CONJ_ELE, *v1, &*v5 );
00439         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00440         
00441         *oss << endl << "v5 = 0.5*op'*v2 + 0.5*v5 ...\n" ;
00442         applyTranspose( op, CONJ_ELE, *v2, &*v5, d_half, d_half );
00443         if(dump_all()) *oss << endl << "v5 =\n" << describe(*v5,verbLevel);
00444         
00445 
00446         std::vector<DomainScalar> sum_v4(num_rhs), sum_v5(num_rhs);
00447         sums(*v4,&sum_v4[0]);
00448         sums(*v5,&sum_v5[0]);
00449         
00450         result = testRelErrors(
00451           num_rhs
00452           ,"sum(v4)", &sum_v4[0]
00453           ,"sum(v5)", &sum_v5[0]
00454           ,"linear_properties_error_tol()", linear_properties_error_tol()
00455           ,"linear_properties_warning_tol()", linear_properties_warning_tol()
00456           ,&*oss
00457           );
00458         if(!result) these_results = false;
00459         
00460       }
00461     }
00462     else {
00463       *oss << endl << "Adjoint operator not supported, skipping check!\n";
00464     }
00465 
00466     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00467 
00468   }
00469   else {
00470     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";
00471   }
00472   
00473   if( check_adjoint() ) {
00474 
00475     if(out.get()) *out << endl << "this->check_adjoint()==true: Checking the agreement of the adjoint and forward operators ... ";
00476 
00477     std::ostringstream ossStore;
00478     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00479     if(out.get()) ossStore.copyfmt(*out);
00480     bool these_results = true;
00481     
00482     *oss << endl << "op.applyTransposeSupports(CONJ_ELE) == true ? ";
00483     result = op.applyTransposeSupports(CONJ_ELE);
00484     if(!result) these_results = false;
00485     *oss << passfail(result) << endl;
00486 
00487     if(result) {
00488     
00489       *oss
00490         << endl << "Checking that the adjoint agrees with the non-adjoint operator as:\n"
00491         << endl << "  <0.5*op'*v2,v1> == <v2,0.5*op*v1>"
00492         << endl << "   \\________/            \\_______/"
00493         << endl << "       v4                   v3"
00494         << endl << ""
00495         << endl << "         <v4,v1>  == <v2,v3>"
00496         << endl;
00497     
00498       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00499       
00500         *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00501 
00502         OSTab tab(oss);
00503       
00504         *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00505         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00506         dRand->randomize(&*v1);
00507         if(dump_all()) *oss << endl << "v1 =\n" << describe(*v1,verbLevel);
00508       
00509         *oss << endl << "v2 = randomize(-1,+1); ...\n" ;
00510         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00511         rRand->randomize(&*v2);
00512         if(dump_all()) *oss << endl << "v2 =\n" << describe(*v2,verbLevel);
00513       
00514         *oss << endl << "v3 = 0.5*op*v1 ...\n" ;
00515         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00516         apply( op, NONCONJ_ELE, *v1, &*v3, r_half );
00517         if(dump_all()) *oss << endl << "v3 =\n" << describe(*v3,verbLevel);
00518       
00519         *oss << endl << "v4 = 0.5*op'*v2 ...\n" ;
00520         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v4 = createMembers(domain,num_rhs);
00521         applyTranspose( op, CONJ_ELE, *v2, &*v4, d_half );
00522         if(dump_all()) *oss << endl << "v4 =\n" << describe(*v4,verbLevel);
00523 
00524         std::vector<DomainScalar> prod_v4_v1(num_rhs);
00525         domain->scalarProds(*v4,*v1,&prod_v4_v1[0]);
00526         std::vector<RangeScalar> prod_v2_v3(num_rhs);
00527         range->scalarProds(*v2,*v3,&prod_v2_v3[0]);
00528         
00529         result = testRelErrors(
00530           num_rhs
00531           ,"<v4,v1>", &prod_v4_v1[0]
00532           ,"<v2,v3>", &prod_v2_v3[0]
00533           ,"adjoint_error_tol()", adjoint_error_tol()
00534           ,"adjoint_warning_tol()", adjoint_warning_tol()
00535           ,&*oss
00536           );
00537         if(!result) these_results = false;
00538         
00539       }
00540     }
00541     else {
00542       *oss << endl << "Adjoint operator not supported, skipping check!\n";
00543     }
00544 
00545     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00546 
00547   }
00548   else {
00549     if(out.get()) *out << endl << "this->check_adjoint()==false: Skipping check for the agreement of the adjoint and forward operators!\n";
00550   }
00551 
00552   if( check_for_symmetry() ) {
00553 
00554     if(out.get()) *out << endl << "this->check_for_symmetry()==true: Performing check of symmetry ... ";
00555 
00556 
00557     std::ostringstream ossStore;
00558     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00559     if(out.get()) ossStore.copyfmt(*out);
00560     bool these_results = true;
00561 
00562     SymmetricLinearOpTester<RangeScalar,DomainScalar>::checkSymmetry(
00563       op,&*dRand,*oss,num_rhs,num_random_vectors(),verbLevel,dump_all(),symmetry_error_tol(),symmetry_warning_tol(),&these_results
00564       );
00565     
00566     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00567     
00568   }
00569   else {
00570     if(out.get()) *out << endl << "this->check_for_symmetry()==false: Skipping check of symmetry ...\n";
00571   }
00572   
00573   if(out.get()) {
00574     if(success)
00575       *out << endl <<"Congratulations, this LinearOpBase object seems to check out!\n";
00576     else
00577       *out << endl <<"Oh no, at least one of the tests performed with this LinearOpBase object failed (see above failures)!\n";
00578     *out << endl << "*** Leaving LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::check(...)\n";
00579   }
00580 
00581   return success;
00582 }
00583 
00584 
00585 template<class RangeScalar, class DomainScalar>
00586 bool LinearOpTester<RangeScalar,DomainScalar>::check(
00587   const LinearOpBase<RangeScalar,DomainScalar>  &op
00588   ,Teuchos::FancyOStream                        *out
00589   ) const
00590 {
00591   return check(op,NULL,NULL,out);
00592 }
00593 
00594 template<class RangeScalar, class DomainScalar>
00595 bool LinearOpTester<RangeScalar,DomainScalar>::compare(
00596   const LinearOpBase<RangeScalar,DomainScalar>  &op1
00597   ,const LinearOpBase<RangeScalar,DomainScalar> &op2
00598   ,MultiVectorRandomizerBase<DomainScalar>      *domainRandomizer
00599   ,Teuchos::FancyOStream                        *out_arg
00600   ) const
00601 {
00602 
00603   using std::endl;
00604   using Teuchos::arrayArg;
00605   using Teuchos::FancyOStream;
00606   using Teuchos::OSTab;
00607   typedef Teuchos::ScalarTraits<RangeScalar>  RST;
00608   typedef Teuchos::ScalarTraits<DomainScalar> DST;
00609   bool success = true, result;
00610   const int num_rhs = this->num_rhs();
00611   const RangeScalar  r_half = RangeScalar(0.5)*RST::one();
00612   Teuchos::RefCountPtr<FancyOStream> out = Teuchos::rcp(out_arg,false);
00613   const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00614 
00615   OSTab tab(out,1,"THYRA");
00616 
00617   if(out.get()) {
00618     *out
00619       << endl << "*** Entering LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::compare(op1,op2,...) ...\n";
00620     if(show_all_tests())
00621       *out << endl << "describe op1:\n" << Teuchos::describe(op1,verbLevel);
00622     else
00623       *out << endl << "describe op1: " << op1.description() << endl;
00624     if(show_all_tests())
00625       *out << endl << "describe op2:\n" << Teuchos::describe(op2,verbLevel);
00626     else
00627       *out << endl << "describe op2: " << op2.description() << endl;
00628   }
00629 
00630   Teuchos::RefCountPtr< MultiVectorRandomizerBase<DomainScalar> > dRand;
00631   if(domainRandomizer)  dRand = Teuchos::rcp(domainRandomizer,false);
00632   else                  dRand = Teuchos::rcp(new UniversalMultiVectorRandomizer<DomainScalar>());
00633 
00634   Teuchos::RefCountPtr<const VectorSpaceBase<RangeScalar> >  range  = op1.range();
00635   Teuchos::RefCountPtr<const VectorSpaceBase<DomainScalar> > domain = op1.domain();
00636 
00637   if(out.get()) *out << endl << "Checking that range and domain spaces are compatible ... ";
00638 
00639   if(1) {
00640 
00641     std::ostringstream ossStore;
00642     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00643     if(out.get()) ossStore.copyfmt(*out);
00644     bool these_results = true;
00645 
00646     *oss << endl << "op1.domain()->isCompatible(*op2.domain()) ? ";
00647     result = op1.domain()->isCompatible(*op2.domain());
00648     if(!result) these_results = false;
00649     *oss << passfail(result) << endl;
00650     
00651     *oss << endl << "op1.range()->isCompatible(*op2.range()) ? ";
00652     result = op1.range()->isCompatible(*op2.range());
00653     if(!result) these_results = false;
00654     *oss << passfail(result) << endl;
00655 
00656     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00657 
00658   }
00659 
00660   if(!success) {
00661     if(out.get()) *out << endl << "Skipping further checks since operators are not compatible!\n";
00662     return success;
00663   }
00664 
00665   if(out.get()) *out << endl << "Checking that op1 == op2 ... ";
00666 
00667   if(1) {
00668 
00669     std::ostringstream ossStore;
00670     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00671     if(out.get()) ossStore.copyfmt(*out);
00672     bool these_results = true;
00673 
00674     *oss
00675       << endl << "Checking that op1 and op2 produce the same results:\n"
00676       << endl << "  0.5*op1*v1 == 0.5*op2*v1"
00677       << endl << "  \\________/    \\________/"
00678       << endl << "      v2            v3"
00679       << endl << ""
00680       << endl << "   |sum(v2)| == |sum(v3)|"
00681       << endl;
00682 
00683     for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00684       
00685       *oss << endl << "Random vector tests = " << rand_vec_i << endl;
00686 
00687       OSTab tab(oss);
00688       
00689       if(dump_all()) *oss << endl << "v1 = randomize(-1,+1); ...\n" ;
00690       Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00691       dRand->randomize(&*v1);
00692       if(dump_all()) *oss << endl << "v1 =\n" << *v1;
00693       
00694       if(dump_all()) *oss << endl << "v2 = 0.5*op1*v1 ...\n" ;
00695       Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00696       apply( op1, NONCONJ_ELE, *v1, &*v2, r_half );
00697       if(dump_all()) *oss << endl << "v2 =\n" << *v2;
00698       
00699       if(dump_all()) *oss << endl << "v3 = 0.5*op2*v1 ...\n" ;
00700       Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00701       apply( op2, NONCONJ_ELE, *v1, &*v3, r_half );
00702       if(dump_all()) *oss << endl << "v3 =\n" << *v3;
00703       
00704       std::vector<RangeScalar> sum_v2(num_rhs), sum_v3(num_rhs);
00705       sums(*v2,&sum_v2[0]);
00706       sums(*v3,&sum_v3[0]);
00707       
00708       result = testRelErrors(
00709         num_rhs
00710         ,"sum(v2)", &sum_v2[0]
00711         ,"sum(v3)", &sum_v3[0]
00712         ,"linear_properties_error_tol()", linear_properties_error_tol()
00713         ,"linear_properties_warning_tol()", linear_properties_warning_tol()
00714         ,&*oss
00715         );
00716       if(!result) these_results = false;
00717       
00718     }
00719 
00720     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00721 
00722   }
00723   
00724   if(out.get()) {
00725     if(success)
00726       *out << endl <<"Congratulations, these two LinearOpBase objects seem to be the same!\n";
00727     else
00728       *out << endl <<"Oh no, these two LinearOpBase objects seem to be different (see above failures)!\n";
00729     *out << endl << "*** Leaving LinearOpTester<"<<RST::name()<<","<<DST::name()<<">::compare(...)\n";
00730   }
00731 
00732   return success;
00733 
00734 }
00735 
00736 template<class RangeScalar, class DomainScalar>
00737 bool LinearOpTester<RangeScalar,DomainScalar>::compare(
00738   const LinearOpBase<RangeScalar,DomainScalar>  &op1
00739   ,const LinearOpBase<RangeScalar,DomainScalar> &op2
00740   ,Teuchos::FancyOStream                        *out
00741   ) const
00742 {
00743   return compare(op1,op2,NULL,out);
00744 }
00745 
00746 } // namespace Thyra
00747 
00748 #endif // THYRA_LINEAR_OP_TESTER_HPP

Generated on Thu Sep 18 12:32:31 2008 for Thyra Operator/Vector Support by doxygen 1.3.9.1