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