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_ScaledAdjointLinearOp.hpp"
00035 #include "Thyra_describeLinearOp.hpp"
00036 #include "Thyra_VectorStdOps.hpp"
00037 #include "Thyra_TestingTools.hpp"
00038
00039 namespace Thyra {
00040
00041 template<class RangeScalar, class DomainScalar>
00042 LinearOpWithSolveTester<RangeScalar,DomainScalar>::LinearOpWithSolveTester(
00043 const bool check_forward_default
00044 ,const RangeScalarMag forward_default_residual_warning_tol
00045 ,const RangeScalarMag forward_default_residual_error_tol
00046 ,const DomainScalarMag forward_default_solution_error_warning_tol
00047 ,const DomainScalarMag forward_default_solution_error_error_tol
00048 ,const bool check_forward_residual
00049 ,const RangeScalarMag forward_residual_solve_tol
00050 ,const RangeScalarMag forward_residual_slack_warning_tol
00051 ,const RangeScalarMag forward_residual_slack_error_tol
00052 ,const bool check_forward_solution_error
00053 ,const RangeScalarMag forward_solution_error_solve_tol
00054 ,const RangeScalarMag forward_solution_error_slack_warning_tol
00055 ,const RangeScalarMag forward_solution_error_slack_error_tol
00056 ,const bool check_adjoint_default
00057 ,const DomainScalarMag adjoint_default_residual_warning_tol
00058 ,const DomainScalarMag adjoint_default_residual_error_tol
00059 ,const RangeScalarMag adjoint_default_solution_error_warning_tol
00060 ,const RangeScalarMag adjoint_default_solution_error_error_tol
00061 ,const bool check_adjoint_residual
00062 ,const DomainScalarMag adjoint_residual_solve_tol
00063 ,const DomainScalarMag adjoint_residual_slack_warning_tol
00064 ,const DomainScalarMag adjoint_residual_slack_error_tol
00065 ,const bool check_adjoint_solution_error
00066 ,const DomainScalarMag adjoint_solution_error_solve_tol
00067 ,const DomainScalarMag adjoint_solution_error_slack_warning_tol
00068 ,const DomainScalarMag adjoint_solution_error_slack_error_tol
00069 ,const int num_random_vectors
00070 ,const bool show_all_tests
00071 ,const bool dump_all
00072 )
00073 :check_forward_default_(check_forward_default)
00074 ,forward_default_residual_warning_tol_(forward_default_residual_warning_tol)
00075 ,forward_default_residual_error_tol_(forward_default_residual_error_tol)
00076 ,forward_default_solution_error_warning_tol_(forward_default_solution_error_warning_tol)
00077 ,forward_default_solution_error_error_tol_(forward_default_solution_error_error_tol)
00078 ,check_forward_residual_(check_forward_residual)
00079 ,forward_residual_solve_tol_(forward_residual_solve_tol)
00080 ,forward_residual_slack_warning_tol_(forward_residual_slack_warning_tol)
00081 ,forward_residual_slack_error_tol_(forward_residual_slack_error_tol)
00082 ,check_forward_solution_error_(check_forward_solution_error)
00083 ,forward_solution_error_solve_tol_(forward_solution_error_solve_tol)
00084 ,forward_solution_error_slack_warning_tol_(forward_solution_error_slack_warning_tol)
00085 ,forward_solution_error_slack_error_tol_(forward_solution_error_slack_error_tol)
00086 ,check_adjoint_default_(check_adjoint_default)
00087 ,adjoint_default_residual_warning_tol_(adjoint_default_residual_warning_tol)
00088 ,adjoint_default_residual_error_tol_(adjoint_default_residual_error_tol)
00089 ,adjoint_default_solution_error_warning_tol_(adjoint_default_solution_error_warning_tol)
00090 ,adjoint_default_solution_error_error_tol_(adjoint_default_solution_error_error_tol)
00091 ,check_adjoint_residual_(check_adjoint_residual)
00092 ,adjoint_residual_solve_tol_(adjoint_residual_solve_tol)
00093 ,adjoint_residual_slack_warning_tol_(adjoint_residual_slack_warning_tol)
00094 ,adjoint_residual_slack_error_tol_(adjoint_residual_slack_error_tol)
00095 ,check_adjoint_solution_error_(check_adjoint_solution_error)
00096 ,adjoint_solution_error_solve_tol_(adjoint_solution_error_solve_tol)
00097 ,adjoint_solution_error_slack_warning_tol_(adjoint_solution_error_slack_warning_tol)
00098 ,adjoint_solution_error_slack_error_tol_(adjoint_solution_error_slack_error_tol)
00099 ,num_random_vectors_(num_random_vectors)
00100 ,show_all_tests_(show_all_tests)
00101 ,dump_all_(dump_all)
00102 {}
00103
00104 template<class RangeScalar, class DomainScalar>
00105 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::turn_off_all_tests()
00106 {
00107 check_forward_default_ = false;
00108 check_forward_residual_ = false;
00109 check_forward_solution_error_ = false;
00110 check_adjoint_default_ = false;
00111 check_adjoint_residual_ = false;
00112 check_adjoint_solution_error_ = false;
00113 }
00114
00115 template<class RangeScalar, class DomainScalar>
00116 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::set_all_solve_tol( const ScalarMag solve_tol )
00117 {
00118 forward_residual_solve_tol_ = solve_tol;
00119 forward_residual_solve_tol_ = solve_tol;
00120 forward_solution_error_solve_tol_ = solve_tol;
00121 adjoint_residual_solve_tol_ = solve_tol;
00122 adjoint_solution_error_solve_tol_ = solve_tol;
00123 }
00124
00125 template<class RangeScalar, class DomainScalar>
00126 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::set_all_slack_warning_tol( const ScalarMag slack_warning_tol )
00127 {
00128 forward_default_residual_warning_tol_ = slack_warning_tol;
00129 forward_default_solution_error_warning_tol_ = slack_warning_tol;
00130 forward_residual_slack_warning_tol_ = slack_warning_tol;
00131 forward_solution_error_slack_warning_tol_ = slack_warning_tol;
00132 adjoint_default_residual_warning_tol_ = slack_warning_tol;
00133 adjoint_default_solution_error_warning_tol_ = slack_warning_tol;
00134 adjoint_residual_slack_warning_tol_ = slack_warning_tol;
00135 adjoint_solution_error_slack_warning_tol_ = slack_warning_tol;
00136 }
00137
00138 template<class RangeScalar, class DomainScalar>
00139 void LinearOpWithSolveTester<RangeScalar,DomainScalar>::set_all_slack_error_tol( const ScalarMag slack_error_tol )
00140 {
00141 forward_default_residual_error_tol_ = slack_error_tol;
00142 forward_default_solution_error_error_tol_ = slack_error_tol;
00143 forward_residual_slack_error_tol_ = slack_error_tol;
00144 forward_solution_error_slack_error_tol_ = slack_error_tol;
00145 adjoint_default_residual_error_tol_ = slack_error_tol;
00146 adjoint_default_solution_error_error_tol_ = slack_error_tol;
00147 adjoint_residual_slack_error_tol_ = slack_error_tol;
00148 adjoint_solution_error_slack_error_tol_ = slack_error_tol;
00149 }
00150
00151 template<class RangeScalar, class DomainScalar>
00152 bool LinearOpWithSolveTester<RangeScalar,DomainScalar>::check(
00153 const LinearOpWithSolveBase<RangeScalar,DomainScalar> &op
00154 ,std::ostream *out
00155 ,const std::string &leadingIndent
00156 ,const std::string &indentSpacer
00157 ) const
00158 {
00159
00160 using std::endl;
00161 using Teuchos::arrayArg;
00162 typedef Teuchos::ScalarTraits<Scalar> ST;
00163 typedef Teuchos::ScalarTraits<RangeScalar> DST;
00164 typedef Teuchos::ScalarTraits<DomainScalar> RST;
00165 bool success = true, result;
00166 const std::string &li = leadingIndent, &is = indentSpacer;
00167 const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);
00168
00169 if(out) {
00170 *out <<endl<<li<< "*** Entering LinearOpWithSolveTester<"<<ST::name()<<">::check(op,...) ...\n";
00171 if(show_all_tests()) {
00172 *out <<endl<<li<< "describe forward op:\n" << Teuchos::describe(op,verbLevel,li,is);
00173 if(op.applyTransposeSupports(CONJ_ELE) && verbLevel==Teuchos::VERB_EXTREME) {
00174 *out <<endl<<li<< "describe adjoint op:\n";
00175 describeLinearOp(
00176 *adjoint(Teuchos::rcp_implicit_cast< const LinearOpBase<RangeScalar,DomainScalar> >(Teuchos::rcp(&op,false)))
00177 ,*out,verbLevel,li,is
00178 );
00179 }
00180 }
00181 else {
00182 *out <<endl<<li<< "describe op: " << op.description() << endl;
00183 }
00184 }
00185
00186 Teuchos::RefCountPtr<const VectorSpaceBase<RangeScalar> > range = op.range();
00187 Teuchos::RefCountPtr<const VectorSpaceBase<DomainScalar> > domain = op.domain();
00188
00189 if( check_forward_default() ) {
00190
00191 if(out) *out <<endl<<li<< "this->check_forward_default()==true: Checking the default forward solve ... ";
00192
00193 std::ostringstream oss;
00194 if(out) oss.copyfmt(*out);
00195 bool these_results = true;
00196
00197 oss <<endl<<li<<is<< "op.solveSupportsConj(NONCONJ_ELE) == true ? ";
00198 result = op.solveSupportsConj(NONCONJ_ELE);
00199 if(!result) these_results = false;
00200 oss << passfail(result) << endl;
00201
00202 if(result) {
00203
00204 oss
00205 <<endl<<li<<is<< "Checking that the forward default solve matches the forward operator:\n"
00206 <<endl<<li<<is<< " inv(Op)*Op*v1 == v1"
00207 <<endl<<li<<is<< " \\___/"
00208 <<endl<<li<<is<< " v2"
00209 <<endl<<li<<is<< " \\___________/"
00210 <<endl<<li<<is<< " v3"
00211 <<endl<<li<<is<< ""
00212 <<endl<<li<<is<< " v4 = v3-v1"
00213 <<endl<<li<<is<< " v5 = Op*v3-v2"
00214 <<endl<<li<<is<< ""
00215 <<endl<<li<<is<< " norm(v4)/norm(v1) <= forward_default_solution_error_error_tol()"
00216 <<endl<<li<<is<< " norm(v5)/norm(v2) <= forward_default_residual_error_tol()"
00217 <<endl;
00218
00219 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00220
00221 oss <<endl<<li<<is<< "Random vector tests = " << rand_vec_i << endl;
00222
00223 oss <<endl<<li<<is<< "v1 = randomize(-1,+1); ...\n" ;
00224 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v1 = createMember(domain);
00225 Thyra::randomize( DomainScalar(-1.0), DomainScalar(+1.0), &*v1 );
00226 if(dump_all()) oss <<endl<<li<<is<< "v1 =\n" << describe(*v1,verbLevel,li,is);
00227
00228 oss <<endl<<li<<is<< "v2 = Op*v1 ...\n" ;
00229 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v2 = createMember(range);
00230 op.apply(NONCONJ_ELE,*v1,&*v2);
00231 if(dump_all()) oss <<endl<<li<<is<< "v2 =\n" << describe(*v2,verbLevel,li,is);
00232
00233 oss <<endl<<li<<is<< "v3 = inv(Op)*v2 ...\n" ;
00234 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v3 = createMember(domain);
00235 assign(&*v3,DST::zero());
00236 SolveStatus<Scalar> solveStatus = solve(op,NONCONJ_ELE,*v2,&*v3,static_cast<const SolveCriteria<Scalar>*>(0));
00237 if(dump_all()) oss <<endl<<li<<is<< "v3 =\n" << describe(*v3,verbLevel,li,is);
00238 oss
00239 <<endl<<li<<is<< "solve status:"
00240 <<endl<<li<<is<<is<< "solveStatus = " << toString(solveStatus.solveStatus)
00241 <<endl<<li<<is<<is<< "achievedTol = " << SolveStatus<RangeScalar>::achievedTolToString(solveStatus.achievedTol)
00242 <<endl<<li<<is<<is<< "iterations = " << solveStatus.iterations
00243 <<endl<<li<<is<<is<< "message = \"" << solveStatus.message << "\""
00244 <<endl;
00245
00246 oss <<endl<<li<<is<< "v4 = v3 - v1 ...\n" ;
00247 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v4 = createMember(domain);
00248 V_VmV( &*v4, *v3, *v1 );
00249 if(dump_all()) oss <<endl<<li<<is<< "v4 =\n" << describe(*v4,verbLevel,li,is);
00250
00251 oss <<endl<<li<<is<< "v5 = Op*v3 - v2 ...\n" ;
00252 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v5 = createMember(range);
00253 assign( &*v5, *v2 );
00254 op.apply(NONCONJ_ELE,*v3,&*v5,Scalar(1.0),Scalar(-1.0));
00255 if(dump_all()) oss <<endl<<li<<is<< "v5 =\n" << describe(*v5,verbLevel,li,is);
00256
00257 const DomainScalarMag
00258 norm_v1 = norm(*v1),
00259 norm_v4 = norm(*v4),
00260 norm_v4_norm_v1 = norm_v4/norm_v1;
00261
00262 result = testMaxErr(
00263 "norm(v4)/norm(v1)", norm_v4_norm_v1
00264 ,"forward_default_solution_error_error_tol()", forward_default_solution_error_error_tol()
00265 ,"forward_default_solution_error_warning_tol()", forward_default_solution_error_warning_tol()
00266 ,&oss,li+is
00267 );
00268 if(!result) these_results = false;
00269
00270 const RangeScalarMag
00271 norm_v2 = norm(*v2),
00272 norm_v5 = norm(*v5),
00273 norm_v5_norm_v2 = norm_v5/norm_v2;
00274
00275 result = testMaxErr(
00276 "norm(v5)/norm(v2)", norm_v5_norm_v2
00277 ,"forward_default_residual_error_tol()", forward_default_residual_error_tol()
00278 ,"forward_default_residual_warning_tol()", forward_default_residual_warning_tol()
00279 ,&oss,li+is
00280 );
00281 if(!result) these_results = false;
00282
00283 }
00284 }
00285 else {
00286 oss <<endl<<li<<is<< "Forward operator not supported, skipping check!\n";
00287 }
00288
00289 printTestResults(these_results,oss.str(),show_all_tests(),&success,out);
00290
00291 }
00292 else {
00293 if(out) *out <<endl<<li<< "this->check_forward_default()==false: Skipping the check of the default forward solve!\n";
00294 }
00295
00296 if( check_forward_residual() ) {
00297
00298 if(out) *out <<endl<<li<< "this->check_forward_residual()==true: Checking the forward solve with a tolerance on the residual ... ";
00299
00300 std::ostringstream oss;
00301 if(out) oss.copyfmt(*out);
00302 bool these_results = true;
00303
00304 oss <<endl<<li<<is<< "op.solveSupportsConj(NONCONJ_ELE) == true ? ";
00305 result = op.solveSupportsConj(NONCONJ_ELE);
00306 if(!result) these_results = false;
00307 oss << passfail(result) << endl;
00308
00309 if(result) {
00310
00311 oss
00312 <<endl<<li<<is<< "Checking that the forward solve matches the forward operator to a residual tolerance:\n"
00313 <<endl<<li<<is<< " v3 = inv(Op)*Op*v1"
00314 <<endl<<li<<is<< " \\___/"
00315 <<endl<<li<<is<< " v2"
00316 <<endl<<li<<is<< ""
00317 <<endl<<li<<is<< " v4 = Op*v3-v2"
00318 <<endl<<li<<is<< ""
00319 <<endl<<li<<is<< " norm(v4)/norm(v2) <= forward_residual_solve_tol() + forward_residual_slack_error_tol()"
00320 <<endl;
00321
00322 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00323
00324 oss <<endl<<li<<is<< "Random vector tests = " << rand_vec_i << endl;
00325
00326 oss <<endl<<li<<is<< "v1 = randomize(-1,+1); ...\n" ;
00327 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v1 = createMember(domain);
00328 Thyra::randomize( DomainScalar(-1.0), DomainScalar(+1.0), &*v1 );
00329 if(dump_all()) oss <<endl<<li<<is<< "v1 =\n" << describe(*v1,verbLevel,li,is);
00330
00331 oss <<endl<<li<<is<< "v2 = Op*v1 ...\n" ;
00332 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v2 = createMember(range);
00333 op.apply(NONCONJ_ELE,*v1,&*v2);
00334 if(dump_all()) oss <<endl<<li<<is<< "v2 =\n" << describe(*v2,verbLevel,li,is);
00335
00336 oss <<endl<<li<<is<< "v3 = inv(Op)*v2 ...\n" ;
00337 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v3 = createMember(domain);
00338 SolveCriteria<Scalar> solveCriteria(SOLVE_TOL_REL_RESIDUAL_NORM,forward_residual_solve_tol());
00339 assign(&*v3,DST::zero());
00340 SolveStatus<Scalar> solveStatus = solve<RangeScalar,DomainScalar>(op,NONCONJ_ELE,*v2,&*v3,&solveCriteria);
00341 if(dump_all()) oss <<endl<<li<<is<< "v3 =\n" << describe(*v3,verbLevel,li,is);
00342 oss
00343 <<endl<<li<<is<< "solve status:"
00344 <<endl<<li<<is<<is<< "solveStatus = " << toString(solveStatus.solveStatus)
00345 <<endl<<li<<is<<is<< "achievedTol = " << SolveStatus<RangeScalar>::achievedTolToString(solveStatus.achievedTol)
00346 <<endl<<li<<is<<is<< "iterations = " << solveStatus.iterations
00347 <<endl<<li<<is<<is<< "message = \"" << solveStatus.message << "\""
00348 <<endl;
00349 oss
00350 <<endl<<li<<is<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00351 << passfail(solveStatus.solveStatus==SOLVE_STATUS_CONVERGED)<<endl;
00352
00353 oss <<endl<<li<<is<< "v4 = Op*v3 - v2 ...\n" ;
00354 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v4 = createMember(range);
00355 assign( &*v4, *v2 );
00356 op.apply(NONCONJ_ELE,*v3,&*v4,Scalar(1.0),Scalar(-1.0));
00357 if(dump_all()) oss <<endl<<li<<is<< "v4 =\n" << describe(*v4,verbLevel,li,is);
00358
00359 const RangeScalarMag
00360 norm_v2 = norm(*v2),
00361 norm_v4 = norm(*v4),
00362 norm_v4_norm_v2 = norm_v4/norm_v2;
00363
00364 result = testMaxErr(
00365 "norm(v4)/norm(v2)", norm_v4_norm_v2
00366 ,"forward_residual_solve_tol()+forward_residual_slack_error_tol()", RangeScalarMag(forward_residual_solve_tol()+forward_residual_slack_error_tol())
00367 ,"forward_residual_solve_tol()_slack_warning_tol()", RangeScalarMag(forward_residual_solve_tol()+forward_residual_slack_warning_tol())
00368 ,&oss,li+is
00369 );
00370 if(!result) these_results = false;
00371
00372 }
00373 }
00374 else {
00375 oss <<endl<<li<<is<< "Forward operator not supported, skipping check!\n";
00376 }
00377
00378 printTestResults(these_results,oss.str(),show_all_tests(),&success,out);
00379
00380 }
00381 else {
00382 if(out) *out <<endl<<li<< "this->check_forward_residual()==false: Skipping the check of the forward solve with a tolerance on the residual!\n";
00383 }
00384
00385 if( check_forward_solution_error() ) {
00386
00387 if(out) *out <<endl<<li<< "this->check_forward_solution_error()==true: Checking the forward solve with a tolerance on the solution error ... ";
00388
00389 std::ostringstream oss;
00390 if(out) oss.copyfmt(*out);
00391 bool these_results = true;
00392
00393 oss <<endl<<li<<is<< "op.solveSupportsConj(NONCONJ_ELE) == true ? ";
00394 result = op.solveSupportsConj(NONCONJ_ELE);
00395 if(!result) these_results = false;
00396 oss << passfail(result) << endl;
00397
00398 if(result) {
00399
00400 oss
00401 <<endl<<li<<is<< "Checking that the forward solve matches the forward operator to a solution error tolerance:\n"
00402 <<endl<<li<<is<< " v3 = inv(Op)*Op*v1"
00403 <<endl<<li<<is<< " \\___/"
00404 <<endl<<li<<is<< " v2"
00405 <<endl<<li<<is<< ""
00406 <<endl<<li<<is<< " v4 = v3-v1"
00407 <<endl<<li<<is<< ""
00408 <<endl<<li<<is<< " norm(v4)/norm(v1) <= forward_solution_error_solve_tol() + forward_solution_error_slack_error_tol()"
00409 <<endl;
00410
00411 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00412
00413 oss <<endl<<li<<is<< "Random vector tests = " << rand_vec_i << endl;
00414
00415 oss <<endl<<li<<is<< "v1 = randomize(-1,+1); ...\n" ;
00416 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v1 = createMember(domain);
00417 Thyra::randomize( DomainScalar(-1.0), DomainScalar(+1.0), &*v1 );
00418 if(dump_all()) oss <<endl<<li<<is<< "v1 =\n" << describe(*v1,verbLevel,li,is);
00419
00420 oss <<endl<<li<<is<< "v2 = Op*v1 ...\n" ;
00421 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v2 = createMember(range);
00422 op.apply(NONCONJ_ELE,*v1,&*v2);
00423 if(dump_all()) oss <<endl<<li<<is<< "v2 =\n" << describe(*v2,verbLevel,li,is);
00424
00425 oss <<endl<<li<<is<< "v3 = inv(Op)*v2 ...\n" ;
00426 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v3 = createMember(domain);
00427 SolveCriteria<Scalar> solveCriteria(SOLVE_TOL_REL_SOLUTION_ERR_NORM,forward_solution_error_solve_tol());
00428 assign(&*v3,DST::zero());
00429 SolveStatus<Scalar> solveStatus = solve(op,NONCONJ_ELE,*v2,&*v3,&solveCriteria);
00430 if(dump_all()) oss <<endl<<li<<is<< "v3 =\n" << describe(*v3,verbLevel,li,is);
00431 oss
00432 <<endl<<li<<is<< "solve status:"
00433 <<endl<<li<<is<<is<< "solveStatus = " << toString(solveStatus.solveStatus)
00434 <<endl<<li<<is<<is<< "achievedTol = " << SolveStatus<RangeScalar>::achievedTolToString(solveStatus.achievedTol)
00435 <<endl<<li<<is<<is<< "iterations = " << solveStatus.iterations
00436 <<endl<<li<<is<<is<< "message = \"" << solveStatus.message << "\""
00437 <<endl;
00438 oss
00439 <<endl<<li<<is<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00440 << passfail(solveStatus.solveStatus==SOLVE_STATUS_CONVERGED)<<endl;
00441
00442 oss <<endl<<li<<is<< "v4 = v3 - v1 ...\n" ;
00443 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v4 = createMember(domain);
00444 V_VmV( &*v4, *v3, *v1 );
00445 if(dump_all()) oss <<endl<<li<<is<< "v4 =\n" << describe(*v4,verbLevel,li,is);
00446
00447 const RangeScalarMag
00448 norm_v1 = norm(*v1),
00449 norm_v4 = norm(*v4),
00450 norm_v4_norm_v1 = norm_v4/norm_v1;
00451
00452 result = testMaxErr(
00453 "norm(v4)/norm(v1)", norm_v4_norm_v1
00454 ,"forward_solution_error_solve_tol()+forward_solution_error_slack_error_tol()", RangeScalarMag(forward_solution_error_solve_tol()+forward_solution_error_slack_error_tol())
00455 ,"forward_solution_error_solve_tol()+forward_solution_error_slack_warning_tol()", RangeScalarMag(forward_solution_error_solve_tol()+forward_solution_error_slack_warning_tol())
00456 ,&oss,li+is
00457 );
00458 if(!result) these_results = false;
00459
00460 }
00461 }
00462 else {
00463 oss <<endl<<li<<is<< "Forward operator not supported, skipping check!\n";
00464 }
00465
00466 printTestResults(these_results,oss.str(),show_all_tests(),&success,out);
00467
00468 }
00469 else {
00470 if(out) *out <<endl<<li<< "this->check_forward_solution_error()==false: Skipping the check of the forward solve with a tolerance on the solution_error!\n";
00471 }
00472
00473 if( check_adjoint_default() ) {
00474
00475 if(out) *out <<endl<<li<< "this->check_adjoint_default()==true: Checking the default adjoint solve ... ";
00476
00477 std::ostringstream oss;
00478 if(out) oss.copyfmt(*out);
00479 bool these_results = true;
00480
00481 oss <<endl<<li<<is<< "op.solveTransposeSupportsConj(CONJ_ELE) == true ? ";
00482 result = op.solveTransposeSupportsConj(CONJ_ELE);
00483 if(!result) these_results = false;
00484 oss << passfail(result) << endl;
00485
00486 if(result) {
00487
00488 oss
00489 <<endl<<li<<is<< "Checking that the adjoint default solve matches the adjoint operator:\n"
00490 <<endl<<li<<is<< " inv(Op')*Op'*v1 == v1"
00491 <<endl<<li<<is<< " \\____/"
00492 <<endl<<li<<is<< " v2"
00493 <<endl<<li<<is<< " \\_____________/"
00494 <<endl<<li<<is<< " v3"
00495 <<endl<<li<<is<< ""
00496 <<endl<<li<<is<< " v4 = v3-v1"
00497 <<endl<<li<<is<< " v5 = Op'*v3-v2"
00498 <<endl<<li<<is<< ""
00499 <<endl<<li<<is<< " norm(v4)/norm(v1) <= adjoint_default_solution_error_error_tol()"
00500 <<endl<<li<<is<< " norm(v5)/norm(v2) <= adjoint_default_residual_error_tol()"
00501 <<endl;
00502
00503 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00504
00505 oss <<endl<<li<<is<< "Random vector tests = " << rand_vec_i << endl;
00506
00507 oss <<endl<<li<<is<< "v1 = randomize(-1,+1); ...\n" ;
00508 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v1 = createMember(range);
00509 Thyra::randomize( RangeScalar(-1.0), RangeScalar(+1.0), &*v1 );
00510 if(dump_all()) oss <<endl<<li<<is<< "v1 =\n" << describe(*v1,verbLevel,li,is);
00511
00512 oss <<endl<<li<<is<< "v2 = Op*v1 ...\n" ;
00513 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v2 = createMember(domain);
00514 op.apply(NONCONJ_ELE,*v1,&*v2);
00515 if(dump_all()) oss <<endl<<li<<is<< "v2 =\n" << describe(*v2,verbLevel,li,is);
00516
00517 oss <<endl<<li<<is<< "v3 = inv(Op)*v2 ...\n" ;
00518 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v3 = createMember(range);
00519 assign(&*v3,DST::zero());
00520 SolveStatus<Scalar> solveStatus = solve(op,NONCONJ_ELE,*v2,&*v3,static_cast<const SolveCriteria<Scalar>*>(0));
00521 if(dump_all()) oss <<endl<<li<<is<< "v3 =\n" << describe(*v3,verbLevel,li,is);
00522 oss
00523 <<endl<<li<<is<< "solve status:"
00524 <<endl<<li<<is<<is<< "solveStatus = " << toString(solveStatus.solveStatus)
00525 <<endl<<li<<is<<is<< "achievedTol = " << SolveStatus<DomainScalar>::achievedTolToString(solveStatus.achievedTol)
00526 <<endl<<li<<is<<is<< "iterations = " << solveStatus.iterations
00527 <<endl<<li<<is<<is<< "message = \"" << solveStatus.message << "\""
00528 <<endl;
00529
00530 oss <<endl<<li<<is<< "v4 = v3 - v1 ...\n" ;
00531 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v4 = createMember(range);
00532 V_VmV( &*v4, *v3, *v1 );
00533 if(dump_all()) oss <<endl<<li<<is<< "v4 =\n" << describe(*v4,verbLevel,li,is);
00534
00535 oss <<endl<<li<<is<< "v5 = Op*v3 - v2 ...\n" ;
00536 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v5 = createMember(domain);
00537 assign( &*v5, *v2 );
00538 op.apply(NONCONJ_ELE,*v3,&*v5,Scalar(1.0),Scalar(-1.0));
00539 if(dump_all()) oss <<endl<<li<<is<< "v5 =\n" << describe(*v5,verbLevel,li,is);
00540
00541 const RangeScalarMag
00542 norm_v1 = norm(*v1),
00543 norm_v4 = norm(*v4),
00544 norm_v4_norm_v1 = norm_v4/norm_v1;
00545
00546 result = testMaxErr(
00547 "norm(v4)/norm(v1)", norm_v4_norm_v1
00548 ,"adjoint_default_solution_error_error_tol()", adjoint_default_solution_error_error_tol()
00549 ,"adjoint_default_solution_error_warning_tol()", adjoint_default_solution_error_warning_tol()
00550 ,&oss,li+is
00551 );
00552 if(!result) these_results = false;
00553
00554 const DomainScalarMag
00555 norm_v2 = norm(*v2),
00556 norm_v5 = norm(*v5),
00557 norm_v5_norm_v2 = norm_v5/norm_v2;
00558
00559 result = testMaxErr(
00560 "norm(v5)/norm(v2)", norm_v5_norm_v2
00561 ,"adjoint_default_residual_error_tol()", adjoint_default_residual_error_tol()
00562 ,"adjoint_default_residual_warning_tol()", adjoint_default_residual_warning_tol()
00563 ,&oss,li+is
00564 );
00565 if(!result) these_results = false;
00566
00567 }
00568 }
00569 else {
00570 oss <<endl<<li<<is<< "Adjoint operator not supported, skipping check!\n";
00571 }
00572
00573 printTestResults(these_results,oss.str(),show_all_tests(),&success,out);
00574
00575 }
00576 else {
00577 if(out) *out <<endl<<li<< "this->check_adjoint_default()==false: Skipping the check of the adjoint solve with a default tolerance!\n";
00578 }
00579
00580 if( check_adjoint_residual() ) {
00581
00582 if(out) *out <<endl<<li<< "this->check_adjoint_residual()==true: Checking the adjoint solve with a tolerance on the residual ... ";
00583
00584 std::ostringstream oss;
00585 if(out) oss.copyfmt(*out);
00586 bool these_results = true;
00587
00588 oss <<endl<<li<<is<< "op.solveTransposeSupportsConj(CONJ_ELE) == true ? ";
00589 result = op.solveTransposeSupportsConj(CONJ_ELE);
00590 if(!result) these_results = false;
00591 oss << passfail(result) << endl;
00592
00593 if(result) {
00594
00595 oss
00596 <<endl<<li<<is<< "Checking that the adjoint solve matches the adjoint operator to a residual tolerance:\n"
00597 <<endl<<li<<is<< " v3 = inv(Op')*Op'*v1"
00598 <<endl<<li<<is<< " \\____/"
00599 <<endl<<li<<is<< " v2"
00600 <<endl<<li<<is<< ""
00601 <<endl<<li<<is<< " v4 = Op'*v3-v2"
00602 <<endl<<li<<is<< ""
00603 <<endl<<li<<is<< " norm(v4)/norm(v2) <= adjoint_residual_solve_tol() + adjoint_residual_slack_error_tol()"
00604 <<endl;
00605
00606 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00607
00608 oss <<endl<<li<<is<< "Random vector tests = " << rand_vec_i << endl;
00609
00610 oss <<endl<<li<<is<< "v1 = randomize(-1,+1); ...\n" ;
00611 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v1 = createMember(range);
00612 Thyra::randomize( RangeScalar(-1.0), RangeScalar(+1.0), &*v1 );
00613 if(dump_all()) oss <<endl<<li<<is<< "v1 =\n" << describe(*v1,verbLevel,li,is);
00614
00615 oss <<endl<<li<<is<< "v2 = Op'*v1 ...\n" ;
00616 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v2 = createMember(domain);
00617 op.applyTranspose(CONJ_ELE,*v1,&*v2);
00618 if(dump_all()) oss <<endl<<li<<is<< "v2 =\n" << describe(*v2,verbLevel,li,is);
00619
00620 oss <<endl<<li<<is<< "v3 = inv(Op')*v2 ...\n" ;
00621 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v3 = createMember(range);
00622 SolveCriteria<Scalar> solveCriteria(SOLVE_TOL_REL_RESIDUAL_NORM,adjoint_residual_solve_tol());
00623 assign(&*v3,RST::zero());
00624 SolveStatus<Scalar> solveStatus = solveTranspose(op,CONJ_ELE,*v2,&*v3,&solveCriteria);
00625 if(dump_all()) oss <<endl<<li<<is<< "v3 =\n" << describe(*v3,verbLevel,li,is);
00626 oss
00627 <<endl<<li<<is<< "solve status:"
00628 <<endl<<li<<is<<is<< "solveStatus = " << toString(solveStatus.solveStatus)
00629 <<endl<<li<<is<<is<< "achievedTol = " << SolveStatus<RangeScalar>::achievedTolToString(solveStatus.achievedTol)
00630 <<endl<<li<<is<<is<< "iterations = " << solveStatus.iterations
00631 <<endl<<li<<is<<is<< "message = \"" << solveStatus.message << "\""
00632 <<endl;
00633 oss
00634 <<endl<<li<<is<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00635 << passfail(solveStatus.solveStatus==SOLVE_STATUS_CONVERGED)<<endl;
00636 if(solveStatus.achievedTol==SolveStatus<RangeScalar>::unknownTolerance())
00637 oss <<endl<<li<<is<<"achievedTol==unknownTolerance(): Setting achievedTol = adjoint_residual_solve_tol() = "<<adjoint_residual_solve_tol()<<endl;
00638
00639 oss <<endl<<li<<is<< "v4 = Op'*v3 - v2 ...\n" ;
00640 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v4 = createMember(domain);
00641 assign( &*v4, *v2 );
00642 op.applyTranspose(CONJ_ELE,*v3,&*v4,Scalar(1.0),Scalar(-1.0));
00643 if(dump_all()) oss <<endl<<li<<is<< "v4 =\n" << describe(*v4,verbLevel,li,is);
00644
00645 const DomainScalarMag
00646 norm_v2 = norm(*v2),
00647 norm_v4 = norm(*v4),
00648 norm_v4_norm_v2 = norm_v4/norm_v2;
00649
00650 result = testMaxErr(
00651 "norm(v4)/norm(v2)", norm_v4_norm_v2
00652 ,"adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol()", DomainScalarMag(adjoint_residual_solve_tol()+adjoint_residual_slack_error_tol())
00653 ,"adjoint_residual_solve_tol()+adjoint_residual_slack_warning_tol()", DomainScalarMag(adjoint_residual_solve_tol()+adjoint_residual_slack_warning_tol())
00654 ,&oss,li+is
00655 );
00656 if(!result) these_results = false;
00657
00658 }
00659 }
00660 else {
00661 oss <<endl<<li<<is<< "Adjoint operator not supported, skipping check!\n";
00662 }
00663
00664 printTestResults(these_results,oss.str(),show_all_tests(),&success,out);
00665
00666 }
00667 else {
00668 if(out) *out <<endl<<li<< "this->check_adjoint_residual()==false: Skipping the check of the adjoint solve with a tolerance on the residual!\n";
00669 }
00670
00671 if( check_adjoint_solution_error() ) {
00672
00673 if(out) *out <<endl<<li<< "this->check_adjoint_solution_error()==true: Checking the adjoint solve with a tolerance on the solution error ... ";
00674
00675 std::ostringstream oss;
00676 if(out) oss.copyfmt(*out);
00677 bool these_results = true;
00678
00679 oss <<endl<<li<<is<< "op.solveTransposeSupportsConj(CONJ_ELE) == true ? ";
00680 result = op.solveTransposeSupportsConj(CONJ_ELE);
00681 if(!result) these_results = false;
00682 oss << passfail(result) << endl;
00683
00684 if(result) {
00685
00686 oss
00687 <<endl<<li<<is<< "Checking that the adjoint solve matches the adjoint operator to a solution error tolerance:\n"
00688 <<endl<<li<<is<< " v3 = inv(Op')*Op'*v1"
00689 <<endl<<li<<is<< " \\____/"
00690 <<endl<<li<<is<< " v2"
00691 <<endl<<li<<is<< ""
00692 <<endl<<li<<is<< " v4 = v3-v1"
00693 <<endl<<li<<is<< ""
00694 <<endl<<li<<is<< " norm(v4)/norm(v1) <= adjoint_solution_error_solve_tol() + adjoint_solution_error_slack_error_tol()"
00695 <<endl;
00696
00697 for( int rand_vec_i = 1; rand_vec_i <= num_random_vectors(); ++rand_vec_i ) {
00698
00699 oss <<endl<<li<<is<< "Random vector tests = " << rand_vec_i << endl;
00700
00701 oss <<endl<<li<<is<< "v1 = randomize(-1,+1); ...\n" ;
00702 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v1 = createMember(range);
00703 Thyra::randomize( RangeScalar(-1.0), RangeScalar(+1.0), &*v1 );
00704 if(dump_all()) oss <<endl<<li<<is<< "v1 =\n" << describe(*v1,verbLevel,li,is);
00705
00706 oss <<endl<<li<<is<< "v2 = Op*v1 ...\n" ;
00707 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v2 = createMember(domain);
00708 op.applyTranspose(CONJ_ELE,*v1,&*v2);
00709 if(dump_all()) oss <<endl<<li<<is<< "v2 =\n" << describe(*v2,verbLevel,li,is);
00710
00711 oss <<endl<<li<<is<< "v3 = inv(Op)*v2 ...\n" ;
00712 Teuchos::RefCountPtr<VectorBase<RangeScalar> > v3 = createMember(range);
00713 SolveCriteria<Scalar> solveCriteria(SOLVE_TOL_REL_SOLUTION_ERR_NORM,adjoint_solution_error_solve_tol());
00714 assign(&*v3,RST::zero());
00715 SolveStatus<Scalar> solveStatus = solveTranspose(op,CONJ_ELE,*v2,&*v3,&solveCriteria);
00716 if(dump_all()) oss <<endl<<li<<is<< "v3 =\n" << describe(*v3,verbLevel,li,is);
00717 oss
00718 <<endl<<li<<is<< "solve status:"
00719 <<endl<<li<<is<<is<< "solveStatus = " << toString(solveStatus.solveStatus)
00720 <<endl<<li<<is<<is<< "achievedTol = " << SolveStatus<DomainScalar>::achievedTolToString(solveStatus.achievedTol)
00721 <<endl<<li<<is<<is<< "iterations = " << solveStatus.iterations
00722 <<endl<<li<<is<<is<< "message = \"" << solveStatus.message << "\""
00723 <<endl;
00724 oss
00725 <<endl<<li<<is<< "check: solveStatus = " << toString(solveStatus.solveStatus) << " == SOLVE_STATUS_CONVERGED : "
00726 << passfail(solveStatus.solveStatus==SOLVE_STATUS_CONVERGED)<<endl;
00727 if(solveStatus.achievedTol==SolveStatus<DomainScalar>::unknownTolerance())
00728 oss <<endl<<li<<is<<"achievedTol==unknownTolerance(): Setting achievedTol = adjoint_solution_error_solve_tol() = "<<adjoint_solution_error_solve_tol()<<endl;
00729
00730 oss <<endl<<li<<is<< "v4 = v3 - v1 ...\n" ;
00731 Teuchos::RefCountPtr<VectorBase<DomainScalar> > v4 = createMember(range);
00732 V_VmV( &*v4, *v3, *v1 );
00733 if(dump_all()) oss <<endl<<li<<is<< "v4 =\n" << describe(*v4,verbLevel,li,is);
00734
00735 const DomainScalarMag
00736 norm_v1 = norm(*v1),
00737 norm_v4 = norm(*v4),
00738 norm_v4_norm_v1 = norm_v4/norm_v1;
00739
00740 result = testMaxErr(
00741 "norm(v4)/norm(v1)", norm_v4_norm_v1
00742 ,"adjoint_solution_error_solve_tol()+adjoint_solution_error_slack_error_tol()", DomainScalarMag(adjoint_solution_error_solve_tol()+adjoint_solution_error_slack_error_tol())
00743 ,"adjoint_solution_error_solve_tol()+adjoint_solution_error_slack_warning_tol()", DomainScalarMag(adjoint_solution_error_solve_tol()+adjoint_solution_error_slack_warning_tol())
00744 ,&oss,li+is
00745 );
00746 if(!result) these_results = false;
00747
00748 }
00749 }
00750 else {
00751 oss <<endl<<li<<is<< "Adjoint operator not supported, skipping check!\n";
00752 }
00753
00754 printTestResults(these_results,oss.str(),show_all_tests(),&success,out);
00755
00756 }
00757 else {
00758 if(out) *out <<endl<<li<< "this->check_adjoint_solution_error()==false: Skipping the check of the adjoint solve with a tolerance on the solution_error!\n";
00759 }
00760
00761 if(out) {
00762 if(success)
00763 *out <<endl<<li<<"Congratulations, this LinearOpWithSolveBase object seems to check out!\n";
00764 else
00765 *out <<endl<<li<<"Oh no, at least one of the tests performed with this LinearOpWithSolveBase object failed (see above failures)!\n";
00766 *out <<endl<<li<< "*** Leaving LinearOpWithSolveTester<"<<ST::name()<<">::check(...)\n";
00767 }
00768
00769 return success;
00770
00771 }
00772
00773 }
00774
00775 #endif // THYRA_LINEAR_OP_WITH_SOLVE_TESTER_HPP
