00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
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 }
00688
00689 #endif // THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
