Thyra_LinearOpWithSolveTester.hpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
00030 #define THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
00031 
00032 #include "Thyra_LinearOpWithSolveTesterDecl.hpp"
00033 #include "Thyra_LinearOpWithSolveBase.hpp"
00034 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00035 #include "Thyra_describeLinearOp.hpp"
00036 #include "Thyra_VectorStdOps.hpp"
00037 #include "Thyra_TestingTools.hpp"
00038 #include "Teuchos_Time.hpp"
00039 
00040 namespace Thyra {
00041 
00042 template<class RangeScalar, class DomainScalar>
00043 LinearOpWithSolveTester<RangeScalar,DomainScalar>::LinearOpWithSolveTester(
00044   const bool                 check_forward_default
00045   ,const RangeScalarMag      forward_default_residual_warning_tol
00046   ,const RangeScalarMag      forward_default_residual_error_tol
00047   ,const DomainScalarMag     forward_default_solution_error_warning_tol
00048   ,const DomainScalarMag     forward_default_solution_error_error_tol
00049   ,const bool                check_forward_residual
00050   ,const RangeScalarMag      forward_residual_solve_tol
00051   ,const RangeScalarMag      forward_residual_slack_warning_tol
00052   ,const RangeScalarMag      forward_residual_slack_error_tol
00053   ,const bool                check_forward_solution_error
00054   ,const RangeScalarMag      forward_solution_error_solve_tol
00055   ,const RangeScalarMag      forward_solution_error_slack_warning_tol
00056   ,const RangeScalarMag      forward_solution_error_slack_error_tol
00057   ,const bool                check_adjoint_default
00058   ,const DomainScalarMag     adjoint_default_residual_warning_tol
00059   ,const DomainScalarMag     adjoint_default_residual_error_tol
00060   ,const RangeScalarMag      adjoint_default_solution_error_warning_tol
00061   ,const RangeScalarMag      adjoint_default_solution_error_error_tol
00062   ,const bool                check_adjoint_residual
00063   ,const DomainScalarMag     adjoint_residual_solve_tol
00064   ,const DomainScalarMag     adjoint_residual_slack_warning_tol
00065   ,const DomainScalarMag     adjoint_residual_slack_error_tol
00066   ,const bool                check_adjoint_solution_error
00067   ,const DomainScalarMag     adjoint_solution_error_solve_tol
00068   ,const DomainScalarMag     adjoint_solution_error_slack_warning_tol
00069   ,const DomainScalarMag     adjoint_solution_error_slack_error_tol
00070   ,const int                 num_random_vectors
00071   ,const bool                show_all_tests
00072   ,const bool                dump_all
00073   ,const int                 num_rhs
00074   )
00075   :check_forward_default_(check_forward_default)
00076   ,forward_default_residual_warning_tol_(forward_default_residual_warning_tol)
00077   ,forward_default_residual_error_tol_(forward_default_residual_error_tol)
00078   ,forward_default_solution_error_warning_tol_(forward_default_solution_error_warning_tol)
00079   ,forward_default_solution_error_error_tol_(forward_default_solution_error_error_tol)
00080   ,check_forward_residual_(check_forward_residual)
00081   ,forward_residual_solve_tol_(forward_residual_solve_tol)
00082   ,forward_residual_slack_warning_tol_(forward_residual_slack_warning_tol)
00083   ,forward_residual_slack_error_tol_(forward_residual_slack_error_tol)
00084   ,check_forward_solution_error_(check_forward_solution_error)
00085   ,forward_solution_error_solve_tol_(forward_solution_error_solve_tol)
00086   ,forward_solution_error_slack_warning_tol_(forward_solution_error_slack_warning_tol)
00087   ,forward_solution_error_slack_error_tol_(forward_solution_error_slack_error_tol)
00088   ,check_adjoint_default_(check_adjoint_default)
00089   ,adjoint_default_residual_warning_tol_(adjoint_default_residual_warning_tol)
00090   ,adjoint_default_residual_error_tol_(adjoint_default_residual_error_tol)
00091   ,adjoint_default_solution_error_warning_tol_(adjoint_default_solution_error_warning_tol)
00092   ,adjoint_default_solution_error_error_tol_(adjoint_default_solution_error_error_tol)
00093   ,check_adjoint_residual_(check_adjoint_residual)
00094   ,adjoint_residual_solve_tol_(adjoint_residual_solve_tol)
00095   ,adjoint_residual_slack_warning_tol_(adjoint_residual_slack_warning_tol)
00096   ,adjoint_residual_slack_error_tol_(adjoint_residual_slack_error_tol)
00097   ,check_adjoint_solution_error_(check_adjoint_solution_error)
00098   ,adjoint_solution_error_solve_tol_(adjoint_solution_error_solve_tol)
00099   ,adjoint_solution_error_slack_warning_tol_(adjoint_solution_error_slack_warning_tol)
00100   ,adjoint_solution_error_slack_error_tol_(adjoint_solution_error_slack_error_tol)
00101   ,num_random_vectors_(num_random_vectors)
00102   ,show_all_tests_(show_all_tests)
00103   ,dump_all_(dump_all)
00104   ,num_rhs_(num_rhs)
00105 {}
00106 
00107 template<class RangeScalar, class DomainScalar>
00108 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::turn_off_all_tests()
00109 {
00110   check_forward_default_         = false;
00111   check_forward_residual_        = false;
00112   check_forward_solution_error_  = false;
00113   check_adjoint_default_         = false;
00114   check_adjoint_residual_        = false;
00115   check_adjoint_solution_error_  = false;
00116 }
00117 
00118 template<class RangeScalar, class DomainScalar>
00119 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::set_all_solve_tol( const ScalarMag solve_tol )
00120 {
00121   forward_residual_solve_tol_ = solve_tol;
00122   forward_residual_solve_tol_ = solve_tol;
00123   forward_solution_error_solve_tol_ = solve_tol;
00124   adjoint_residual_solve_tol_ = solve_tol;
00125   adjoint_solution_error_solve_tol_ = solve_tol;
00126 }
00127 
00128 template<class RangeScalar, class DomainScalar>
00129 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::set_all_slack_warning_tol( const ScalarMag slack_warning_tol )
00130 {
00131   forward_default_residual_warning_tol_ = slack_warning_tol;
00132   forward_default_solution_error_warning_tol_ = slack_warning_tol;
00133   forward_residual_slack_warning_tol_ = slack_warning_tol;
00134   forward_solution_error_slack_warning_tol_ = slack_warning_tol;
00135   adjoint_default_residual_warning_tol_ = slack_warning_tol;
00136   adjoint_default_solution_error_warning_tol_ = slack_warning_tol;
00137   adjoint_residual_slack_warning_tol_ = slack_warning_tol;
00138   adjoint_solution_error_slack_warning_tol_ = slack_warning_tol;
00139 }
00140 
00141 template<class RangeScalar, class DomainScalar>
00142 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::set_all_slack_error_tol( const ScalarMag slack_error_tol )
00143 {
00144   forward_default_residual_error_tol_ = slack_error_tol;
00145   forward_default_solution_error_error_tol_ = slack_error_tol;
00146   forward_residual_slack_error_tol_ = slack_error_tol;
00147   forward_solution_error_slack_error_tol_ = slack_error_tol;
00148   adjoint_default_residual_error_tol_ = slack_error_tol;
00149   adjoint_default_solution_error_error_tol_ = slack_error_tol;
00150   adjoint_residual_slack_error_tol_ = slack_error_tol;
00151   adjoint_solution_error_slack_error_tol_ = slack_error_tol;
00152 }
00153   
00154 template<class RangeScalar, class DomainScalar>
00155 bool LinearOpWithSolveTester<RangeScalar,DomainScalar>::check(
00156   const LinearOpWithSolveBase<RangeScalar,DomainScalar>   &op
00157   ,Teuchos::FancyOStream                                  *out_arg
00158   ) const
00159 {
00160 
00161   using std::endl;
00162   using Teuchos::arrayArg;
00163   using Teuchos::FancyOStream;
00164   using Teuchos::OSTab;
00165   typedef Teuchos::VerboseObjectTempState<LinearOpWithSolveBase<RangeScalar,DomainScalar> > VOTS;
00166   typedef Teuchos::ScalarTraits<Scalar>        ST;
00167   typedef Teuchos::ScalarTraits<RangeScalar>   DST;
00168   typedef Teuchos::ScalarTraits<DomainScalar>  RST;
00169   bool success = true, result;
00170   const int num_rhs = this->num_rhs();
00171   Teuchos::RefCountPtr<FancyOStream> out = Teuchos::rcp(out_arg,false);
00172   const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00173   Teuchos::Time timer("");
00174 
00175   OSTab tab(out,1,"THYRA");
00176 
00177   if(out.get()) {
00178     *out <<endl<< "*** Entering LinearOpWithSolveTester<"<<ST::name()<<">::check(op,...) ...\n";
00179     if(show_all_tests()) {
00180       *out <<endl<< "describe forward op:\n" << Teuchos::describe(op,verbLevel);
00181       if(op.applyTransposeSupports(CONJ_ELE) && verbLevel==Teuchos::VERB_EXTREME) {
00182         *out <<endl<< "describe adjoint op:\n";
00183         describeLinearOp<RangeScalar,DomainScalar>(
00184           *adjoint(Teuchos::rcp_implicit_cast< const LinearOpBase<RangeScalar,DomainScalar> >(Teuchos::rcp(&op,false)))
00185           ,*out,verbLevel
00186           );
00187       }
00188     }
00189     else {
00190       *out <<endl<< "describe op: " << op.description() << endl;
00191     }
00192   }
00193   
00194   Teuchos::RefCountPtr<const VectorSpaceBase<RangeScalar> >   range  = op.range();
00195   Teuchos::RefCountPtr<const VectorSpaceBase<DomainScalar> >  domain = op.domain();
00196   
00197   if( check_forward_default() ) {
00198 
00199     if(out.get()) *out <<endl<< "this->check_forward_default()==true: Checking the default forward solve ... ";
00200 
00201     std::ostringstream ossStore;
00202     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00203     if(out.get()) ossStore.copyfmt(*out);
00204     bool these_results = true;
00205 
00206     *oss <<endl<< "op.solveSupportsConj(NONCONJ_ELE) == true ? ";
00207     result = op.solveSupportsConj(NONCONJ_ELE);
00208     if(!result) these_results = false;
00209     *oss << passfail(result) << endl;
00210 
00211     if(result) {
00212     
00213       *oss
00214         <<endl<< "Checking that the forward default solve matches the forward operator:\n"
00215         <<endl<< "  inv(Op)*Op*v1 == v1"
00216         <<endl<< "          \\___/"
00217         <<endl<< "           v2"
00218         <<endl<< "  \\___________/"
00219         <<endl<< "         v3"
00220         <<endl<< ""
00221         <<endl<< "  v4 = v3-v1"
00222         <<endl<< "  v5 = Op*v3-v2"
00223         <<endl<< ""
00224         <<endl<< "  norm(v4)/norm(v1) <= forward_default_solution_error_error_tol()"
00225         <<endl<< "  norm(v5)/norm(v2) <= forward_default_residual_error_tol()"
00226         <<endl;
00227 
00228       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00229 
00230         OSTab tab(oss);
00231 
00232         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00233 
00234         tab.incrTab();
00235       
00236         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00237         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00238         Thyra::randomize( DomainScalar(-1.0), DomainScalar(+1.0), &*v1 );
00239         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00240       
00241         *oss <<endl<< "v2 = Op*v1 ...\n" ;
00242         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00243         timer.start(true);
00244         op.apply(NONCONJ_ELE,*v1,&*v2);
00245         timer.stop();
00246         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00247         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00248 
00249         *oss <<endl<< "v3 = inv(Op)*v2 ...\n" ;
00250         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v3 = createMembers(domain,num_rhs);
00251         assign(&*v3,DST::zero());
00252         SolveStatus<Scalar> solveStatus;
00253         if(1){
00254           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00255           timer.start(true);
00256           solveStatus = solve(op,NONCONJ_ELE,*v2,&*v3,static_cast<const SolveCriteria<Scalar>*>(0));
00257           timer.stop();
00258           *OSTab(oss).getOStream() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00259         }
00260         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00261         *oss
00262           <<endl<< "solve status:\n";
00263         *OSTab(oss).getOStream() << solveStatus;
00264 
00265         *oss <<endl<< "v4 = v3 - v1 ...\n" ;
00266         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v4 = createMembers(domain,num_rhs);
00267         V_VmV( &*v4, *v3, *v1 );
00268         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00269       
00270         *oss <<endl<< "v5 = Op*v3 - v2 ...\n" ;
00271         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v5 = createMembers(range,num_rhs);
00272         assign( &*v5, *v2 );
00273         timer.start(true);
00274         op.apply(NONCONJ_ELE,*v3,&*v5,Scalar(1.0),Scalar(-1.0));
00275         timer.stop();
00276         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00277         if(dump_all()) *oss <<endl<< "v5 =\n" << describe(*v5,verbLevel);
00278 
00279         std::vector<DomainScalarMag> norms_v1(num_rhs), norms_v4(num_rhs), norms_v4_norms_v1(num_rhs);
00280         norms(*v1,&norms_v1[0]);
00281         norms(*v4,&norms_v4[0]);
00282         std::transform(
00283           norms_v4.begin(),norms_v4.end(),norms_v1.begin(),norms_v4_norms_v1.begin()
00284           ,std::divides<DomainScalarMag>()
00285           );
00286 
00287         result = testMaxErrors(
00288           num_rhs, "norm(v4)/norm(v1)", &norms_v4_norms_v1[0]
00289           ,"forward_default_solution_error_error_tol()", forward_default_solution_error_error_tol()
00290           ,"forward_default_solution_error_warning_tol()", forward_default_solution_error_warning_tol()
00291           ,&*oss
00292           );
00293         if(!result) these_results = false;
00294 
00295         std::vector<RangeScalarMag> norms_v2(num_rhs), norms_v5(num_rhs), norms_v5_norms_v2(num_rhs);
00296         norms(*v2,&norms_v2[0]);
00297         norms(*v5,&norms_v5[0]);
00298         std::transform(
00299           norms_v5.begin(),norms_v5.end(),norms_v2.begin(),norms_v5_norms_v2.begin()
00300           ,std::divides<RangeScalarMag>()
00301           );
00302 
00303         result = testMaxErrors(
00304           num_rhs, "norm(v5)/norm(v2)", &norms_v5_norms_v2[0]
00305           ,"forward_default_residual_error_tol()", forward_default_residual_error_tol()
00306           ,"forward_default_residual_warning_tol()", forward_default_residual_warning_tol()
00307           ,&*oss
00308           );
00309         if(!result) these_results = false;
00310         
00311       }
00312     }
00313     else {
00314       *oss <<endl<< "Forward operator not supported, skipping check!\n";
00315     }
00316 
00317     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00318 
00319   }
00320   else {
00321     if(out.get()) *out <<endl<< "this->check_forward_default()==false: Skipping the check of the default forward solve!\n";
00322   }
00323   
00324   if( check_forward_residual() ) {
00325 
00326     if(out.get()) *out <<endl<< "this->check_forward_residual()==true: Checking the forward solve with a tolerance on the residual ... ";
00327 
00328     std::ostringstream ossStore;
00329     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00330     if(out.get()) ossStore.copyfmt(*out);
00331     bool these_results = true;
00332 
00333     *oss <<endl<< "op.solveSupportsConj(NONCONJ_ELE) == true ? ";
00334     result = op.solveSupportsConj(NONCONJ_ELE);
00335     if(!result) these_results = false;
00336     *oss << passfail(result) << endl;
00337 
00338     if(result) {
00339     
00340       *oss
00341         <<endl<< "Checking that the forward solve matches the forward operator to a residual tolerance:\n"
00342         <<endl<< "  v3 = inv(Op)*Op*v1"
00343         <<endl<< "               \\___/"
00344         <<endl<< "                 v2"
00345         <<endl<< ""
00346         <<endl<< "  v4 = Op*v3-v2"
00347         <<endl<< ""
00348         <<endl<< "  norm(v4)/norm(v2) <= forward_residual_solve_tol() + forward_residual_slack_error_tol()"
00349         <<endl;
00350 
00351       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00352 
00353         OSTab tab(oss);
00354 
00355         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00356 
00357         tab.incrTab();
00358       
00359         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00360         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v1 = createMembers(domain,num_rhs);
00361         Thyra::randomize( DomainScalar(-1.0), DomainScalar(+1.0), &*v1 );
00362         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00363       
00364         *oss <<endl<< "v2 = Op*v1 ...\n" ;
00365         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v2 = createMembers(range,num_rhs);
00366         timer.start(true);
00367         op.apply(NONCONJ_ELE,*v1,&*v2);
00368         timer.stop();
00369         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00370         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00371 
00372         *oss <<endl<< "v3 = inv(Op)*v2 ...\n" ;
00373         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v3 = createMembers(domain,num_rhs);
00374         SolveCriteria<Scalar> solveCriteria(
00375           SolveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00376           ,forward_residual_solve_tol()
00377           );
00378         assign(&*v3,DST::zero());
00379         SolveStatus<Scalar> solveStatus;
00380         if(1){
00381           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00382           timer.start(true);
00383           solveStatus = solve<RangeScalar,DomainScalar>(op,NONCONJ_ELE,*v2,&*v3,&solveCriteria);
00384           timer.stop();
00385           *OSTab(oss).getOStream() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00386         }
00387         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00388         *oss
00389           <<endl<< "solve status:\n";
00390         *OSTab(oss).getOStream() << solveStatus;
00391         result = solveStatus.solveStatus==SOLVE_STATUS_CONVERGED;
00392         if(!result) these_results = false;
00393         *oss
00394           <<endl<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00395           << passfail(result)<<endl;
00396       
00397         *oss <<endl<< "v4 = Op*v3 - v2 ...\n" ;
00398         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v4 = createMembers(range,num_rhs);
00399         assign( &*v4, *v2 );
00400         timer.start(true);
00401         op.apply(NONCONJ_ELE,*v3,&*v4,Scalar(1.0),Scalar(-1.0));
00402         timer.stop();
00403         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00404         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00405 
00406         std::vector<DomainScalarMag> norms_v2(num_rhs), norms_v4(num_rhs), norms_v4_norms_v2(num_rhs);
00407         norms(*v2,&norms_v2[0]);
00408         norms(*v4,&norms_v4[0]);
00409         std::transform(
00410           norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.begin()
00411           ,std::divides<DomainScalarMag>()
00412           );
00413 
00414         result = testMaxErrors(
00415           num_rhs, "norm(v4)/norm(v2)", &norms_v4_norms_v2[0]
00416           ,"forward_residual_solve_tol()+forward_residual_slack_error_tol()", DomainScalarMag(forward_residual_solve_tol()+forward_residual_slack_error_tol())
00417           ,"forward_residual_solve_tol()_slack_warning_tol()", DomainScalarMag(forward_residual_solve_tol()+forward_residual_slack_warning_tol())
00418           ,&*oss
00419           );
00420         if(!result) these_results = false;
00421 
00422       }
00423     }
00424     else {
00425       *oss <<endl<< "Forward operator not supported, skipping check!\n";
00426     }
00427 
00428     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00429 
00430   }
00431   else {
00432     if(out.get()) *out <<endl<< "this->check_forward_residual()==false: Skipping the check of the forward solve with a tolerance on the residual!\n";
00433   }
00434   
00435   if( check_adjoint_default() ) {
00436 
00437     if(out.get()) *out <<endl<< "this->check_adjoint_default()==true: Checking the default adjoint solve ... ";
00438 
00439     std::ostringstream ossStore;
00440     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00441     if(out.get()) ossStore.copyfmt(*out);
00442     bool these_results = true;
00443 
00444     *oss <<endl<< "op.solveTransposeSupportsConj(CONJ_ELE) == true ? ";
00445     result = op.solveTransposeSupportsConj(CONJ_ELE);
00446     if(!result) these_results = false;
00447     *oss << passfail(result) << endl;
00448 
00449     if(result) {
00450     
00451       *oss
00452         <<endl<< "Checking that the adjoint default solve matches the adjoint operator:\n"
00453         <<endl<< "  inv(Op')*Op'*v1 == v1"
00454         <<endl<< "           \\____/"
00455         <<endl<< "             v2"
00456         <<endl<< "  \\_____________/"
00457         <<endl<< "         v3"
00458         <<endl<< ""
00459         <<endl<< "  v4 = v3-v1"
00460         <<endl<< "  v5 = Op'*v3-v2"
00461         <<endl<< ""
00462         <<endl<< "  norm(v4)/norm(v1) <= adjoint_default_solution_error_error_tol()"
00463         <<endl<< "  norm(v5)/norm(v2) <= adjoint_default_residual_error_tol()"
00464         <<endl;
00465 
00466       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00467 
00468         OSTab tab(oss);
00469 
00470         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00471         
00472         tab.incrTab();
00473       
00474         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00475         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v1 = createMembers(range,num_rhs);
00476         Thyra::randomize( RangeScalar(-1.0), RangeScalar(+1.0), &*v1 );
00477         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00478       
00479         *oss <<endl<< "v2 = Op'*v1 ...\n" ;
00480         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v2 = createMembers(domain,num_rhs);
00481         timer.start(true);
00482         op.applyTranspose(CONJ_ELE,*v1,&*v2);
00483         timer.stop();
00484         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00485         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00486 
00487         *oss <<endl<< "v3 = inv(Op')*v2 ...\n" ;
00488         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00489         assign(&*v3,DST::zero());
00490         SolveStatus<Scalar> solveStatus;
00491         if(1){
00492           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00493           timer.start(true);
00494           solveStatus = solveTranspose(op,CONJ_ELE,*v2,&*v3,static_cast<const SolveCriteria<Scalar>*>(0));
00495           timer.stop();
00496           *OSTab(oss).getOStream() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00497         }
00498         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00499         *oss
00500           <<endl<< "solve status:\n";
00501         *OSTab(oss).getOStream() << solveStatus;
00502 
00503         *oss <<endl<< "v4 = v3 - v1 ...\n" ;
00504         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v4 = createMembers(range,num_rhs);
00505         V_VmV( &*v4, *v3, *v1 );
00506         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00507       
00508         *oss <<endl<< "v5 = Op'*v3 - v2 ...\n" ;
00509         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v5 = createMembers(domain,num_rhs);
00510         assign( &*v5, *v2 );
00511         timer.start(true);
00512         op.applyTranspose(CONJ_ELE,*v3,&*v5,Scalar(1.0),Scalar(-1.0));
00513         timer.stop();
00514         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00515         if(dump_all()) *oss <<endl<< "v5 =\n" << describe(*v5,verbLevel);
00516 
00517         std::vector<RangeScalarMag> norms_v1(num_rhs), norms_v4(num_rhs), norms_v4_norms_v1(num_rhs);
00518         norms(*v1,&norms_v1[0]);
00519         norms(*v4,&norms_v4[0]);
00520         std::transform(
00521           norms_v4.begin(),norms_v4.end(),norms_v1.begin(),norms_v4_norms_v1.begin()
00522           ,std::divides<RangeScalarMag>()
00523           );
00524 
00525         result = testMaxErrors(
00526           num_rhs, "norm(v4)/norm(v1)", &norms_v4_norms_v1[0]
00527           ,"adjoint_default_solution_error_error_tol()", adjoint_default_solution_error_error_tol()
00528           ,"adjoint_default_solution_error_warning_tol()", adjoint_default_solution_error_warning_tol()
00529           ,&*oss
00530           );
00531         if(!result) these_results = false;
00532 
00533         std::vector<DomainScalarMag> norms_v2(num_rhs), norms_v5(num_rhs), norms_v5_norms_v2(num_rhs);
00534         norms(*v2,&norms_v2[0]);
00535         norms(*v5,&norms_v5[0]);
00536         std::transform(
00537           norms_v5.begin(),norms_v5.end(),norms_v2.begin(),norms_v5_norms_v2.begin()
00538           ,std::divides<DomainScalarMag>()
00539           );
00540 
00541         result = testMaxErrors(
00542           num_rhs, "norm(v5)/norm(v2)", &norms_v5_norms_v2[0]
00543           ,"adjoint_default_residual_error_tol()", adjoint_default_residual_error_tol()
00544           ,"adjoint_default_residual_warning_tol()", adjoint_default_residual_warning_tol()
00545           ,&*oss
00546           );
00547         if(!result) these_results = false;
00548 
00549       }
00550     }
00551     else {
00552       *oss <<endl<< "Adjoint operator not supported, skipping check!\n";
00553     }
00554 
00555     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00556 
00557   }
00558   else {
00559     if(out.get()) *out <<endl<< "this->check_adjoint_default()==false: Skipping the check of the adjoint solve with a default tolerance!\n";
00560   }
00561   
00562   if( check_adjoint_residual() ) {
00563 
00564     if(out.get()) *out <<endl<< "this->check_adjoint_residual()==true: Checking the adjoint solve with a tolerance on the residual ... ";
00565     
00566     std::ostringstream ossStore;
00567     Teuchos::RefCountPtr<FancyOStream> oss = Teuchos::rcp(new FancyOStream(Teuchos::rcp(&ossStore,false)));
00568     if(out.get()) ossStore.copyfmt(*out);
00569     bool these_results = true;
00570 
00571     *oss <<endl<< "op.solveTransposeSupportsConj(CONJ_ELE) == true ? ";
00572     result = op.solveTransposeSupportsConj(CONJ_ELE);
00573     if(!result) these_results = false;
00574     *oss << passfail(result) << endl;
00575 
00576     if(result) {
00577     
00578       *oss
00579         <<endl<< "Checking that the adjoint solve matches the adjoint operator to a residual tolerance:\n"
00580         <<endl<< "  v3 = inv(Op')*Op'*v1"
00581         <<endl<< "                \\____/"
00582         <<endl<< "                  v2"
00583         <<endl<< ""
00584         <<endl<< "  v4 = Op'*v3-v2"
00585         <<endl<< ""
00586         <<endl<< "  norm(v4)/norm(v2) <= adjoint_residual_solve_tol() + adjoint_residual_slack_error_tol()"
00587         <<endl;
00588 
00589       for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00590 
00591         OSTab tab(oss);
00592 
00593         *oss <<endl<< "Random vector tests = " << rand_vec_i << endl;
00594 
00595         tab.incrTab();
00596       
00597         *oss <<endl<< "v1 = randomize(-1,+1); ...\n" ;
00598         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v1 = createMembers(range,num_rhs);
00599         Thyra::randomize( RangeScalar(-1.0), RangeScalar(+1.0), &*v1 );
00600         if(dump_all()) *oss <<endl<< "v1 =\n" << describe(*v1,verbLevel);
00601       
00602         *oss <<endl<< "v2 = Op'*v1 ...\n" ;
00603         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v2 = createMembers(domain,num_rhs);
00604         timer.start(true);
00605         op.applyTranspose(CONJ_ELE,*v1,&*v2);
00606         timer.stop();
00607         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00608         if(dump_all()) *oss <<endl<< "v2 =\n" << describe(*v2,verbLevel);
00609 
00610         *oss <<endl<< "v3 = inv(Op')*v2 ...\n" ;
00611         Teuchos::RefCountPtr<MultiVectorBase<RangeScalar> > v3 = createMembers(range,num_rhs);
00612         SolveCriteria<Scalar> solveCriteria(
00613           SolveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00614           ,adjoint_residual_solve_tol()
00615           );
00616         assign(&*v3,RST::zero());
00617         SolveStatus<Scalar> solveStatus;
00618         if(1){
00619           VOTS lowsTempState(Teuchos::rcp(&op,false),oss,verbLevel);
00620           timer.start(true);
00621           solveStatus = solveTranspose(op,CONJ_ELE,*v2,&*v3,&solveCriteria);
00622           timer.stop();
00623           *OSTab(oss).getOStream() <<"\n=> Solve time = "<<timer.totalElapsedTime()<<" sec\n";
00624         }
00625         if(dump_all()) *oss <<endl<< "v3 =\n" << describe(*v3,verbLevel);
00626         *oss
00627           <<endl<< "solve status:\n";
00628         *OSTab(oss).getOStream() << solveStatus;
00629         result = solveStatus.solveStatus==SOLVE_STATUS_CONVERGED;
00630         if(!result) these_results = false;
00631         *oss
00632           <<endl<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00633           << passfail(result)<<endl;
00634         if(solveStatus.achievedTol==SolveStatus<RangeScalar>::unknownTolerance())
00635           *oss <<endl<<"achievedTol==unknownTolerance(): Setting achievedTol = adjoint_residual_solve_tol() = "<<adjoint_residual_solve_tol()<<endl;
00636       
00637         *oss <<endl<< "v4 = Op'*v3 - v2 ...\n" ;
00638         Teuchos::RefCountPtr<MultiVectorBase<DomainScalar> > v4 = createMembers(domain,num_rhs);
00639         assign( &*v4, *v2 );
00640         timer.start(true);
00641         op.applyTranspose(CONJ_ELE,*v3,&*v4,Scalar(1.0),Scalar(-1.0));
00642         timer.stop();
00643         *OSTab(oss).getOStream() <<"\n=> Apply time = "<<timer.totalElapsedTime()<<" sec\n";
00644         if(dump_all()) *oss <<endl<< "v4 =\n" << describe(*v4,verbLevel);
00645 
00646         std::vector<RangeScalarMag> norms_v2(num_rhs), norms_v4(num_rhs), norms_v4_norms_v2(num_rhs);
00647         norms(*v2,&norms_v2[0]);
00648         norms(*v4,&norms_v4[0]);
00649         std::transform(
00650           norms_v4.begin(),norms_v4.end(),norms_v2.begin(),norms_v4_norms_v2.begin()
00651           ,std::divides<RangeScalarMag>()
00652           );
00653 
00654         result = testMaxErrors(
00655           num_rhs, "norm(v4)/norm(v2)", &norms_v4_norms_v2[0]
00656           ,"adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol()", RangeScalarMag(adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol())
00657           ,"adjoint_residual_solve_tol()_slack_warning_tol()", RangeScalarMag(adjoint_residual_solve_tol()+adjoint_residual_slack_warning_tol())
00658           ,&*oss
00659           );
00660         if(!result) these_results = false;
00661 
00662       }
00663     }
00664     else {
00665       *oss <<endl<< "Adjoint operator not supported, skipping check!\n";
00666     }
00667 
00668     printTestResults(these_results,ossStore.str(),show_all_tests(),&success,OSTab(out).getOStream().get());
00669 
00670   }
00671   else {
00672     if(out.get()) *out <<endl<< "this->check_adjoint_residual()==false: Skipping the check of the adjoint solve with a tolerance on the residual!\n";
00673   }
00674   
00675   if(out.get()) {
00676     if(success)
00677       *out <<endl<<"Congratulations, this LinearOpWithSolveBase object seems to check out!\n";
00678     else
00679       *out <<endl<<"Oh no, at least one of the tests performed with this LinearOpWithSolveBase object failed (see above failures)!\n";
00680     *out <<endl<< "*** Leaving LinearOpWithSolveTester<"<<ST::name()<<">::check(...)\n";
00681   }
00682   
00683   return success;
00684   
00685 }
00686 
00687 } // namespace Thyra
00688 
00689 #endif // THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP

Generated on Thu Sep 18 12:32:45 2008 for Thyra Operator Solve Support by doxygen 1.3.9.1