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