Thyra Version of the Day
Thyra_LinearOpWithSolveTester_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 Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 
00043 #ifndef THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
00044 #define THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
00045 
00046 
00047 #include "Thyra_LinearOpWithSolveTester_decl.hpp"
00048 #include "Thyra_LinearOpWithSolveBase.hpp"
00049 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00050 #include "Thyra_describeLinearOp.hpp"
00051 #include "Thyra_VectorStdOps.hpp"
00052 #include "Thyra_MultiVectorStdOps.hpp"
00053 #include "Thyra_TestingTools.hpp"
00054 #include "Teuchos_Time.hpp"
00055 #include "Teuchos_TestingHelpers.hpp"
00056 
00057 
00058 namespace Thyra {
00059 
00060 
00061 // Constructors/initializers
00062 
00063 
00064 template<class Scalar>
00065 LinearOpWithSolveTester<Scalar>::LinearOpWithSolveTester()
00066   :check_forward_default_(check_forward_default_default_),
00067    forward_default_residual_warning_tol_(warning_tol_default_),
00068    forward_default_residual_error_tol_(error_tol_default_),
00069    forward_default_solution_error_warning_tol_(warning_tol_default_),
00070    forward_default_solution_error_error_tol_(error_tol_default_),
00071    check_forward_residual_(check_forward_residual_default_),
00072    forward_residual_solve_tol_(solve_tol_default_),
00073    forward_residual_slack_warning_tol_(slack_warning_tol_default_),
00074    forward_residual_slack_error_tol_(slack_error_tol_default_),
00075    check_adjoint_default_(check_adjoint_default_default_),
00076    adjoint_default_residual_warning_tol_(warning_tol_default_),
00077    adjoint_default_residual_error_tol_(error_tol_default_),
00078    adjoint_default_solution_error_warning_tol_(warning_tol_default_),
00079    adjoint_default_solution_error_error_tol_(error_tol_default_),
00080    check_adjoint_residual_(check_adjoint_residual_default_),
00081    adjoint_residual_solve_tol_(solve_tol_default_),
00082    adjoint_residual_slack_warning_tol_(slack_warning_tol_default_),
00083    adjoint_residual_slack_error_tol_(slack_error_tol_default_),
00084    num_random_vectors_(num_random_vectors_default_),
00085    show_all_tests_(show_all_tests_default_),
00086    dump_all_(dump_all_default_),
00087    num_rhs_(num_rhs_default_)
00088 {}
00089 
00090 
00091 template<class Scalar>
00092 void LinearOpWithSolveTester<Scalar>::turn_off_all_tests()
00093 {
00094   check_forward_default_ = false;
00095   check_forward_residual_ = false;
00096   check_adjoint_default_ = false;
00097   check_adjoint_residual_ = false;
00098 }
00099 
00100 
00101 template<class Scalar>
00102 void
00103 LinearOpWithSolveTester<Scalar>::set_all_solve_tol(
00104   const ScalarMag solve_tol )
00105 {
00106   forward_residual_solve_tol_ = solve_tol;
00107   forward_residual_solve_tol_ = solve_tol;
00108   adjoint_residual_solve_tol_ = solve_tol;
00109 }
00110 
00111 
00112 template<class Scalar>
00113 void
00114 LinearOpWithSolveTester<Scalar>::set_all_slack_warning_tol(
00115   const ScalarMag slack_warning_tol )
00116 {
00117   forward_default_residual_warning_tol_ = slack_warning_tol;
00118   forward_default_solution_error_warning_tol_ = slack_warning_tol;
00119   forward_residual_slack_warning_tol_ = slack_warning_tol;
00120   adjoint_default_residual_warning_tol_ = slack_warning_tol;
00121   adjoint_default_solution_error_warning_tol_ = slack_warning_tol;
00122   adjoint_residual_slack_warning_tol_ = slack_warning_tol;
00123 }
00124 
00125 
00126 template<class Scalar>
00127 void
00128 LinearOpWithSolveTester<Scalar>::set_all_slack_error_tol(
00129   const ScalarMag slack_error_tol )
00130 {
00131   forward_default_residual_error_tol_ = slack_error_tol;
00132   forward_default_solution_error_error_tol_ = slack_error_tol;
00133   forward_residual_slack_error_tol_ = slack_error_tol;
00134   adjoint_default_residual_error_tol_ = slack_error_tol;
00135   adjoint_default_solution_error_error_tol_ = slack_error_tol;
00136   adjoint_residual_slack_error_tol_ = slack_error_tol;
00137 }
00138 
00139 
00140 // Overridden from ParameterListAcceptor
00141 
00142 
00143 template<class Scalar>
00144 void LinearOpWithSolveTester<Scalar>::setParameterList(
00145   const RCP<ParameterList>& paramList )
00146 {
00147   using Teuchos::getParameter;
00148   ParameterList &pl = *paramList;
00149   this->setMyParamList(paramList);
00150   paramList->validateParametersAndSetDefaults(*getValidParameters());
00151   set_all_solve_tol(getParameter<ScalarMag>(pl, AllSolveTol_name_));
00152   set_all_slack_warning_tol(getParameter<ScalarMag>(pl, AllSlackWarningTol_name_));
00153   set_all_slack_error_tol(getParameter<ScalarMag>(pl, AllSlackErrorTol_name_));
00154   show_all_tests(getParameter<bool>(pl, ShowAllTests_name_));
00155   dump_all(getParameter<bool>(pl, DumpAll_name_));
00156   // ToDo: Add more parameters as you need them
00157 }
00158 
00159 
00160 template<class Scalar>
00161 RCP<const ParameterList>
00162 LinearOpWithSolveTester<Scalar>::getValidParameters() const
00163 {
00164   static RCP<const ParameterList> validPL;
00165   if (is_null(validPL) ) {
00166     RCP<ParameterList> pl = Teuchos::parameterList();
00167     pl->set(AllSolveTol_name_, solve_tol_default_,
00168       "Sets all of the solve tolerances to the same value.  Note: This option\n"
00169       "is applied before any specific tolerance which may override it.");
00170     pl->set(AllSlackWarningTol_name_, slack_warning_tol_default_,
00171       "Sets all of the slack warning values to the same value.  Note: This option\n"
00172       "is applied before any specific tolerance which may override it.");
00173     pl->set(AllSlackErrorTol_name_, slack_error_tol_default_,
00174       "Sets all of the slack error values to the same value.  Note: This option\n"
00175       "is applied before any specific tolerance which may override it.");
00176     pl->set(ShowAllTests_name_, false,
00177       "If true, then all tests be traced to the output stream.");
00178     pl->set(DumpAll_name_, false,
00179       "If true, then all qualtities will be dumped to the output stream.  Warning!\n"
00180       "only do this to debug smaller problems as this can create a lot of output");
00181     // ToDo: Add more parameters as you need them (don't forget to test them)
00182     validPL = pl;
00183   }
00184   return validPL;
00185 }
00186 
00187 
00188 // LOWS testing
00189 
00190 
00191 template<class Scalar>
00192 bool LinearOpWithSolveTester<Scalar>::check(
00193   const LinearOpWithSolveBase<Scalar> &op,
00194   Teuchos::FancyOStream *out_arg ) const
00195 {
00196 
00197   using std::endl;
00198   using Teuchos::optInArg;
00199   using Teuchos::FancyOStream;
00200   using Teuchos::OSTab;
00201   typedef Teuchos::VerboseObjectTempState<LinearOpWithSolveBase<Scalar> > VOTS;
00202   typedef Teuchos::ScalarTraits<Scalar> ST;
00203   bool success = true, result;
00204   const int l_num_rhs = this->num_rhs();
00205   Teuchos::RCP<FancyOStream> out = Teuchos::rcp(out_arg,false);
00206   const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00207   Teuchos::Time timer("");
00208 
00209   OSTab tab(out,1,"THYRA");
00210 
00211   if(out.get()) {
00212     *out <<endl<< "*** Entering LinearOpWithSolveTester<"<<ST::name()<<">::check(op,...) ...\n";
00213     if(show_all_tests()) {
00214       *out <<endl<< "describe forward op:\n" << Teuchos::describe(op,verbLevel);
00215       if(opSupported(op, CONJTRANS) && verbLevel==Teuchos::VERB_EXTREME) {
00216         *out <<endl<< "describe adjoint op:\n";
00217         describeLinearOp<Scalar>(
00218           *adjoint(RCP<const LinearOpBase<Scalar> >(Teuchos::rcp(&op,false))),
00219           *out, verbLevel
00220           );
00221       }
00222     }
00223     else {
00224       *out <<endl<< "describe op: " << op.description() << endl;
00225     }
00226   }
00227   
00228   Teuchos::RCP<const VectorSpaceBase<Scalar> >   range  = op.range();
00229   Teuchos::RCP<const VectorSpaceBase<Scalar> >  domain = op.domain();
00230   
00231   if( check_forward_default() ) {
00232 
00233     if(out.get()) *out <<endl<< "this->check_forward_default()==true: Checking the default forward solve ... ";
00234 
00235     std::ostringstream ossStore;
00236     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00237     if(out.get()) ossStore.copyfmt(*out);
00238     bool these_results = true;
00239 
00240     result = true;
00241     TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(NOTRANS), true, *oss, result);
00242     if (!result) these_results = false;
00243 
00244     if(result) {
00245     
00246       *oss
00247         <<endl<< "Checking that the forward default solve matches the forward operator:\n"
00248         <<endl<< "  inv(Op)*Op*v1 == v1"
00249         <<endl<< "          \\___/"
00250         <<endl<< "           v2"
00251         <<endl<< "  \\___________/"
00252         <<endl<< "         v3"
00253         <<endl<< ""
00254         <<endl<< "  v4 = v3-v1"
00255         <<endl<< "  v5 = Op*v3-v2"
00256         <<endl<< ""
00257         <<endl<< "  norm(v4)/norm(v1) <= forward_default_solution_error_error_tol()"
00258         <<endl<< "  norm(v5)/norm(v2) <= forward_default_residual_error_tol()"
00259         <<endl;
00260 
00261       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00262 
00263         OSTab tab2(oss);
00264 
00265         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00266 
00267         tab.incrTab();
00268       
00269         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00270         Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,l_num_rhs);
00271         Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*v1 );
00272         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00273       
00274         *oss <<endl<< "v2 = Op*v1 ...\n" ;
00275         Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,l_num_rhs);
00276         timer.start(true);
00277         Thyra::apply(op, NOTRANS, *v1, v2.ptr());
00278         timer.stop();
00279         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00280         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00281 
00282         *oss <<endl<< "v3 = inv(Op)*v2 ...\n" ;
00283         Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,l_num_rhs);
00284         assign(&*v3, ST::zero());
00285         SolveStatus<Scalar> solveStatus;
00286         {
00287           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00288           timer.start(true);
00289           solveStatus = solve<Scalar>(op, NOTRANS, *v2, v3.ptr());
00290           timer.stop();
00291           OSTab(oss).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00292         }
00293         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00294         *oss
00295           <<endl<< "solve status:\n";
00296         OSTab(oss).o() << solveStatus;
00297 
00298         *oss <<endl<< "v4 = v3 - v1 ...\n" ;
00299         Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,l_num_rhs);
00300         V_VmV( &*v4, *v3, *v1 );
00301         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00302       
00303         *oss <<endl<< "v5 = Op*v3 - v2 ...\n" ;
00304         Teuchos::RCP<MultiVectorBase<Scalar> > v5 = createMembers(range,l_num_rhs);
00305         assign( &*v5, *v2 );
00306         timer.start(true);
00307         Thyra::apply(op, NOTRANS, *v3, v5.ptr(), Scalar(1.0) ,Scalar(-1.0));
00308         timer.stop();
00309         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00310         if(dump_all()) *oss <<endl<< "v5 =\n" << describe(*v5,verbLevel);
00311 
00312         std::vector<ScalarMag> norms_v1(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v1(l_num_rhs);
00313         norms(*v1,&norms_v1[0]);
00314         norms(*v4,&norms_v4[0]);
00315         std::transform(
00316           norms_v4.begin(),norms_v4.end(),norms_v1.begin(),norms_v4_norms_v1.begin()
00317           ,std::divides<ScalarMag>()
00318           );
00319 
00320         result = testMaxErrors(
00321           l_num_rhs, "norm(v4)/norm(v1)", &norms_v4_norms_v1[0]
00322           ,"forward_default_solution_error_error_tol()", forward_default_solution_error_error_tol()
00323           ,"forward_default_solution_error_warning_tol()", forward_default_solution_error_warning_tol()
00324           ,&*oss
00325           );
00326         if(!result) these_results = false;
00327 
00328         std::vector<ScalarMag> norms_v2(l_num_rhs), norms_v5(l_num_rhs), norms_v5_norms_v2(l_num_rhs);
00329         norms(*v2,&norms_v2[0]);
00330         norms(*v5,&norms_v5[0]);
00331         std::transform(
00332           norms_v5.begin(),norms_v5.end(),norms_v2.begin(),norms_v5_norms_v2.begin()
00333           ,std::divides<ScalarMag>()
00334           );
00335 
00336         result = testMaxErrors(
00337           l_num_rhs, "norm(v5)/norm(v2)", &norms_v5_norms_v2[0]
00338           ,"forward_default_residual_error_tol()", forward_default_residual_error_tol()
00339           ,"forward_default_residual_warning_tol()", forward_default_residual_warning_tol()
00340           ,&*oss
00341           );
00342         if(!result) these_results = false;
00343         
00344       }
00345     }
00346     else {
00347       *oss <<endl<< "Forward operator not supported, skipping check!\n";
00348     }
00349 
00350     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00351 
00352   }
00353   else {
00354     if(out.get()) *out <<endl<< "this->check_forward_default()==false: Skipping the check of the default forward solve!\n";
00355   }
00356   
00357   if( check_forward_residual() ) {
00358 
00359     if(out.get()) *out <<endl<< "this->check_forward_residual()==true: Checking the forward solve with a tolerance on the residual ... ";
00360 
00361     std::ostringstream ossStore;
00362     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00363     if(out.get()) ossStore.copyfmt(*out);
00364     bool these_results = true;
00365 
00366     result = true;
00367     TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(NOTRANS), true, *oss, result);
00368     if (!result) these_results = false;
00369 
00370     if(result) {
00371     
00372       *oss
00373         <<endl<< "Checking that the forward solve matches the forward operator to a residual tolerance:\n"
00374         <<endl<< "  v3 = inv(Op)*Op*v1"
00375         <<endl<< "               \\___/"
00376         <<endl<< "                 v2"
00377         <<endl<< ""
00378         <<endl<< "  v4 = Op*v3-v2"
00379         <<endl<< ""
00380         <<endl<< "  norm(v4)/norm(v2) <= forward_residual_solve_tol() + forward_residual_slack_error_tol()"
00381         <<endl;
00382 
00383       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00384 
00385         OSTab tab2(oss);
00386 
00387         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00388 
00389         tab.incrTab();
00390       
00391         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00392         Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(domain,l_num_rhs);
00393         Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*v1 );
00394         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00395       
00396         *oss <<endl<< "v2 = Op*v1 ...\n" ;
00397         Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(range,l_num_rhs);
00398         timer.start(true);
00399         Thyra::apply(op, NOTRANS, *v1, v2.ptr());
00400         timer.stop();
00401         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00402         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00403 
00404         *oss <<endl<< "v3 = inv(Op)*v2 ...\n" ;
00405         Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(domain,l_num_rhs);
00406         SolveCriteria<Scalar> solveCriteria(
00407           SolveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00408           ,forward_residual_solve_tol()
00409           );
00410         assign(&*v3, ST::zero());
00411         SolveStatus<Scalar> solveStatus;
00412         {
00413           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00414           timer.start(true);
00415           solveStatus = solve<Scalar>(op, NOTRANS, *v2, v3.ptr(),
00416             optInArg(solveCriteria));
00417           timer.stop();
00418           OSTab(oss).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00419         }
00420         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00421         *oss
00422           <<endl<< "solve status:\n";
00423         OSTab(oss).o() << solveStatus;
00424         result = solveStatus.solveStatus==SOLVE_STATUS_CONVERGED;
00425         if(!result) these_results = false;
00426         *oss
00427           <<endl<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00428           << passfail(result)<<endl;
00429       
00430         *oss <<endl<< "v4 = Op*v3 - v2 ...\n" ;
00431         Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(range,l_num_rhs);
00432         assign( &*v4, *v2 );
00433         timer.start(true);
00434         Thyra::apply(op, NOTRANS, *v3, v4.ptr(), Scalar(1.0), Scalar(-1.0));
00435         timer.stop();
00436         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00437         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00438 
00439         std::vector<ScalarMag> norms_v2(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v2(l_num_rhs);
00440         norms(*v2,&norms_v2[0]);
00441         norms(*v4,&norms_v4[0]);
00442         std::transform(
00443           norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.begin()
00444           ,std::divides<ScalarMag>()
00445           );
00446 
00447         result = testMaxErrors(
00448           l_num_rhs, "norm(v4)/norm(v2)", &norms_v4_norms_v2[0]
00449           ,"forward_residual_solve_tol()+forward_residual_slack_error_tol()", ScalarMag(forward_residual_solve_tol()+forward_residual_slack_error_tol())
00450           ,"forward_residual_solve_tol()_slack_warning_tol()", ScalarMag(forward_residual_solve_tol()+forward_residual_slack_warning_tol())
00451           ,&*oss
00452           );
00453         if(!result) these_results = false;
00454 
00455       }
00456     }
00457     else {
00458       *oss <<endl<< "Forward operator not supported, skipping check!\n";
00459     }
00460 
00461     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00462 
00463   }
00464   else {
00465     if(out.get()) *out <<endl<< "this->check_forward_residual()==false: Skipping the check of the forward solve with a tolerance on the residual!\n";
00466   }
00467   
00468   if( check_adjoint_default() ) {
00469 
00470     if(out.get()) *out <<endl<< "this->check_adjoint_default()==true: Checking the default adjoint solve ... ";
00471 
00472     std::ostringstream ossStore;
00473     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00474     if(out.get()) ossStore.copyfmt(*out);
00475     bool these_results = true;
00476 
00477     result = true;
00478     TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(CONJTRANS), true, *oss, result);
00479     if (!result) these_results = false;
00480 
00481     if(result) {
00482     
00483       *oss
00484         <<endl<< "Checking that the adjoint default solve matches the adjoint operator:\n"
00485         <<endl<< "  inv(Op')*Op'*v1 == v1"
00486         <<endl<< "           \\____/"
00487         <<endl<< "             v2"
00488         <<endl<< "  \\_____________/"
00489         <<endl<< "         v3"
00490         <<endl<< ""
00491         <<endl<< "  v4 = v3-v1"
00492         <<endl<< "  v5 = Op'*v3-v2"
00493         <<endl<< ""
00494         <<endl<< "  norm(v4)/norm(v1) <= adjoint_default_solution_error_error_tol()"
00495         <<endl<< "  norm(v5)/norm(v2) <= adjoint_default_residual_error_tol()"
00496         <<endl;
00497 
00498       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00499 
00500         OSTab tab2(oss);
00501 
00502         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00503         
00504         tab.incrTab();
00505       
00506         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00507         Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(range,l_num_rhs);
00508         Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*v1 );
00509         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00510       
00511         *oss <<endl<< "v2 = Op'*v1 ...\n" ;
00512         Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,l_num_rhs);
00513         timer.start(true);
00514         Thyra::apply(op, CONJTRANS, *v1, v2.ptr());
00515         timer.stop();
00516         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00517         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00518 
00519         *oss <<endl<< "v3 = inv(Op')*v2 ...\n" ;
00520         Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,l_num_rhs);
00521         assign(&*v3, ST::zero());
00522         SolveStatus<Scalar> solveStatus;
00523         {
00524           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00525           timer.start(true);
00526           solveStatus = solve<Scalar>(op, CONJTRANS, *v2, v3.ptr());
00527           timer.stop();
00528           OSTab(oss).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00529         }
00530         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00531         *oss
00532           <<endl<< "solve status:\n";
00533         OSTab(oss).o() << solveStatus;
00534 
00535         *oss <<endl<< "v4 = v3 - v1 ...\n" ;
00536         Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(range,l_num_rhs);
00537         V_VmV( &*v4, *v3, *v1 );
00538         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00539       
00540         *oss <<endl<< "v5 = Op'*v3 - v2 ...\n" ;
00541         Teuchos::RCP<MultiVectorBase<Scalar> > v5 = createMembers(domain,l_num_rhs);
00542         assign( &*v5, *v2 );
00543         timer.start(true);
00544         Thyra::apply(op, CONJTRANS, *v3, v5.ptr(), Scalar(1.0), Scalar(-1.0));
00545         timer.stop();
00546         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00547         if(dump_all()) *oss <<endl<< "v5 =\n" << describe(*v5,verbLevel);
00548 
00549         std::vector<ScalarMag> norms_v1(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v1(l_num_rhs);
00550         norms(*v1,&norms_v1[0]);
00551         norms(*v4,&norms_v4[0]);
00552         std::transform(
00553           norms_v4.begin(),norms_v4.end(),norms_v1.begin(),norms_v4_norms_v1.begin()
00554           ,std::divides<ScalarMag>()
00555           );
00556 
00557         result = testMaxErrors(
00558           l_num_rhs, "norm(v4)/norm(v1)", &norms_v4_norms_v1[0]
00559           ,"adjoint_default_solution_error_error_tol()", adjoint_default_solution_error_error_tol()
00560           ,"adjoint_default_solution_error_warning_tol()", adjoint_default_solution_error_warning_tol()
00561           ,&*oss
00562           );
00563         if(!result) these_results = false;
00564 
00565         std::vector<ScalarMag> norms_v2(l_num_rhs), norms_v5(l_num_rhs), norms_v5_norms_v2(l_num_rhs);
00566         norms(*v2,&norms_v2[0]);
00567         norms(*v5,&norms_v5[0]);
00568         std::transform(
00569           norms_v5.begin(),norms_v5.end(),norms_v2.begin(),norms_v5_norms_v2.begin()
00570           ,std::divides<ScalarMag>()
00571           );
00572 
00573         result = testMaxErrors(
00574           l_num_rhs, "norm(v5)/norm(v2)", &norms_v5_norms_v2[0]
00575           ,"adjoint_default_residual_error_tol()", adjoint_default_residual_error_tol()
00576           ,"adjoint_default_residual_warning_tol()", adjoint_default_residual_warning_tol()
00577           ,&*oss
00578           );
00579         if(!result) these_results = false;
00580 
00581       }
00582     }
00583     else {
00584       *oss <<endl<< "Adjoint operator not supported, skipping check!\n";
00585     }
00586 
00587     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00588 
00589   }
00590   else {
00591     if(out.get()) *out <<endl<< "this->check_adjoint_default()==false: Skipping the check of the adjoint solve with a default tolerance!\n";
00592   }
00593   
00594   if( check_adjoint_residual() ) {
00595 
00596     if(out.get()) *out <<endl<< "this->check_adjoint_residual()==true: Checking the adjoint solve with a tolerance on the residual ... ";
00597     
00598     std::ostringstream ossStore;
00599     Teuchos::RCP<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00600     if(out.get()) ossStore.copyfmt(*out);
00601     bool these_results = true;
00602 
00603     result = true;
00604     TEUCHOS_TEST_EQUALITY_CONST(op.solveSupports(CONJTRANS), true, *oss, result);
00605     if (!result) these_results = false;
00606 
00607     if(result) {
00608     
00609       *oss
00610         <<endl<< "Checking that the adjoint solve matches the adjoint operator to a residual tolerance:\n"
00611         <<endl<< "  v3 = inv(Op')*Op'*v1"
00612         <<endl<< "                \\____/"
00613         <<endl<< "                  v2"
00614         <<endl<< ""
00615         <<endl<< "  v4 = Op'*v3-v2"
00616         <<endl<< ""
00617         <<endl<< "  norm(v4)/norm(v2) <= adjoint_residual_solve_tol() + adjoint_residual_slack_error_tol()"
00618         <<endl;
00619 
00620       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00621 
00622         OSTab tab2(oss);
00623 
00624         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00625 
00626         tab.incrTab();
00627       
00628         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00629         Teuchos::RCP<MultiVectorBase<Scalar> > v1 = createMembers(range,l_num_rhs);
00630         Thyra::randomize( Scalar(-1.0), Scalar(+1.0), &*v1 );
00631         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00632       
00633         *oss <<endl<< "v2 = Op'*v1 ...\n" ;
00634         Teuchos::RCP<MultiVectorBase<Scalar> > v2 = createMembers(domain,l_num_rhs);
00635         timer.start(true);
00636         Thyra::apply(op, CONJTRANS, *v1, v2.ptr());
00637         timer.stop();
00638         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00639         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00640 
00641         *oss <<endl<< "v3 = inv(Op')*v2 ...\n" ;
00642         Teuchos::RCP<MultiVectorBase<Scalar> > v3 = createMembers(range,l_num_rhs);
00643         SolveCriteria<Scalar> solveCriteria(
00644           SolveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00645           ,adjoint_residual_solve_tol()
00646           );
00647         assign(&*v3, ST::zero());
00648         SolveStatus<Scalar> solveStatus;
00649         {
00650           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00651           timer.start(true);
00652           solveStatus = solve<Scalar>(op, CONJTRANS, *v2, v3.ptr(), optInArg(solveCriteria));
00653           timer.stop();
00654           OSTab(oss).o() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00655         }
00656         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00657         *oss
00658           <<endl<< "solve status:\n";
00659         OSTab(oss).o() << solveStatus;
00660         result = solveStatus.solveStatus==SOLVE_STATUS_CONVERGED;
00661         if(!result) these_results = false;
00662         *oss
00663           <<endl<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00664           << passfail(result)<<endl;
00665         if(solveStatus.achievedTol==SolveStatus<Scalar>::unknownTolerance())
00666           *oss <<endl<<"achievedTol==unknownTolerance(): Setting achievedTol = adjoint_residual_solve_tol() = "<<adjoint_residual_solve_tol()<<endl;
00667       
00668         *oss <<endl<< "v4 = Op'*v3 - v2 ...\n" ;
00669         Teuchos::RCP<MultiVectorBase<Scalar> > v4 = createMembers(domain,l_num_rhs);
00670         assign( &*v4, *v2 );
00671         timer.start(true);
00672         Thyra::apply(op, CONJTRANS, *v3, v4.ptr(), Scalar(1.0), Scalar(-1.0));
00673         timer.stop();
00674         OSTab(oss).o() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00675         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00676 
00677         std::vector<ScalarMag> norms_v2(l_num_rhs), norms_v4(l_num_rhs), norms_v4_norms_v2(l_num_rhs);
00678         norms(*v2,&norms_v2[0]);
00679         norms(*v4,&norms_v4[0]);
00680         std::transform(
00681           norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.begin()
00682           ,std::divides<ScalarMag>()
00683           );
00684 
00685         result = testMaxErrors(
00686           l_num_rhs, "norm(v4)/norm(v2)", &norms_v4_norms_v2[0],
00687           "adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol()",
00688           ScalarMag(adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol()),
00689           "adjoint_residual_solve_tol()_slack_warning_tol()",
00690           ScalarMag(adjoint_residual_solve_tol()+adjoint_residual_slack_warning_tol()),
00691           &*oss
00692           );
00693         if(!result) these_results = false;
00694 
00695       }
00696     }
00697     else {
00698       *oss <<endl<< "Adjoint operator not supported, skipping check!\n";
00699     }
00700 
00701     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).get());
00702 
00703   }
00704   else {
00705     if(out.get())
00706       *out <<endl<< "this->check_adjoint_residual()==false: Skipping the check of the adjoint solve with a tolerance on the residual!\n";
00707   }
00708   
00709   if(out.get()) {
00710     if(success)
00711       *out <<endl<<"Congratulations, this LinearOpWithSolveBase object seems to check out!\n";
00712     else
00713       *out <<endl<<"Oh no, at least one of the tests performed with this LinearOpWithSolveBase object failed (see above failures)!\n";
00714     *out <<endl<< "*** Leaving LinearOpWithSolveTester<"<<ST::name()<<">::check(...)\n";
00715   }
00716   
00717   return success;
00718   
00719 }
00720 
00721 
00722 // private static members
00723 
00724 
00725 // Define local macros used in the next few lines and then undefined
00726 
00727 #define LOWST_DEFINE_RAW_STATIC_MEMBER( DATA_TYPE, MEMBER_NAME, DEFAULT_VALUE ) \
00728 template<class Scalar> \
00729 const DATA_TYPE \
00730 LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
00731 
00732 #define LOWST_DEFINE_MTD_STATIC_MEMBER( DATA_TYPE, MEMBER_NAME, DEFAULT_VALUE ) \
00733 template<class Scalar> \
00734 const typename LinearOpWithSolveTester<Scalar>::DATA_TYPE \
00735 LinearOpWithSolveTester<Scalar>::MEMBER_NAME = DEFAULT_VALUE
00736 
00737 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_default_default_, true);
00738 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_forward_residual_default_, true);
00739 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_default_default_, true);
00740 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, check_adjoint_residual_default_, true);
00741 
00742 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, warning_tol_default_, 1e-6);
00743 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, error_tol_default_, 1e-5);
00744 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, solve_tol_default_, 1e-5);
00745 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_warning_tol_default_, 1e-6);
00746 LOWST_DEFINE_MTD_STATIC_MEMBER(ScalarMag, slack_error_tol_default_, 1e-5);
00747 
00748 LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_random_vectors_default_, 1);
00749 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, show_all_tests_default_, false);
00750 LOWST_DEFINE_RAW_STATIC_MEMBER(bool, dump_all_default_, false);
00751 LOWST_DEFINE_RAW_STATIC_MEMBER(int, num_rhs_default_, 1);
00752 
00753 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSolveTol_name_, "All Solve Tol");
00754 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackWarningTol_name_, "All Slack Warning Tol");
00755 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, AllSlackErrorTol_name_, "All Slack Error Tol");
00756 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, ShowAllTests_name_, "Show All Tests");
00757 LOWST_DEFINE_RAW_STATIC_MEMBER(std::string, DumpAll_name_, "Dump All");
00758 
00759 #undef LOWST_DEFINE_MTD_STATIC_MEMBER
00760 
00761 #undef LOWST_DEFINE_RAW_STATIC_MEMBER
00762 
00763 
00764 } // namespace Thyra
00765 
00766 
00767 #endif // THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines