MOOCHO  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DiagonalTransientInversionMain.cpp
1 /*
2 // @HEADER
3 // ***********************************************************************
4 //
5 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
6 // Copyright (2003) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 */
43 
44 
46 #include "MoochoPack_MoochoThyraSolver.hpp"
47 #include "Rythmos_BackwardEulerStepper.hpp"
48 #include "Rythmos_ImplicitBDFStepper.hpp"
49 #include "Rythmos_ForwardSensitivityStepper.hpp"
50 #include "Rythmos_TimeStepNonlinearSolver.hpp"
51 #include "Rythmos_DefaultIntegrator.hpp"
52 #include "Rythmos_SimpleIntegrationControlStrategy.hpp"
53 #include "Rythmos_StepperAsModelEvaluator.hpp"
54 #include "Rythmos_ForwardSensitivityIntegratorAsModelEvaluator.hpp"
55 #include "Thyra_EpetraModelEvaluator.hpp"
56 #include "Thyra_DefaultSpmdMultiVectorFileIO.hpp"
57 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
58 #include "Thyra_DefaultInverseModelEvaluator.hpp"
59 #include "Thyra_ModelEvaluatorHelpers.hpp"
60 #include "Thyra_DefaultFiniteDifferenceModelEvaluator.hpp"
61 #include "Thyra_TestingTools.hpp"
62 #include "Teuchos_GlobalMPISession.hpp"
63 #include "Teuchos_CommandLineProcessor.hpp"
64 #include "Teuchos_VerbosityLevelCommandLineProcessorHelpers.hpp"
65 #include "Teuchos_StandardCatchMacros.hpp"
66 #include "Teuchos_VerboseObject.hpp"
67 #include "Teuchos_XMLParameterListHelpers.hpp"
68 #ifdef HAVE_MPI
69 # include "Epetra_MpiComm.h"
70 #else
71 # include "Epetra_SerialComm.h"
72 #endif
73 
74 
75 namespace {
76 
77 
78 typedef AbstractLinAlgPack::value_type Scalar;
79 
80 
81 } // namespace
82 
83 
84 int main( int argc, char* argv[] )
85 {
86 
87  using std::endl;
88  using Teuchos::rcp;
89  using Teuchos::RCP;
90  using Teuchos::rcp_dynamic_cast;
91  using Teuchos::OSTab;
92  using Teuchos::Array;
97  typedef Thyra::ModelEvaluatorBase MEB;
101 
102  Teuchos::GlobalMPISession mpiSession(&argc,&argv);
103 
104  RCP<Epetra_Comm> epetra_comm;
105 #ifdef HAVE_MPI
106  epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
107 #else
108  epetra_comm = rcp( new Epetra_SerialComm );
109 #endif // HAVE_MPI
110 
111  bool success = true, result;
112 
115 
116  // Rythmos uses long names so enlarge the size of this!
117  out->setMaxLenLinePrefix(20);
118 
119  try {
120 
121  // Create the solver objects
123  MoochoThyraSolver moochoThyraSolver;
124 
125  //
126  // A) Get options from the command line and get the parameter list
127  //
128 
129  CommandLineProcessor clp;
130  clp.throwExceptions(false);
131  clp.addOutputSetupOptions(true);
132 
133  moochoThyraSolver.setupCLP(&clp);
134 
135  std::string paramsFileName = "";
136  clp.setOption( "params-file", &paramsFileName,
137  "File name for varaious parameter lists in xml format.",
138  true );
139 
140  bool useBDF = false;
141  clp.setOption( "use-BDF", "use-BE", &useBDF,
142  "Use BDF or Backward Euler (BE)" );
143 
144  double finalTime = 1e-3;
145  clp.setOption( "final-time", &finalTime,
146  "Final integration time (initial time is 0.0)" );
147 
148  int numObservationPoints= 10;
149  clp.setOption( "num-observations", &numObservationPoints,
150  "Number of state observations to take as basis for inversion" );
151 
152  double scaleParamsBy = 1e-2;
153  clp.setOption( "scale-params-by", &scaleParamsBy,
154  "Amount to scale the parameters by to perturb them" );
155 
157  setVerbosityLevelOption( "verb-level", &verbLevel,
158  "Top-level verbosity level. By default, this gets deincremented"
159  " as you go deeper into numerical objects.",
160  &clp );
161 
162  bool printValidOptions = false;
163  clp.setOption( "print-valid-options", "no-print-valid-options",
164  &printValidOptions, "Print the valid options or not" );
165 
166  bool testTransientModelEvaluator = true;
167  clp.setOption( "test-transient-model", "no-test-transient-model", &testTransientModelEvaluator,
168  "Test the Rythmos::ForwardSensitivityIntegratorAsModelEvaluator object or not." );
169 
170  double maxSensError = 1e-4;
171  clp.setOption( "max-sens-error", &maxSensError,
172  "The maximum allowed error in the integrated sensitivity in relation to"
173  " the finite-difference sensitivity" );
174 
175  bool doInverseProblem = true;
176  clp.setOption( "do-inv-prob", "no-do-inv-prob", &doInverseProblem,
177  "Run and test the inverse problem or not" );
178 
179  double p_inv_err_tol = 1e-8;
180  clp.setOption( "max-p-inv-error", &p_inv_err_tol,
181  "The maximum allowed error in the inverted for parameters" );
182 
183  CommandLineProcessor::EParseCommandLineReturn
184  parse_return = clp.parse(argc,argv,&std::cerr);
185 
186  if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
187  return parse_return;
188 
189  moochoThyraSolver.readParameters(&*out);
190 
191  RCP<ParameterList>
192  paramList = Teuchos::parameterList();
193  if (paramsFileName.length())
194  updateParametersFromXmlFile( paramsFileName, paramList.ptr() );
195  *out << "\nRead in parameters list:\n";
196  paramList->print(*out, PLPO().indent(2).showTypes(true).showDoc(true));
197 
198  //
199  // B) Create the Thyra::ModelEvaluator object for the ODE
200  //
201 
202  *out << "\nCreate the EpetraExt::DiagonalTransientModel wrapper object ...\n";
203 
204  RCP<EpetraExt::DiagonalTransientModel>
205  epetraStateModel = EpetraExt::diagonalTransientModel(
206  epetra_comm,
207  sublist(paramList,"DiagonalTransientModel",true)
208  );
209 
210  if (printValidOptions) {
211  *out <<"\nepetraStateModel valid options:\n";
212  epetraStateModel->getValidParameters()->print(
213  *out, PLPO().indent(2).showTypes(true).showDoc(true)
214  );
215  }
216 
217  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
218  linearSolverBuilder.setParameterList(sublist(paramList,"Stratimikos",true));
219  RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> >
220  W_factory = createLinearSolveStrategy(linearSolverBuilder);
221 
222  RCP<Thyra::ModelEvaluator<Scalar> >
223  stateModel = epetraModelEvaluator(epetraStateModel,W_factory);
224 
225  //
226  // C) Create and initialize the state stepper and state integrator objects
227  //
228 
229  RCP<ParameterList>
230  nonlinearSolverPL = sublist(paramList,"TimeStepNonlinearSolver",true);
231  //nonlinearSolverPL->get("Default Tol",1e-3*maxStateError); // Set default if not set
232  RCP<Rythmos::TimeStepNonlinearSolver<Scalar> >
233  nonlinearSolver = rcp(new Rythmos::TimeStepNonlinearSolver<Scalar>());
234  nonlinearSolver->setParameterList(nonlinearSolverPL);
235 
236  RCP<Rythmos::StepperBase<Scalar> > stateStepper;
237  if (useBDF) {
238  stateStepper = rcp(
239  new Rythmos::ImplicitBDFStepper<Scalar>(
240  stateModel, nonlinearSolver
241  )
242  );
243  }
244  else {
245  stateStepper = rcp(
246  new Rythmos::BackwardEulerStepper<Scalar>(
247  stateModel, nonlinearSolver
248  )
249  );
250  }
251  stateStepper->setParameterList(sublist(paramList,"Rythmos Stepper",true));
252 
253  *out <<"\nstateStepper: " << describe(*stateStepper,verbLevel);
254  if (printValidOptions) {
255  *out <<"\nstateStepper valid options:\n";
256  stateStepper->getValidParameters()->print(
257  *out, PLPO().indent(2).showTypes(true).showDoc(true)
258  );
259  }
260 
261  RCP<Rythmos::IntegratorBase<Scalar> > integrator;
262  {
263  integrator = Rythmos::controlledDefaultIntegrator<Scalar>(
264  Rythmos::simpleIntegrationControlStrategy<Scalar>(
265  sublist(paramList,"Rythmos Integration Control",true)
266  )
267  );
268  }
269 
270  const MEB::InArgs<Scalar>
271  state_ic = stateModel->getNominalValues();
272  *out << "\nstate_ic: " << describe(state_ic,verbLevel);
273 
274  RCP<const Thyra::VectorBase<Scalar> >
275  p = state_ic.get_p(0);
276 
277  stateStepper->setInitialCondition(state_ic);
278  integrator->setStepper(stateStepper,finalTime);
279 
280  if (printValidOptions) {
281  *out <<"\nintegrator valid options:\n";
282  integrator->getValidParameters()->print(
283  *out, PLPO().indent(2).showTypes(true).showDoc(true)
284  );
285  }
286 
287  //
288  // D) Generate the matching vectors by running the forward problem through
289  // the state integrator
290  //
291 
292  Array<Scalar> observationTimes;
293  Array<RCP<const Thyra::VectorBase<Scalar> > > stateObservations;
294  {
295  const double obs_dt = (finalTime - state_ic.get_t())/numObservationPoints;
296  *out << "\nCollecting observations at intervals of obs_dt = "
297  << obs_dt << " ...\n";
298  OSTab tab(out);
299  int obs_idx;
300  double curr_t;
301  for (
302  obs_idx = 0, curr_t = state_ic.get_t() + obs_dt;
303  obs_idx < numObservationPoints;
304  ++obs_idx, curr_t += obs_dt
305  )
306  {
307  // Make sure we don't step over final time
308  curr_t = std::min(curr_t,finalTime);
309 
310  *out << "\nCollecting state observation at curr_t = " << curr_t << " ...\n";
311  observationTimes.push_back(curr_t);
312  stateObservations.push_back(get_fwd_x(*integrator,curr_t)->clone_v());
313  if (includesVerbLevel(verbLevel,Teuchos::VERB_EXTREME)) {
314  *out << "\nx("<<curr_t<<") = "
315  << Teuchos::describe(*stateObservations.back(),Teuchos::VERB_EXTREME);
316  }
317  }
318  }
319 
320  //
321  // F) Create the forward sensitivity stepper and reduced sensitivity model
322  // evaluator
323  //
324 
325  //
326  // F.1) Create the forward sensitivity stepper
327  //
328 
329  RCP<Rythmos::ForwardSensitivityStepper<Scalar> >
330  stateAndSensStepper = Rythmos::forwardSensitivityStepper<Scalar>(
331  stateModel, 0, stateModel->getNominalValues(),
332  stateStepper, nonlinearSolver
333  );
334  // The above call will result in stateStepper and nonlinearSolver being
335  // cloned. This helps to ensure consistency between the state and
336  // sensitivity computations!
337 
338  const RCP<const DMVPVS>
339  s_bar_space = rcp_dynamic_cast<const DMVPVS>(
340  stateAndSensStepper->getFwdSensModel()->get_x_space()
341  );
342 
343  //
344  // F.2) Set the initial condition for the state and forward sensitivities
345  //
346 
347  RCP<Thyra::VectorBase<Scalar> > s_bar_init
348  = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
349  assign( s_bar_init.ptr(), 0.0 );
350  RCP<Thyra::VectorBase<Scalar> > s_bar_dot_init
351  = createMember(stateAndSensStepper->getFwdSensModel()->get_x_space());
352  assign( s_bar_dot_init.ptr(), 0.0 );
353  // Above, I believe that these are the correct initial conditions for
354  // s_bar and s_bar_dot given how the EpetraExt::DiagonalTransientModel
355  // is currently implemented!
356 
357  RCP<const Rythmos::StateAndForwardSensitivityModelEvaluator<Scalar> >
358  stateAndSensModel = stateAndSensStepper->getStateAndFwdSensModel();
359 
360  MEB::InArgs<Scalar>
361  state_and_sens_ic = stateAndSensStepper->getModel()->createInArgs();
362 
363  // Copy time, parameters etc.
364  state_and_sens_ic.setArgs(state_ic);
365  // Set initial condition for x_bar = [ x; s_bar ]
366  state_and_sens_ic.set_x(
367  stateAndSensModel->create_x_bar_vec(state_ic.get_x(),s_bar_init)
368  );
369  // Set initial condition for x_bar_dot = [ x_dot; s_bar_dot ]
370  state_and_sens_ic.set_x_dot(
371  stateAndSensModel->create_x_bar_vec(state_ic.get_x_dot(),s_bar_dot_init)
372  );
373 
374  *out << "\nstate_and_sens_ic:\n" << describe(state_and_sens_ic,verbLevel);
375 
376  stateAndSensStepper->setInitialCondition(state_and_sens_ic);
377 
378  //
379  // F.3) Setup the state-matching response functions
380  //
381 
382  Array<RCP<const Thyra::ModelEvaluator<Scalar> > > matchingFuncs;
383  Array<Thyra::ModelEvaluatorBase::InArgs<Scalar> > matchingFuncBasePoints;
384 
385  for (int obs_idx = 0; obs_idx < numObservationPoints; ++obs_idx ) {
386  RCP<Thyra::DefaultInverseModelEvaluator<Scalar> >
387  matchingFunc = Thyra::defaultInverseModelEvaluator(stateModel);
388  matchingFunc->setParameterList(
389  sublist(paramList,"DefaultInverseModelEvalautor",true));
390  matchingFunc->set_observationTarget(stateObservations[obs_idx]);
391  matchingFuncs.push_back(matchingFunc);
392  matchingFuncBasePoints.push_back(state_ic);
393  }
394 
395  //
396  // F.4) Setup the final integrator that will evaluate the reduced response
397  // function and it's derivatives
398 
399  namespace FSIAMET = Rythmos::ForwardSensitivityIntegratorAsModelEvaluatorTypes;
400  const RCP<Rythmos::IntegratorBase<Scalar> > nullIntegrator;
401  RCP<Rythmos::ForwardSensitivityIntegratorAsModelEvaluator<Scalar> >
402  sensIntegatorAsModelEvaluator =
403  Rythmos::forwardSensitivityIntegratorAsModelEvaluator<Scalar>(
404  stateStepper, integrator, stateAndSensStepper, nullIntegrator,
405  state_and_sens_ic, observationTimes, matchingFuncs,
406  matchingFuncBasePoints, FSIAMET::RESPONSE_TYPE_SUM
407  );
408  sensIntegatorAsModelEvaluator->setParameterList(
409  sublist(paramList,"ForwardSensitivityIntegratorAsModelEvaluator",true) );
410 
411  //
412  // G) Test the ForwardSensitivityIntegratorAsModelEvaluator object
413  //
414 
415  if (testTransientModelEvaluator) {
416 
417  //
418  // G.1) Test that the base evaluation computes a zero response function.
419  //
420 
421  // Here, we check that if we evaluation p -> g(p) for the base p, then we
422  // should get zero since the state should match the read observations
423  // exactly. We also check that the computed state is exactly the same
424  // also. We can cheat and get the final state from stateStepper which is
425  // used internally.
426  //
427 
428  *out << "\nTest that sensIntegatorAsModelEvaluator computes the right base point ... \n";
429  {
430 
431  RCP<Thyra::VectorBase<Scalar> >
432  d_hat = createMember(sensIntegatorAsModelEvaluator->get_g_space(0));
433 
434  eval_g(
435  *sensIntegatorAsModelEvaluator,
436  0, *state_ic.get_p(0),
437  0, d_hat.ptr()
438  );
439 
440  *out
441  << "\nd_hat(p) evaluated using sensIntegatorAsModelEvaluator:\n"
442  << *d_hat;
443 
444  // Test that d_hat == 0.0 to all digits since
445 
446  const Scalar d_hat_nrm = norm(*d_hat);
447 
448  const bool d_hat_norm_is_zero = (d_hat_nrm == 0.0);
449 
450  *out << "\nCheck: ||d_hat|| = " << d_hat_nrm << " == 0 : "
451  << Thyra::passfail(d_hat_norm_is_zero) << endl;
452 
453  if (!d_hat_norm_is_zero)
454  success = false;
455 
456  // Test that:
457  //
458  // get_x(stateStepper,finalTime) == stateObservations[numObservationPoints-1]
459  //
460  // to the given precision.
461 
462  *out
463  << "\nChecking that get_x(stateStepper,finalTime) == stateObservations[numObservationPoints-1] ... \n";
464 
465  RCP<const Thyra::VectorBase<Scalar> >
466  x_final = get_x(*stateStepper,observationTimes[numObservationPoints-1]);
467 
468  *out << endl;
469  result = Thyra::testRelNormDiffErr(
470  "get_x(stateStepper,finalTime)", *x_final,
471  "stateObservations[numObservationPoints-1]", *stateObservations[numObservationPoints-1],
472  "errorTol", 0.0, "warningTol", 1e-8,
473  &OSTab(out).o(), verbLevel
474  );
475  if (!result) success = false;
476 
477  }
478 
479  *out << "\nBase parameters p = " << *p;
480 
481  RCP<Thyra::VectorBase<Scalar> >
482  p_perturb = createMember(p->space());
483 
484  V_StV( p_perturb.ptr(), scaleParamsBy, *p );
485 
486  *out << "\nPerturbed parameters p_perturb = " << *p_perturb;
487 
488  //
489  // G.2) Test the sensitivity D(d_hat)/d(p) by finite differences
490  //
491 
492  //
493  *out << "\nCompute analytic D(d_hat)/d(p) at p_perturb ...\n";
494  //
495 
496  MEB::DerivativeMultiVector<Scalar> D_d_hat_D_p_exact;
497  {
498 
499  D_d_hat_D_p_exact = Thyra::create_DgDp_mv(
500  *sensIntegatorAsModelEvaluator, 0, 0, MEB::DERIV_TRANS_MV_BY_ROW );
501 
502  MEB::InArgs<Scalar>
503  inArgs = sensIntegatorAsModelEvaluator->createInArgs();
504  inArgs.set_p(0,p_perturb);
505 
506  MEB::OutArgs<Scalar>
507  outArgs = sensIntegatorAsModelEvaluator->createOutArgs();
508  outArgs.set_DgDp(0,0,D_d_hat_D_p_exact);
509 
510  sensIntegatorAsModelEvaluator->evalModel(inArgs,outArgs);
511 
512  *out
513  << "\nAnalytic D(d_hat)/d(p)^T = "
514  << describe(*D_d_hat_D_p_exact.getMultiVector(),verbLevel);
515 
516  }
517 
518  //
519  *out << "\nCompute finite-diff D(d_hat)/d(p) at p_perturb ...\n";
520  //
521 
522  MEB::DerivativeMultiVector<Scalar> D_d_hat_D_p_fd;
523  {
524 
525  RCP<Thyra::DirectionalFiniteDiffCalculator<Scalar> > fdCalc =
526  Thyra::directionalFiniteDiffCalculator<Scalar>(
527  sublist(paramList,"DirectionalFiniteDiffCalculator",true)
528  );
529 
530  RCP<Thyra::DefaultFiniteDifferenceModelEvaluator<Scalar> > fdModelEval =
531  Thyra::defaultFiniteDifferenceModelEvaluator<Scalar>(
532  sensIntegatorAsModelEvaluator,
533  fdCalc
534  );
535 
536  MEB::InArgs<Scalar>
537  inArgs = fdModelEval->createInArgs();
538  inArgs.set_p(0,p_perturb);
539 
540  D_d_hat_D_p_fd = Thyra::create_DgDp_mv(
541  *fdModelEval, 0, 0, MEB::DERIV_TRANS_MV_BY_ROW
542  );
543 
544  MEB::OutArgs<Scalar>
545  outArgs = fdModelEval->createOutArgs();
546  outArgs.set_DgDp( 0, 0, D_d_hat_D_p_fd );
547 
548  // Silence the model evaluators that are called. The fdCal object
549  // will show all of the inputs and outputs for each call.
550  stateStepper->setVerbLevel(Teuchos::VERB_NONE);
551  sensIntegatorAsModelEvaluator->setVerbLevel(Teuchos::VERB_NONE);
552 
553  fdModelEval->evalModel(inArgs,outArgs);
554 
555  *out
556  << "\nFinite difference D(d_hat)/d(p)^T = "
557  << describe(*D_d_hat_D_p_fd.getMultiVector(),verbLevel);
558 
559  }
560 
561  //
562  *out << "\nCompare analytic and finite-diff D(d_hat)/d(p)^T ...\n";
563  //
564 
565  RCP<const Thyra::VectorBase<Scalar> >
566  D_d_hat_D_p_exact_vec = D_d_hat_D_p_exact.getMultiVector()->col(0),
567  D_d_hat_D_p_fd_vec = D_d_hat_D_p_fd.getMultiVector()->col(0);
568 
569  result = Thyra::testRelNormDiffErr(
570  "D_d_hat_D_p_exact_vec", *D_d_hat_D_p_exact_vec,
571  "D_d_hat_D_p_fd_vec", *D_d_hat_D_p_fd_vec,
572  "maxSensError", maxSensError,
573  "warningTol", 1.0, // Don't warn
574  &*out, verbLevel
575  );
576  if (!result) success = false;
577 
578  }
579 
580  if (doInverseProblem) {
581 
582  //
583  // H) Setup the MoochoThyraSolver object
584  //
585 
586  *out << "\nSetup the MoochoThyraSolver object ...\n";
587 
588  moochoThyraSolver.setParameterList(
589  sublist(paramList,"MoochoThyraSolver",true) );
590 
591  moochoThyraSolver.setModel(sensIntegatorAsModelEvaluator);
592 
593  //
594  // I) Solve the transient inverse problem using MOOCHO.
595  //
596 
597  *out << "\nSetup the transient inverse problem ...\n";
598 
599  moochoThyraSolver.readInitialGuess(&*out);
600 
601  const MoochoSolver::ESolutionStatus
602  solution_status = moochoThyraSolver.solve();
603 
604  if ( solution_status != MoochoSolver::SOLVE_RETURN_SOLVED )
605  success = false;
606 
607  moochoThyraSolver.writeFinalSolution(&*out);
608 
609  //
610  // J) Check the inverse solution against the true solution
611  //
612 
613  RCP<const Thyra::VectorBase<Scalar> >
614  p_inv = moochoThyraSolver.getFinalPoint().get_p(0);
615 
616  result = Thyra::testRelNormDiffErr(
617  "p_inv", *p_inv,
618  "p", *p,
619  "errorTol", p_inv_err_tol,
620  "warningTol", p_inv_err_tol * 1e+2,
621  &OSTab(out).o(), Teuchos::VERB_EXTREME
622  );
623  if (!result) success = false;
624 
625  }
626 
627  }
628  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
629 
630  if (success)
631  *out << "\nEnd Result: TEST PASSED\n";
632  else
633  *out << "\nEnd Result: TEST FAILED\n";
634 
635  return ( success ? 0 : 1 );
636 
637 }
basic_OSTab< char > OSTab
void setupCLP(Teuchos::CommandLineProcessor *clp)
basic_FancyOStream & setMaxLenLinePrefix(const int maxLenLinePrefix)
const std::string passfail(const bool result)
bool testRelNormDiffErr(const std::string &v1_name, const VectorBase< Scalar > &v1, const std::string &v2_name, const VectorBase< Scalar > &v2, const std::string &maxRelErr_error_name, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType &maxRelErr_error, const std::string &maxRelErr_warning_name, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType &maxRelErr_warning, std::ostream *out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_LOW, const std::string &leadingIndent=std::string(""))
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
void setParameterList(RCP< ParameterList > const &paramList)
void assign(MatrixOp *M_lhs, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< FancyOStream > getDefaultOStream()
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)