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