DiagonalTransientInversionMain.cpp

Go to the documentation of this file.
00001 
00002 #include "EpetraExt_DiagonalTransientModel.hpp"
00003 #include "MoochoPack_MoochoThyraSolver.hpp"
00004 #include "Rythmos_BackwardEulerStepper.hpp"
00005 #include "Rythmos_ImplicitBDFStepper.hpp"
00006 #include "Rythmos_ForwardSensitivityStepper.hpp"
00007 #include "Rythmos_TimeStepNonlinearSolver.hpp"
00008 #include "Rythmos_DefaultIntegrator.hpp"
00009 #include "Rythmos_SimpleIntegrationControlStrategy.hpp"
00010 #include "Rythmos_StepperAsModelEvaluator.hpp"
00011 #include "Rythmos_ForwardSensitivityIntegratorAsModelEvaluator.hpp"
00012 #include "Thyra_EpetraModelEvaluator.hpp"
00013 #include "Thyra_DefaultSpmdMultiVectorFileIO.hpp"
00014 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00015 #include "Thyra_DefaultInverseModelEvaluator.hpp"
00016 #include "Thyra_ModelEvaluatorHelpers.hpp"
00017 #include "Thyra_DefaultFiniteDifferenceModelEvaluator.hpp"
00018 #include "Thyra_TestingTools.hpp"
00019 #include "Teuchos_GlobalMPISession.hpp"
00020 #include "Teuchos_CommandLineProcessor.hpp"
00021 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
00022 #include "Teuchos_StandardCatchMacros.hpp"
00023 #include "Teuchos_VerboseObject.hpp"
00024 #include "Teuchos_XMLParameterListHelpers.hpp"
00025 #ifdef HAVE_MPI
00026 #  include "Epetra_MpiComm.h"
00027 #else
00028 #  include "Epetra_SerialComm.h"
00029 #endif
00030 
00031 
00032 namespace {
00033 
00034 
00035 typedef AbstractLinAlgPack::value_type  Scalar;
00036 
00037 
00038 } // namespace
00039 
00040 
00041 int main( int argc, char* argv[] )
00042 {
00043 
00044   using std::endl;
00045   using Teuchos::rcp;
00046   using Teuchos::RCP;
00047   using Teuchos::rcp_dynamic_cast;
00048   using Teuchos::OSTab;
00049   using Teuchos::Array;
00050   using Teuchos::includesVerbLevel;
00051   using Teuchos::ParameterList;
00052   using Teuchos::CommandLineProcessor;
00053   typedef Teuchos::ParameterList::PrintOptions PLPO;
00054   typedef Thyra::ModelEvaluatorBase MEB;
00055   typedef Thyra::DefaultMultiVectorProductVectorSpace<Scalar> DMVPVS;
00056   using MoochoPack::MoochoSolver;
00057   using MoochoPack::MoochoThyraSolver;
00058 
00059   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00060 
00061   RCP<Epetra_Comm> epetra_comm;
00062 #ifdef HAVE_MPI
00063   epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
00064 #else
00065   epetra_comm = rcp( new Epetra_SerialComm );
00066 #endif // HAVE_MPI
00067   
00068   bool success = true, result;
00069 
00070   Teuchos::RCP<Teuchos::FancyOStream>
00071     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00072   
00073   // Rythmos uses long names so enlarge the size of this!
00074   out->setMaxLenLinePrefix(20);
00075 
00076   try {
00077   
00078     // Create the solver objects
00079     Stratimikos::DefaultLinearSolverBuilder lowsfCreator;
00080     MoochoThyraSolver moochoThyraSolver;
00081 
00082     //
00083     // A) Get options from the command line and get the parameter list
00084     //
00085 
00086     CommandLineProcessor  clp;
00087     clp.throwExceptions(false);
00088     clp.addOutputSetupOptions(true);
00089 
00090     moochoThyraSolver.setupCLP(&clp);
00091 
00092     std::string paramsFileName = "";
00093     clp.setOption( "params-file", &paramsFileName,
00094       "File name for varaious parameter lists in xml format.",
00095       true );
00096 
00097     bool useBDF = false;
00098     clp.setOption( "use-BDF", "use-BE", &useBDF,
00099       "Use BDF or Backward Euler (BE)" );
00100 
00101     double finalTime = 1e-3;
00102     clp.setOption( "final-time", &finalTime,
00103       "Final integration time (initial time is 0.0)" );
00104     
00105     int numObservationPoints= 10;
00106     clp.setOption( "num-observations", &numObservationPoints,
00107       "Number of state observations to take as basis for inversion" );
00108 
00109     double scaleParamsBy = 1e-2;
00110     clp.setOption( "scale-params-by", &scaleParamsBy,
00111       "Amount to scale the parameters by to perturb them" );
00112 
00113     Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
00114     setVerbosityLevelOption( "verb-level", &verbLevel,
00115       "Top-level verbosity level.  By default, this gets deincremented"
00116       " as you go deeper into numerical objects.",
00117       &clp );
00118 
00119     bool printValidOptions = false;
00120     clp.setOption( "print-valid-options", "no-print-valid-options",
00121       &printValidOptions, "Print the valid options or not" );
00122 
00123     bool testTransientModelEvaluator = true;
00124     clp.setOption( "test-transient-model", "no-test-transient-model", &testTransientModelEvaluator,
00125       "Test the Rythmos::ForwardSensitivityIntegratorAsModelEvaluator object or not." );
00126 
00127     double maxSensError = 1e-4;
00128     clp.setOption( "max-sens-error", &maxSensError,
00129       "The maximum allowed error in the integrated sensitivity in relation to"
00130       " the finite-difference sensitivity" );
00131 
00132     bool doInverseProblem = true;
00133     clp.setOption( "do-inv-prob", "no-do-inv-prob", &doInverseProblem,
00134       "Run and test the inverse problem or not" );
00135 
00136     double p_inv_err_tol = 1e-8;
00137     clp.setOption( "max-p-inv-error", &p_inv_err_tol,
00138       "The maximum allowed error in the inverted for parameters" );
00139     
00140     CommandLineProcessor::EParseCommandLineReturn
00141       parse_return = clp.parse(argc,argv,&std::cerr);
00142 
00143     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00144       return parse_return;
00145 
00146     moochoThyraSolver.readParameters(&*out);
00147 
00148     RCP<ParameterList>
00149       paramList = Teuchos::parameterList();
00150     if (paramsFileName.length())
00151       updateParametersFromXmlFile( paramsFileName, &*paramList );
00152     *out << "\nRead in parameters list:\n";
00153     paramList->print(*out, PLPO().indent(2).showTypes(true).showDoc(true));
00154     
00155     //
00156     // B) Create the Thyra::ModelEvaluator object for the ODE
00157     //
00158     
00159     *out << "\nCreate the EpetraExt::DiagonalTransientModel wrapper object ...\n";
00160 
00161     RCP<EpetraExt::DiagonalTransientModel>
00162       epetraStateModel = EpetraExt::diagonalTransientModel(
00163         epetra_comm,
00164         sublist(paramList,"DiagonalTransientModel",true)
00165         );
00166 
00167     if (printValidOptions) {
00168       *out <<"\nepetraStateModel valid options:\n";
00169       epetraStateModel->getValidParameters()->print(
00170         *out, PLPO().indent(2).showTypes(true).showDoc(true)
00171         );
00172     }
00173 
00174     Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
00175     linearSolverBuilder.setParameterList(sublist(paramList,"Stratimikos",true));
00176     RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
00177       W_factory = createLinearSolveStrategy(linearSolverBuilder);
00178 
00179     RCP<Thyra::ModelEvaluator<Scalar> >
00180       stateModel = epetraModelEvaluator(epetraStateModel,W_factory);
00181     
00182     //
00183     // C) Create and initialize the state stepper and state integrator objects
00184     //
00185 
00186     RCP<ParameterList>
00187       nonlinearSolverPL = sublist(paramList,"TimeStepNonlinearSolver",true);
00188     //nonlinearSolverPL->get("Default Tol",1e-3*maxStateError); // Set default if not set
00189     RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
00190       nonlinearSolver = rcp(new Rythmos::TimeStepNonlinearSolver<Scalar>());
00191     nonlinearSolver->setParameterList(nonlinearSolverPL);
00192 
00193     RCP<Rythmos::StepperBase<Scalar> > stateStepper;
00194     if (useBDF) {
00195       stateStepper = rcp(
00196         new Rythmos::ImplicitBDFStepper<Scalar>(
00197           stateModel, nonlinearSolver
00198           )
00199         );
00200     }
00201     else {
00202       stateStepper = rcp(
00203         new Rythmos::BackwardEulerStepper<Scalar>(
00204           stateModel, nonlinearSolver
00205           )
00206         );
00207     }
00208     stateStepper->setParameterList(sublist(paramList,"Rythmos Stepper",true));
00209 
00210     *out <<"\nstateStepper: " << describe(*stateStepper,verbLevel);
00211     if (printValidOptions) {
00212       *out <<"\nstateStepper valid options:\n";
00213       stateStepper->getValidParameters()->print(
00214         *out, PLPO().indent(2).showTypes(true).showDoc(true)
00215         );
00216     }
00217       
00218     RCP<Rythmos::IntegratorBase<Scalar> > integrator;
00219     {
00220       integrator = Rythmos::controlledDefaultIntegrator<Scalar>(
00221         Rythmos::simpleIntegrationControlStrategy<Scalar>(
00222           sublist(paramList,"Rythmos Integration Control",true)
00223           )
00224         );
00225     }
00226     
00227     const MEB::InArgs<Scalar>
00228       state_ic = stateModel->getNominalValues();
00229     *out << "\nstate_ic: " << describe(state_ic,verbLevel);
00230 
00231     RCP<const Thyra::VectorBase<Scalar> >
00232       p = state_ic.get_p(0);
00233     
00234     stateStepper->setInitialCondition(state_ic);
00235     integrator->setStepper(stateStepper,finalTime);
00236       
00237     if (printValidOptions) {
00238       *out <<"\nintegrator valid options:\n";
00239       integrator->getValidParameters()->print(
00240         *out, PLPO().indent(2).showTypes(true).showDoc(true)
00241         );
00242     }
00243 
00244     //
00245     // D) Generate the matching vectors by running the forward problem through
00246     //    the state integrator
00247     //
00248 
00249     Array<Scalar> observationTimes;
00250     Array<RCP<const Thyra::VectorBase<Scalar> > > stateObservations;
00251     {
00252       const double obs_dt = (finalTime - state_ic.get_t())/numObservationPoints;
00253       *out << "\nCollecting observations at intervals of obs_dt = "
00254            << obs_dt << " ...\n";
00255       OSTab tab(out);
00256       int obs_idx;
00257       double curr_t;
00258       for (
00259         obs_idx = 0, curr_t = state_ic.get_t() + obs_dt; 
00260         obs_idx < numObservationPoints;
00261         ++obs_idx, curr_t += obs_dt
00262         )
00263       {
00264         // Make sure we don't step over final time
00265         curr_t = std::min(curr_t,finalTime);
00266 
00267         *out << "\nCollecting state observation at curr_t = " << curr_t << " ...\n";
00268         observationTimes.push_back(curr_t);
00269         stateObservations.push_back(get_fwd_x(*integrator,curr_t)->clone_v());
00270         if (includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME)) {
00271           *out << "\nx("<<curr_t<<") = "
00272                << Teuchos::describe(*stateObservations.back(),Teuchos::VERB_EXTREME);
00273         }
00274       }
00275     }
00276     
00277     //
00278     // F) Create the forward sensitivity stepper and reduced sensitivity model
00279     //    evaluator
00280     //
00281     
00282     //
00283     // F.1) Create the forward sensitivity stepper
00284     //
00285     
00286     RCP<Rythmos::ForwardSensitivityStepper<Scalar> >
00287       stateAndSensStepper = Rythmos::forwardSensitivityStepper<Scalar>(
00288         stateModel, 0, stateModel->getNominalValues(),
00289         stateStepper, nonlinearSolver
00290         );
00291     // The above call will result in stateStepper and nonlinearSolver being
00292     // cloned.  This helps to ensure consistency between the state and
00293     // sensitivity computations!
00294 
00295     const RCP<const DMVPVS>
00296       s_bar_space = rcp_dynamic_cast<const DMVPVS>(
00297         stateAndSensStepper->getFwdSensModel()->get_x_space()
00298         );
00299 
00300     //
00301     // F.2) Set the initial condition for the state and forward sensitivities
00302     //
00303 
00304     RCP<Thyra::VectorBase<Scalar> > s_bar_init
00305       = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
00306     assign( &*s_bar_init, 0.0 );
00307     RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init
00308       = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
00309     assign( &*s_bar_dot_init, 0.0 );
00310     // Above, I believe that these are the correct initial conditions for
00311     // s_bar and s_bar_dot given how the EpetraExt::DiagonalTransientModel
00312     // is currently implemented!
00313 
00314     RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
00315       stateAndSensModel = stateAndSensStepper->getStateAndFwdSensModel();
00316 
00317     MEB::InArgs<Scalar>
00318       state_and_sens_ic = stateAndSensStepper->getModel()->createInArgs();
00319 
00320     // Copy time, parameters etc.
00321     state_and_sens_ic.setArgs(state_ic);
00322     // Set initial condition for x_bar = [ x; s_bar ]
00323     state_and_sens_ic.set_x(
00324       stateAndSensModel->create_x_bar_vec(state_ic.get_x(),s_bar_init)
00325       );
00326     // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
00327     state_and_sens_ic.set_x_dot(
00328       stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(),s_bar_dot_init)
00329       );
00330     
00331     *out << "\nstate_and_sens_ic:\n" << describe(state_and_sens_ic,verbLevel);
00332     
00333     stateAndSensStepper->setInitialCondition(state_and_sens_ic);
00334 
00335     //
00336     // F.3) Setup the state-matching response functions
00337     //
00338 
00339     Array<RCP<const Thyra::ModelEvaluator<Scalar> > > matchingFuncs;
00340     Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > matchingFuncBasePoints;
00341 
00342     for (int obs_idx = 0; obs_idx < numObservationPoints; ++obs_idx ) {
00343       RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
00344         matchingFunc = Thyra::defaultInverseModelEvaluator(stateModel);
00345       matchingFunc->setParameterList(
00346         sublist(paramList,"DefaultInverseModelEvalautor",true));
00347       matchingFunc->set_observationTarget(stateObservations[obs_idx]);
00348       matchingFuncs.push_back(matchingFunc);
00349       matchingFuncBasePoints.push_back(state_ic);
00350     }
00351 
00352     //
00353     // F.4) Setup the final integrator that will evaluate the reduced response
00354     //      function and it's derivatives
00355     
00356     namespace FSIAMET = Rythmos::ForwardSensitivityIntegratorAsModelEvaluatorTypes;
00357     const RCP<Rythmos::IntegratorBase<Scalar> > nullIntegrator;
00358     RCP<Rythmos::ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
00359       sensIntegatorAsModelEvaluator =
00360       Rythmos::forwardSensitivityIntegratorAsModelEvaluator<Scalar>(
00361         stateStepper, integrator, stateAndSensStepper, nullIntegrator,
00362         state_and_sens_ic, observationTimes, matchingFuncs,
00363         matchingFuncBasePoints, FSIAMET::RESPONSE_TYPE_SUM
00364         );
00365     sensIntegatorAsModelEvaluator->setParameterList(
00366       sublist(paramList,"ForwardSensitivityIntegratorAsModelEvaluator",true) );
00367     
00368     //
00369     // G) Test the ForwardSensitivityIntegratorAsModelEvaluator object
00370     //
00371 
00372     if (testTransientModelEvaluator) {
00373 
00374       //
00375       // G.1) Test that the base evaluation computes a zero response function.
00376       //
00377 
00378       // Here, we check that if we evaluation p -> g(p) for the base p, then we
00379       // should get zero since the state should match the read observations
00380       // exactly.  We also check that the computed state is exactly the same
00381       // also.  We can cheat and get the final state from stateStepper which is
00382       // used internally.
00383       //
00384 
00385       *out << "\nTest that sensIntegatorAsModelEvaluator computes the right base point ... \n";
00386       {
00387       
00388         RCP<Thyra::VectorBase<Scalar> >
00389           d_hat = createMember(sensIntegatorAsModelEvaluator->get_g_space(0));
00390 
00391         eval_g(
00392           *sensIntegatorAsModelEvaluator,
00393           0, *state_ic.get_p(0),
00394           0, &*d_hat
00395           );
00396 
00397         *out
00398           << "\nd_hat(p) evaluated using sensIntegatorAsModelEvaluator:\n"
00399           << *d_hat;
00400 
00401         // Test that d_hat == 0.0 to all digits since 
00402 
00403         const Scalar d_hat_nrm = norm(*d_hat);
00404 
00405         const bool d_hat_norm_is_zero = (d_hat_nrm == 0.0);
00406 
00407         *out << "\nCheck: ||d_hat|| = " << d_hat_nrm << " == 0 : "
00408              << Thyra::passfail(d_hat_norm_is_zero) << endl;
00409 
00410         if (!d_hat_norm_is_zero)
00411           success = false;
00412 
00413         // Test that:
00414         //
00415         //  get_x(stateStepper,finalTime) == stateObservations[numObservationPoints-1]
00416         //
00417         // to the given precision.
00418 
00419         *out
00420           << "\nChecking that get_x(stateStepper,finalTime) == stateObservations[numObservationPoints-1] ... \n";
00421 
00422         RCP<const Thyra::VectorBase<Scalar> >
00423           x_final = get_x(*stateStepper,observationTimes[numObservationPoints-1]);
00424       
00425         *out << endl;
00426         result = Thyra::testRelNormDiffErr(
00427           "get_x(stateStepper,finalTime)", *x_final,
00428           "stateObservations[numObservationPoints-1]", *stateObservations[numObservationPoints-1],
00429           "errorTol", 0.0, "warningTol", 1e-8,
00430           &OSTab(out).o(), verbLevel
00431           );
00432         if (!result) success = false;
00433 
00434       }
00435 
00436       *out << "\nBase parameters p = " << *p;
00437 
00438       RCP<Thyra::VectorBase<Scalar> >
00439         p_perturb = createMember(p->space());
00440 
00441       V_StV( &*p_perturb, scaleParamsBy, *p );
00442 
00443       *out << "\nPerturbed parameters p_perturb = " << *p_perturb;
00444 
00445       //
00446       // G.2) Test the sensitivity D(d_hat)/d(p) by finite differences
00447       //
00448 
00449       //
00450       *out << "\nCompute analytic D(d_hat)/d(p) at p_perturb ...\n";
00451       //
00452 
00453       MEB::DerivativeMultiVector<Scalar> D_d_hat_D_p_exact;
00454       {
00455 
00456         D_d_hat_D_p_exact = Thyra::create_DgDp_mv(
00457           *sensIntegatorAsModelEvaluator, 0, 0, MEB::DERIV_TRANS_MV_BY_ROW );
00458 
00459         MEB::InArgs<Scalar>
00460           inArgs = sensIntegatorAsModelEvaluator->createInArgs();
00461         inArgs.set_p(0,p_perturb);
00462 
00463         MEB::OutArgs<Scalar>
00464           outArgs = sensIntegatorAsModelEvaluator->createOutArgs();
00465         outArgs.set_DgDp(0,0,D_d_hat_D_p_exact);
00466 
00467         sensIntegatorAsModelEvaluator->evalModel(inArgs,outArgs);
00468       
00469         *out
00470           << "\nAnalytic D(d_hat)/d(p)^T = "
00471           << describe(*D_d_hat_D_p_exact.getMultiVector(),verbLevel);
00472 
00473       }
00474  
00475       //
00476       *out << "\nCompute finite-diff D(d_hat)/d(p) at p_perturb ...\n";
00477       //
00478 
00479       MEB::DerivativeMultiVector<Scalar> D_d_hat_D_p_fd;
00480       {
00481 
00482         RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > fdCalc =
00483           Thyra::directionalFiniteDiffCalculator<Scalar>(
00484             sublist(paramList,"DirectionalFiniteDiffCalculator",true)
00485             );
00486 
00487         RCP<Thyra::DefaultFiniteDifferenceModelEvaluator<Scalar> > fdModelEval =
00488           Thyra::defaultFiniteDifferenceModelEvaluator<Scalar>(
00489             sensIntegatorAsModelEvaluator,
00490             fdCalc
00491             );
00492 
00493         MEB::InArgs<Scalar>
00494           inArgs = fdModelEval->createInArgs();
00495         inArgs.set_p(0,p_perturb);
00496       
00497         D_d_hat_D_p_fd = Thyra::create_DgDp_mv(
00498           *fdModelEval, 0, 0, MEB::DERIV_TRANS_MV_BY_ROW
00499           );
00500       
00501         MEB::OutArgs<Scalar>
00502           outArgs = fdModelEval->createOutArgs();
00503         outArgs.set_DgDp( 0, 0, D_d_hat_D_p_fd );
00504       
00505         // Silence the model evaluators that are called.  The fdCal object
00506         // will show all of the inputs and outputs for each call.
00507         stateStepper->setVerbLevel(Teuchos::VERB_NONE);
00508         sensIntegatorAsModelEvaluator->setVerbLevel(Teuchos::VERB_NONE);
00509 
00510         fdModelEval->evalModel(inArgs,outArgs);
00511       
00512         *out
00513           << "\nFinite difference D(d_hat)/d(p)^T = "
00514           << describe(*D_d_hat_D_p_fd.getMultiVector(),verbLevel);
00515       
00516       }
00517     
00518       //
00519       *out << "\nCompare analytic and finite-diff D(d_hat)/d(p)^T ...\n";
00520       //
00521 
00522       RCP<const Thyra::VectorBase<Scalar> >
00523         D_d_hat_D_p_exact_vec = D_d_hat_D_p_exact.getMultiVector()->col(0),
00524         D_d_hat_D_p_fd_vec = D_d_hat_D_p_fd.getMultiVector()->col(0);
00525     
00526       result = Thyra::testRelNormDiffErr(
00527         "D_d_hat_D_p_exact_vec", *D_d_hat_D_p_exact_vec,
00528         "D_d_hat_D_p_fd_vec", *D_d_hat_D_p_fd_vec,
00529         "maxSensError", maxSensError,
00530         "warningTol", 1.0, // Don't warn
00531         &*out, verbLevel
00532         );
00533       if (!result) success = false;
00534 
00535     }
00536 
00537     if (doInverseProblem) {
00538     
00539       //
00540       // H) Setup the MoochoThyraSolver object
00541       //
00542       
00543       *out << "\nSetup the MoochoThyraSolver object ...\n";
00544       
00545       moochoThyraSolver.setParameterList(
00546         sublist(paramList,"MoochoThyraSolver",true) );
00547       
00548       moochoThyraSolver.setModel(sensIntegatorAsModelEvaluator);
00549       
00550       //
00551       // I) Solve the transient inverse problem using MOOCHO.
00552       //
00553       
00554       *out << "\nSetup the transient inverse problem ...\n";
00555       
00556       moochoThyraSolver.readInitialGuess(&*out);
00557       
00558       const MoochoSolver::ESolutionStatus
00559         solution_status = moochoThyraSolver.solve();
00560       
00561       if ( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED )
00562         success = false;
00563       
00564       moochoThyraSolver.writeFinalSolution(&*out);
00565       
00566       //
00567       // J) Check the inverse solution against the true solution
00568       //
00569       
00570       RCP<const Thyra::VectorBase<Scalar> >
00571         p_inv = moochoThyraSolver.getFinalPoint().get_p(0);
00572       
00573       result = Thyra::testRelNormDiffErr(
00574         "p_inv", *p_inv,
00575         "p", *p,
00576         "errorTol", p_inv_err_tol,
00577         "warningTol", p_inv_err_tol * 1e+2,
00578         &OSTab(out).o(), Teuchos::VERB_EXTREME
00579         );
00580       if (!result) success = false;
00581 
00582     }
00583       
00584   }
00585   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
00586 
00587   if (success)
00588     *out << "\nEnd Result: TEST PASSED\n";
00589   else
00590     *out << "\nEnd Result: TEST FAILED\n";
00591   
00592   return ( success ? 0 : 1 );
00593 
00594 }

Generated on Wed May 12 21:52:29 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7