void DiagonalROME<Scalar>::setDiagonalBarVector(
    const RCP<const Thyra::VectorBase<Scalar> > &diag_bar)
{

    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::rcp_dynamic_cast;
    using Thyra::createMember;
    using Thyra::ele_wise_divide;
    using Thyra::V_S;

    diag_bar_ = diag_bar.assert_not_null();

    // Reset the scalar product for p_space!

    RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);

    // s_bar[i] = diag[i] / diag_bar[i]
    V_S( s_bar.ptr(), ST::zero() );
    ele_wise_divide( ST::one(), *diag_, *diag_bar_, s_bar.ptr() );
    s_bar_ = s_bar;

    const RCP<Thyra::ScalarProdVectorSpaceBase<Scalar> > sp_p_space =
        rcp_dynamic_cast<Thyra::ScalarProdVectorSpaceBase<Scalar> >(p_space_, true);
    //sp_p_space->setScalarProd(diagonalScalarProd<Scalar>(s_bar_));

}

bool sillyCgSolve(
  const Thyra::LinearOpBase<Scalar> &A,
  const Thyra::VectorBase<Scalar> &b,
  const int maxNumIters,
  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance,
  const Teuchos::Ptr<Thyra::VectorBase<Scalar> > &x,
  std::ostream &out
  )
{

  // Create some typedefs and some other stuff to make the code cleaner
  typedef Teuchos::ScalarTraits<Scalar> ST; typedef typename ST::magnitudeType ScalarMag;
  const Scalar one = ST::one(), zero = ST::zero();  using Teuchos::as;
  using Teuchos::RCP; using Thyra::VectorSpaceBase; using Thyra::VectorBase;
  using Thyra::NOTRANS; using Thyra::V_V; using Thyra::apply;
  

  // Validate input
  THYRA_ASSERT_LINEAR_OP_VEC_APPLY_SPACES("sillyCgSolve()", A, Thyra::NOTRANS, *x, &b);
  Teuchos::EVerbosityLevel vl = Teuchos::VERB_MEDIUM;

  out << "\nStarting CG solver ...\n" << std::scientific << "\ndescribe A:\n"<<describe(A, vl)
      << "\ndescribe b:\n"<<describe(b, vl)<<"\ndescribe x:\n"<<describe(*x, vl)<<"\n";

  // Initialization
  const RCP<const VectorSpaceBase<Scalar> > space = A.domain();
  const RCP<VectorBase<Scalar> > r = createMember(space);
  // r = -A*x + b
  V_V(r.ptr(), b); apply<Scalar>(A, NOTRANS, *x, r.ptr(), -one, one);
  const ScalarMag r0_nrm = norm(*r);
  if (r0_nrm==zero) return true;
  const RCP<VectorBase<Scalar> > p = createMember(space), q = createMember(space);
  Scalar rho_old = -one;

  // Perform the iterations
  for( int iter = 0; iter <= maxNumIters; ++iter ) {

    // Check convergence and output iteration
    const ScalarMag r_nrm = norm(*r);
    const bool isConverged = r_nrm/r0_nrm <= tolerance;
    if( iter%(maxNumIters/10+1) == 0 || iter == maxNumIters || isConverged ) {
      out << "Iter = " << iter << ", ||b-A*x||/||b-A*x0|| = " << (r_nrm/r0_nrm) << std::endl;
      if( r_nrm/r0_nrm < tolerance ) return true; // Success!
    }

    // Compute iteration
    const Scalar rho = inner(*r, *r);        // <r,r>              -> rho
    if (iter==0) V_V(p.ptr(), *r);           // r                  -> p   (iter == 0)
    else Vp_V( p.ptr(), *r, rho/rho_old );   // r+(rho/rho_old)*p  -> p   (iter  > 0)
    apply<Scalar>(A, NOTRANS, *p, q.ptr());  // A*p                -> q
    const Scalar alpha = rho/inner(*p, *q);  // rho/<p,q>          -> alpha
    Vp_StV( x, +alpha, *p );                 // +alpha*p + x       -> x
    Vp_StV( r.ptr(), -alpha, *q );           // -alpha*q + r       -> r
    rho_old = rho;                           // rho                -> rho_old (for next iter)

  }

  return false; // Failure

} // end sillyCgSolve

void ImplicitBDFStepperRampingStepControl<Scalar>::initialize(
  const StepperBase<Scalar>& stepper)
{
  // Initialize can be called from the stepper when setInitialCondition
  // is called.
  using Teuchos::as;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Thyra::createMember;

  // Set initial time:
  TimeRange<Scalar> stepperRange = stepper.getTimeRange();
  TEUCHOS_TEST_FOR_EXCEPTION(
      !stepperRange.isValid(),
      std::logic_error,
      "Error, Stepper does not have valid time range for initialization "
      "of ImplicitBDFStepperRampingStepControl!\n");

  if (is_null(parameterList_)) {
    RCP<Teuchos::ParameterList> emptyParameterList =
      Teuchos::rcp(new Teuchos::ParameterList);
    this->setParameterList(emptyParameterList);
  }

  if (is_null(errWtVecCalc_)) {
    RCP<ImplicitBDFStepperErrWtVecCalc<Scalar> > IBDFErrWtVecCalc =
      rcp(new ImplicitBDFStepperErrWtVecCalc<Scalar>());
    errWtVecCalc_ = IBDFErrWtVecCalc;
  }

  stepControlState_ = UNINITIALIZED;

  requestedStepSize_ = Scalar(-1.0);
  currentStepSize_ = initialStepSize_;
  currentOrder_ = 1;
  nextStepSize_ = initialStepSize_;
  nextOrder_ = 1;
  numberOfSteps_ = 0;
  totalNumberOfFailedSteps_ = 0;
  countOfConstantStepsAfterFailure_ = 0;

  if (is_null(delta_)) {
    delta_ = createMember(stepper.get_x_space());
  }
  if (is_null(errWtVec_)) {
    errWtVec_ = createMember(stepper.get_x_space());
  }
  V_S(delta_.ptr(),ST::zero());

  if ( doOutput_(Teuchos::VERB_HIGH) ) {
    RCP<Teuchos::FancyOStream> out = this->getOStream();
    Teuchos::OSTab ostab(out,1,"initialize");
    *out << "currentOrder_ = " << currentOrder_ << std::endl;
    *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
  }

  setStepControlState_(BEFORE_FIRST_STEP);

}

const Teuchos::RCP<TriKota::DiagonalROME<Scalar> >
createModel(
    const int globalDim,
    const typename Teuchos::ScalarTraits<Scalar>::magnitudeType &g_offset
)
{
    using Teuchos::RCP;

    const RCP<const Teuchos::Comm<Thyra::Ordinal> > comm =
        Teuchos::DefaultComm<Thyra::Ordinal>::getComm();

    const int numProcs = comm->getSize();
    TEUCHOS_TEST_FOR_EXCEPT_MSG( numProcs > globalDim,
                                 "Error, the number of processors can not be greater than the global"
                                 " dimension of the vectors!." );
    const int localDim = globalDim / numProcs;
    const int localDimRemainder = globalDim % numProcs;
    TEUCHOS_TEST_FOR_EXCEPT_MSG( localDimRemainder != 0,
                                 "Error, the number of processors must divide into the global number"
                                 " of elements exactly for now!." );

    const RCP<TriKota::DiagonalROME<Scalar> > model =
        Teuchos::rcp(new TriKota::DiagonalROME<Scalar>(localDim));
    const RCP<const Thyra::VectorSpaceBase<Scalar> > p_space = model->get_p_space(0);
    const RCP<Thyra::VectorBase<Scalar> > ps = createMember(p_space);
    const Scalar ps_val = 2.0;
    Thyra::V_S(ps.ptr(), ps_val);
    model->setSolutionVector(ps);
    model->setScalarOffset(g_offset);

    return model;
}

void DefaultPolyLineSearchPointEvaluator<Scalar>::computePoint( const ScalarMag &alpha,
  const Ptr<Thyra::VectorBase<Scalar> > &p
  ) const
{
  typedef ScalarTraits<Scalar> ST;
  using Teuchos::as;
  using Thyra::V_V;
  using Thyra::Vp_StV;
  V_V( p, *vecs_[0] );
  if (alpha != ST::zero()) {
    ScalarMag alpha_i = alpha;
    const int n = vecs_.size();
    for (int i = 1; i < n; ++i, alpha_i *= alpha) {
      Vp_StV(p, alpha_i, *vecs_[i]); 
    }
  }
}

void ExplicitModelEvaluator<Scalar>::
buildInverseMassMatrix() const
{
  typedef Thyra::ModelEvaluatorBase MEB;
  using Teuchos::RCP;
  using Thyra::createMember;
  
  RCP<const Thyra::ModelEvaluator<Scalar> > me = this->getUnderlyingModel();

  // first allocate space for the mass matrix
  RCP<Thyra::LinearOpBase<Scalar> > mass = me->create_W_op();

  // intialize a zero to get rid of the x-dot 
  if(zero_==Teuchos::null) {
    zero_ = Thyra::createMember(*me->get_x_space());
    Thyra::assign(zero_.ptr(),0.0);
  }
  
  // request only the mass matrix from the physics
  // Model evaluator builds: alpha*u_dot + beta*F(u) = 0
  MEB::InArgs<Scalar>  inArgs  = me->createInArgs();
  inArgs.set_x(createMember(me->get_x_space()));
  inArgs.set_x_dot(zero_);
  inArgs.set_alpha(-1.0);
  inArgs.set_beta(0.0);

  // set the one time beta to ensure dirichlet conditions
  // are correctly included in the mass matrix: do it for
  // both epetra and Tpetra. If a panzer model evaluator has
  // not been passed in...oh well you get what you asked for!
  if(panzerModel_!=Teuchos::null)
    panzerModel_->setOneTimeDirichletBeta(-1.0);
  else if(panzerEpetraModel_!=Teuchos::null)
    panzerEpetraModel_->setOneTimeDirichletBeta(-1.0);

  // set only the mass matrix
  MEB::OutArgs<Scalar> outArgs = me->createOutArgs();
  outArgs.set_W_op(mass);

  // this will fill the mass matrix operator 
  me->evalModel(inArgs,outArgs);

  if(!massLumping_) {
    invMassMatrix_ = Thyra::inverse<Scalar>(*me->get_W_factory(),mass);
  }
  else {
    // build lumped mass matrix (assumes all positive mass entries, does a simple sum)
    Teuchos::RCP<Thyra::VectorBase<Scalar> > ones = Thyra::createMember(*mass->domain());
    Thyra::assign(ones.ptr(),1.0);

    RCP<Thyra::VectorBase<Scalar> > invLumpMass = Thyra::createMember(*mass->range());
    Thyra::apply(*mass,Thyra::NOTRANS,*ones,invLumpMass.ptr());
    Thyra::reciprocal(*invLumpMass,invLumpMass.ptr());

    invMassMatrix_ = Thyra::diagonal(invLumpMass);
  }
}

DiagonalROME<Scalar>::DiagonalROME(
    const int localDim,
    const RCP<const Teuchos::Comm<Thyra::Ordinal> > &comm
)
    :Np_(1), Ng_(1), comm_(comm), localDim_(localDim),
     nonlinearTermFactor_(0.0), g_offset_(0.0)
{

    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Thyra::createMember;

    TEUCHOS_ASSERT( localDim > 0 );

    // Get the comm
    if (is_null(comm_)) {
        comm_ = Teuchos::DefaultComm<Thyra::Ordinal>::getComm();
    }

    // Locally replicated space for g
    g_space_ = Thyra::locallyReplicatedDefaultSpmdVectorSpace<Scalar>(comm_, 1);

    // Distributed space for p
    p_space_ = Thyra::defaultSpmdVectorSpace<Scalar>(comm_, localDim, -1);

    // Default solution
    const RCP<Thyra::VectorBase<Scalar> > ps = createMember<Scalar>(p_space_);
    V_S(ps.ptr(), ST::zero());
    ps_ = ps;

    // Default diagonal
    const RCP<Thyra::VectorBase<Scalar> > diag = createMember<Scalar>(p_space_);
    V_S(diag.ptr(), ST::one());
    diag_ = diag;
    diag_bar_ = diag;

    // Default inner product
    const RCP<Thyra::VectorBase<Scalar> > s_bar = createMember<Scalar>(p_space_);
    V_S(s_bar.ptr(), ST::one());
    s_bar_ = s_bar;

    // Default response offset
    g_offset_ = ST::zero();

}

bool DefaultBlockedTriangularLinearOpWithSolve<Scalar>::opSupportedImpl(
  EOpTransp M_trans
  ) const
{
  using Thyra::opSupported;
  assertBlockFillIsActive(false);
  for ( int k = 0; k < numDiagBlocks_; ++k ) {
    if ( !opSupported(*diagonalBlocks_[k].getConstObj(),M_trans) )
      return false;
  }
  return true;
  // ToDo: To be safe we really should do a collective reduction of all
  // clusters of processes.  However, for the typical use case, every block
  // will return the same info and we should be safe!
}

Simple2DModelEvaluator<Scalar>::Simple2DModelEvaluator()
  : x_space_(Thyra::defaultSpmdVectorSpace<Scalar>(2)),
    f_space_(x_space_),
    W_factory_(Thyra::defaultSerialDenseLinearOpWithSolveFactory<Scalar>()),
    d_(0.0),
    p_(x_space_->dim(), Teuchos::ScalarTraits<Scalar>::zero()),
    showGetInvalidArg_(false)
{

  using Teuchos::RCP;
  using Thyra::VectorBase;
  using Thyra::createMember;
  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  
  MEB::InArgsSetup<Scalar> inArgs;
  inArgs.setModelEvalDescription(this->description());
  inArgs.setSupports(MEB::IN_ARG_x);
  prototypeInArgs_ = inArgs;
  
  MEB::OutArgsSetup<Scalar> outArgs;
  outArgs.setModelEvalDescription(this->description());
  outArgs.setSupports(MEB::OUT_ARG_f);
  outArgs.setSupports(MEB::OUT_ARG_W_op);
  outArgs.setSupports(MEB::OUT_ARG_W_prec);
  prototypeOutArgs_ = outArgs;

  nominalValues_ = inArgs;
  x0_ = createMember(x_space_);
  V_S(x0_.ptr(), ST::zero());
  nominalValues_.set_x(x0_);

  set_d(10.0);
  set_p(Teuchos::tuple<Scalar>(2.0, 0.0)());
  set_x0(Teuchos::tuple<Scalar>(1.0, 1.0)());

}

NonlinearCGUtils::ESolveReturn
NonlinearCG<Scalar>::doSolve(
  const Ptr<Thyra::VectorBase<Scalar> > &p_inout,
  const Ptr<ScalarMag> &g_opt_out,
  const Ptr<const ScalarMag> &g_reduct_tol_in,
  const Ptr<const ScalarMag> &g_grad_tol_in,
  const Ptr<const ScalarMag> &alpha_init_in,
  const Ptr<int> &numIters_out
  )
{

  typedef ScalarTraits<Scalar> ST;
  typedef ScalarTraits<ScalarMag> SMT;
  
  using Teuchos::null;
  using Teuchos::as;
  using Teuchos::tuple;
  using Teuchos::rcpFromPtr;
  using Teuchos::optInArg;
  using Teuchos::inOutArg;
  using GlobiPack::computeValue;
  using GlobiPack::PointEval1D;
  using Thyra::VectorSpaceBase;
  using Thyra::VectorBase;
  using Thyra::MultiVectorBase;
  using Thyra::scalarProd;
  using Thyra::createMember;
  using Thyra::createMembers;
  using Thyra::get_ele;
  using Thyra::norm;
  using Thyra::V_StV;
  using Thyra::Vt_S;
  using Thyra::eval_g_DgDp;
  typedef Thyra::Ordinal Ordinal;
  typedef Thyra::ModelEvaluatorBase MEB;
  namespace NCGU = NonlinearCGUtils;
  using std::max;

  // Validate input

  g_opt_out.assert_not_null();

  // Set streams

  const RCP<Teuchos::FancyOStream> out = this->getOStream();
  linesearch_->setOStream(out);

  // Determine what step constants will be computed

  const bool compute_beta_PR =
    (
      solverType_ == NCGU::NONLINEAR_CG_PR_PLUS
      ||
      solverType_ == NCGU::NONLINEAR_CG_FR_PR
      );

  const bool compute_beta_HS = (solverType_ == NCGU::NONLINEAR_CG_HS);

  //
  // A) Set up the storage for the algorithm
  //
  
  const RCP<DefaultPolyLineSearchPointEvaluator<Scalar> >
    pointEvaluator = defaultPolyLineSearchPointEvaluator<Scalar>();

  const RCP<UnconstrainedOptMeritFunc1D<Scalar> >
    meritFunc = unconstrainedOptMeritFunc1D<Scalar>(
      model_, paramIndex_, responseIndex_ );

  const RCP<const VectorSpaceBase<Scalar> >
    p_space = model_->get_p_space(paramIndex_),
    g_space = model_->get_g_space(responseIndex_);

  // Stoarge for current iteration
  RCP<VectorBase<Scalar> >
    p_k = rcpFromPtr(p_inout),        // Current solution for p
    p_kp1 = createMember(p_space),    // Trial point for p (in line search)
    g_vec = createMember(g_space),    // Vector (size 1) form of objective g(p) 
    g_grad_k = createMember(p_space), // Gradient of g DgDp^T
    d_k = createMember(p_space),      // Search direction
    g_grad_k_diff_km1 = null;         // g_grad_k - g_grad_km1 (if needed)

  // Storage for previous iteration
  RCP<VectorBase<Scalar> >
    g_grad_km1 = null, // Will allocate if we need it!
    d_km1 = null; // Will allocate if we need it!
  ScalarMag
    alpha_km1 = SMT::zero(),
    g_km1 = SMT::zero(),
    g_grad_km1_inner_g_grad_km1 = SMT::zero(),
    g_grad_km1_inner_d_km1 = SMT::zero();
  
  if (compute_beta_PR || compute_beta_HS) {
    g_grad_km1 = createMember(p_space);
    g_grad_k_diff_km1 = createMember(p_space);
  }
  
  if (compute_beta_HS) {
    d_km1 = createMember(p_space);
  }
//.........这里部分代码省略.........

bool runCgSolveExample(
  const int dim,
  const Scalar diagScale,
  const bool symOp,
  const bool showAllTests,
  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType tolerance,
  const int maxNumIters
  )
{

  using Teuchos::as;
  using Teuchos::null;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::OSTab;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Thyra::multiply;
  using Thyra::scale;
  typedef typename ST::magnitudeType  ScalarMag;
  bool success = true;
  bool result;
  Teuchos::RCP<Teuchos::FancyOStream> out =
    Teuchos::VerboseObjectBase::getDefaultOStream();
  *out << "\n***\n*** Running silly CG solver using scalar type = \'"
       << ST::name() << "\' ...\n***\n";
  Teuchos::Time timer("");
  timer.start(true);

  //
  // (A) Setup a simple linear system with tridiagonal operator:
  //
  //       [   a*2   -1                         ]
  //       [ -r(1)  a*2       -1                ]
  //  A =  [          .        .        .       ]
  //       [             -r(n-2)      a*2    -1 ]
  //       [                      -r(n-1)   a*2 ]
  //

  // (A.1) Create the tridiagonal matrix operator
  *out << "\nConstructing tridiagonal matrix A of dimension = " << dim
       << " and diagonal multiplier = " << diagScale << " ...\n";
  Teuchos::Array<Scalar> lower(dim-1), diag(dim), upper(dim-1);
  const Scalar
    up = -ST::one(),
    diagTerm = as<Scalar>(2.0) * diagScale * ST::one(),
    low = -(symOp ? ST::one() : ST::random());
  int k = 0;
  // First row
  diag[k] = diagTerm; upper[k] = up;
  // Middle rows
  for( k = 1; k < dim - 1; ++k ) {
    lower[k-1] = low; diag[k] = diagTerm; upper[k] = up;
  }
  // Last row
  lower[k-1] = low; diag[k] = diagTerm;
  RCP<const Thyra::LinearOpBase<Scalar> > A =
    rcp(new ExampleTridiagSerialLinearOp<Scalar>(dim, lower, diag, upper));

  // (A.2) Testing the linear operator constructed linear operator
  *out << "\nTesting the constructed linear operator A ...\n";
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.enable_all_tests(false);
  linearOpTester.check_linear_properties(true);
  linearOpTester.set_all_error_tol(tolerance);
  linearOpTester.set_all_warning_tol(1e-2*tolerance);
  linearOpTester.show_all_tests(showAllTests);
  result = linearOpTester.check(*A, out.ptr());
  if(!result) success = false;

  // (A.3) Create RHS vector b and set to a random value
  RCP<Thyra::VectorBase<Scalar> > b = createMember(A->range());
  Thyra::seed_randomize<Scalar>(0);
  Thyra::randomize( -ST::one(), +ST::one(), b.ptr() );

  // (A.4) Create LHS vector x and set to zero
  RCP<Thyra::VectorBase<Scalar> > x = createMember(A->domain());
  Thyra::V_S( x.ptr(), ST::zero() );

  // (A.5) Create the final linear system
  if(!symOp) {
    *out << "\nSetting up normal equations for unsymmetric system A^H*(A*x-b) => new A*x = b ...\n";
    // A^H*A
    RCP<const Thyra::LinearOpBase<Scalar> > AtA = multiply(adjoint(A), A);
    // A^H*b
    RCP<Thyra::VectorBase<Scalar> > nb = createMember(AtA->range());
    Thyra::apply<Scalar>(*A, Thyra::CONJTRANS, *b, nb.ptr());
    A = AtA;
    b = nb;
  }

  // (A.6) Testing the linear operator used with the solve
  *out << "\nTesting the linear operator used with the solve ...\n";
  linearOpTester.check_for_symmetry(true);
  result = linearOpTester.check(*A, out.ptr());
  if(!result) success = false;

  //
  // (B) Solve the linear system with the silly CG solver
  //
  *out << "\nSolving the linear system with sillyCgSolve(...) ...\n";
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{

    using std::endl;
    typedef double Scalar;
    typedef double ScalarMag;
    using Teuchos::describe;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcp_implicit_cast;
    using Teuchos::rcp_dynamic_cast;
    using Teuchos::as;
    using Teuchos::ParameterList;
    using Teuchos::CommandLineProcessor;
    typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
    typedef Thyra::ModelEvaluatorBase MEB;
    using Thyra::createMember;
    using Thyra::createMembers;

    bool success = true;

    Teuchos::GlobalMPISession mpiSession(&argc,&argv);

    RCP<Epetra_Comm> epetra_comm;
#ifdef HAVE_MPI
    epetra_comm = rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
#else
    epetra_comm = rcp( new Epetra_SerialComm );
#endif // HAVE_MPI

    RCP<Teuchos::FancyOStream>
    out = Teuchos::VerboseObjectBase::getDefaultOStream();

    try {

        //
        // A) Read commandline options
        //

        CommandLineProcessor clp;
        clp.throwExceptions(false);
        clp.addOutputSetupOptions(true);

        std::string paramsFileName = "";
        clp.setOption( "params-file", &paramsFileName,
                       "File name for XML parameters" );

        std::string extraParamsString = "";
        clp.setOption( "extra-params", &extraParamsString,
                       "Extra XML parameter string" );

        Teuchos::EVerbosityLevel verbLevel = Teuchos::VERB_DEFAULT;
        setVerbosityLevelOption( "verb-level", &verbLevel,
                                 "Top-level verbosity level.  By default, this gets deincremented as you go deeper into numerical objects.",
                                 &clp );

        double finalTime = 1.0;
        clp.setOption( "final-time", &finalTime, "Final time (the inital time)" );

        int numTimeSteps = 2;
        clp.setOption( "num-time-steps", &numTimeSteps, "Number of time steps" );

        bool dumpFinalSolutions = false;
        clp.setOption(
            "dump-final-solutions", "no-dump-final-solutions", &dumpFinalSolutions,
            "Determine if the final solutions are dumpped or not." );

        double maxStateError = 1e-6;
        clp.setOption( "max-state-error", &maxStateError,
                       "The maximum allowed error in the integrated state in relation to the exact state solution" );

        // ToDo: Read in more parameters

        CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
        if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;

        if ( Teuchos::VERB_DEFAULT == verbLevel )
            verbLevel = Teuchos::VERB_LOW;

        const Teuchos::EVerbosityLevel
        solnVerbLevel = ( dumpFinalSolutions ? Teuchos::VERB_EXTREME : verbLevel );

        //
        // B) Get the base parameter list that all other parameter lists will be
        // read from.
        //

        RCP<ParameterList> paramList = Teuchos::parameterList();
        if (paramsFileName.length())
            updateParametersFromXmlFile( paramsFileName, &*paramList );
        if (extraParamsString.length())
            updateParametersFromXmlString( extraParamsString, &*paramList );

        paramList->validateParameters(*getValidParameters());

        //
        // C) Create the Stratimikos linear solver factories.
        //

        // Get the linear solve strategy that will be used to solve for the linear
//.........这里部分代码省略.........

bool run_composite_linear_ops_tests(
  const Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> > comm,
  const int n,
  const bool useSpmd,
  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType &tol,
  const bool dumpAll,
  Teuchos::FancyOStream *out_arg
  )
{

  using Teuchos::as;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType    ScalarMag;
  typedef Teuchos::ScalarTraits<ScalarMag> STM;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::null;
  using Teuchos::rcp_const_cast;
  using Teuchos::rcp_dynamic_cast;
  using Teuchos::dyn_cast;
  using Teuchos::OSTab;
  using Thyra::relErr;
  using Thyra::passfail;

  RCP<Teuchos::FancyOStream>
    out = rcp(new Teuchos::FancyOStream(rcp(out_arg,false)));

  const Teuchos::EVerbosityLevel
    verbLevel = dumpAll?Teuchos::VERB_EXTREME:Teuchos::VERB_HIGH;

  if (nonnull(out)) *out
    << "\n*** Entering run_composite_linear_ops_tests<"<<ST::name()<<">(...) ...\n";

  bool success = true, result;

  const ScalarMag warning_tol = ScalarMag(1e-2)*tol, error_tol = tol;
  Thyra::LinearOpTester<Scalar> linearOpTester;
  linearOpTester.linear_properties_warning_tol(warning_tol);
  linearOpTester.linear_properties_error_tol(error_tol);
  linearOpTester.adjoint_warning_tol(warning_tol);
  linearOpTester.adjoint_error_tol(error_tol);
  linearOpTester.dump_all(dumpAll);
  Thyra::LinearOpTester<Scalar> symLinearOpTester(linearOpTester);
  symLinearOpTester.check_for_symmetry(true);
  symLinearOpTester.symmetry_warning_tol(STM::squareroot(warning_tol));
  symLinearOpTester.symmetry_error_tol(STM::squareroot(error_tol));

  RCP<const Thyra::VectorSpaceBase<Scalar> > space;
  if(useSpmd) space = Thyra::defaultSpmdVectorSpace<Scalar>(comm,n,-1);
  else space = Thyra::defaultSpmdVectorSpace<Scalar>(n);
  if (nonnull(out)) *out
    << "\nUsing a basic vector space described as " << describe(*space,verbLevel) << " ...\n";

  if (nonnull(out)) *out << "\nCreating random n x (n/2) multi-vector origA ...\n";
  RCP<Thyra::MultiVectorBase<Scalar> >
    mvOrigA = createMembers(space,n/2,"origA");
  Thyra::seed_randomize<Scalar>(0);
  //RTOpPack::show_spmd_apply_op_dump = true;
  Thyra::randomize( as<Scalar>(as<Scalar>(-1)*ST::one()), as<Scalar>(as<Scalar>(+1)*ST::one()),
    mvOrigA.ptr() );
  RCP<const Thyra::LinearOpBase<Scalar> >
    origA = mvOrigA;
  if (nonnull(out)) *out << "\norigA =\n" << describe(*origA,verbLevel);
  //RTOpPack::show_spmd_apply_op_dump = false;

  if (nonnull(out)) *out << "\nTesting origA ...\n";
  Thyra::seed_randomize<Scalar>(0);
  result = linearOpTester.check(*origA, out.ptr());
  if(!result) success = false;

  if (nonnull(out)) *out
    << "\nCreating implicit scaled linear operator A1 = scale(0.5,origA) ...\n";
  RCP<const Thyra::LinearOpBase<Scalar> >
    A1 = scale(as<Scalar>(0.5),origA);
  if (nonnull(out)) *out << "\nA1 =\n" << describe(*A1,verbLevel);

  if (nonnull(out)) *out << "\nTesting A1 ...\n";
  Thyra::seed_randomize<Scalar>(0);
  result = linearOpTester.check(*A1,out.ptr());
  if(!result) success = false;

  if (nonnull(out)) *out << "\nTesting that A1.getOp() == origA ...\n";
  Thyra::seed_randomize<Scalar>(0);
  result = linearOpTester.compare(
    *dyn_cast<const Thyra::DefaultScaledAdjointLinearOp<Scalar> >(*A1).getOp(),
    *origA,out.ptr());
  if(!result) success = false;

  {

    if (nonnull(out)) *out
      << "\nUnwrapping origA to get non-persisting pointer to origA_1, scalar and transp ...\n";
    Scalar  scalar;
    Thyra::EOpTransp transp;
    const Thyra::LinearOpBase<Scalar> *origA_1 = NULL;
    unwrap( *origA, &scalar, &transp, &origA_1 );
    TEUCHOS_TEST_FOR_EXCEPT( origA_1 == NULL );

    if (nonnull(out)) *out << "\nscalar = " << scalar << " == 1 ? ";
    result = (scalar == ST::one());
//.........这里部分代码省略.........

void CubicSplineInterpolator<Scalar>::interpolate(
  const Array<Scalar> &t_values,
  typename DataStore<Scalar>::DataStoreVector_t *data_out
  ) const
{
  using Teuchos::as;
  using Teuchos::outArg;
  typedef Teuchos::ScalarTraits<Scalar> ST;

  TEUCHOS_TEST_FOR_EXCEPTION( nodesSet_ == false, std::logic_error,
      "Error!, setNodes must be called before interpolate"
      );
#ifdef HAVE_RYTHMOS_DEBUG
  // Check that our nodes_ have not changed between the call to setNodes and interpolate
  assertNodesUnChanged<Scalar>(*nodes_,*nodes_copy_);
  // Assert that the base interpolator preconditions are satisfied
  assertBaseInterpolatePreconditions(*nodes_,t_values,data_out);
#endif // HAVE_RYTHMOS_DEBUG
  
  // Output info
  const RCP<FancyOStream> out = this->getOStream();
  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel(); 
  Teuchos::OSTab ostab(out,1,"CSI::interpolator");
  if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
    *out << "nodes_:" << std::endl;
    for (Teuchos::Ordinal i=0 ; i<(*nodes_).size() ; ++i) {
      *out << "nodes_[" << i << "] = " << std::endl;
      (*nodes_)[i].describe(*out,Teuchos::VERB_EXTREME);
    }
    *out << "t_values = " << std::endl;
    for (Teuchos::Ordinal i=0 ; i<t_values.size() ; ++i) {
      *out << "t_values[" << i << "] = " << t_values[i] << std::endl;
    }
  }

  data_out->clear();

  // Return immediately if no time points are requested ...
  if (t_values.size() == 0) {
    return;
  }

  if ((*nodes_).size() == 1) {
    // trivial case of one node.  Preconditions assert that t_values[0] ==
    // (*nodes_)[0].time so we can just pass it out
    DataStore<Scalar> DS((*nodes_)[0]);
    data_out->push_back(DS);
  }
  else { // (*nodes_).size() >= 2
    int n = 0; // index into t_values
    // Loop through all of the time interpolation points in the buffer and
    // satisfiy all of the requested time points that you find.  NOTE: The
    // loop will be existed once all of the time points are satisified (see
    // return below).
    for (Teuchos::Ordinal i=0 ; i < (*nodes_).size()-1; ++i) {
      const Scalar& ti = (*nodes_)[i].time;
      const Scalar& tip1 = (*nodes_)[i+1].time;
      const TimeRange<Scalar> range_i(ti,tip1);
      // For the interpolation range of [ti,tip1], satisify all of the
      // requested points in this range.
      while ( range_i.isInRange(t_values[n]) ) {
        // First we check for exact node matches:
        if (compareTimeValues(t_values[n],ti)==0) {
          DataStore<Scalar> DS((*nodes_)[i]);
          data_out->push_back(DS);
        }
        else if (compareTimeValues(t_values[n],tip1)==0) {
          DataStore<Scalar> DS((*nodes_)[i+1]);
          data_out->push_back(DS);
        } else {
          if (!splineCoeffComputed_) {
            computeCubicSplineCoeff<Scalar>(*nodes_,outArg(splineCoeff_));
            splineCoeffComputed_ = true;
          }
          DataStore<Scalar> DS;
          RCP<Thyra::VectorBase<Scalar> > x = createMember((*nodes_)[i].x->space());
          evaluateCubicSpline<Scalar>( splineCoeff_, i, t_values[n], outArg(*x) );
          DS.time = t_values[n];
          DS.x = x;
          DS.accuracy = ST::zero();
          data_out->push_back(DS);
        }
        // Move to the next user time point to consider!
        n++;
        if (n == as<int>(t_values.size())) {
          // WE ARE ALL DONE!  MOVE OUT!
          return;
        }
      }
      // Move on the the next interpolation time range
    }
  } // (*nodes_).size() == 1
}

bool tBlockJacobiPreconditionerFactory::test_initializePrec(int verbosity,std::ostream & os)
{
   using Thyra::zero;

   bool status = false;
   bool allPassed = true;

   std::string constrType[3] = {
       std::string("Static"),
       std::string("2x2 Static Strategy"),
       std::string("3x3 Static Strategy") };

   // three by three bloock diagonal 
   std::vector<RCP<const Thyra::LinearOpBase<double> > > invD;
   invD.push_back(invF_); invD.push_back(invC_); invD.push_back(invF_);

   // allocate new linear operator
   const RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blkOp
        = Thyra::defaultBlockedLinearOp<double>();
   blkOp->beginBlockFill(3,3);
   blkOp->setBlock(0,0,F_);  blkOp->setBlock(0,1,Bt_);
   blkOp->setBlock(1,0,B_);  blkOp->setBlock(1,1,C_);  blkOp->setBlock(1,2,B_);
   blkOp->setBlock(2,1,Bt_); blkOp->setBlock(2,2,F_);
   blkOp->endBlockFill();

   const RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > invBlkOp
        = Thyra::defaultBlockedLinearOp<double>();
   invBlkOp->beginBlockFill(3,3);
   invBlkOp->setBlock(0,0,invF_);  
   invBlkOp->setBlock(1,1,invC_);
   invBlkOp->setBlock(2,2,invF_);
   invBlkOp->endBlockFill();

   // build factory array
   RCP<JacobiPreconditionerFactory> fact_array[3] 
         = { rcp(new JacobiPreconditionerFactory(invF_,invC_)),
             rcp(new JacobiPreconditionerFactory(rcp(new StaticInvDiagStrategy(invF_,invC_)))),
             rcp(new JacobiPreconditionerFactory(rcp(new StaticInvDiagStrategy(invD)))) };

   RCP<const Thyra::LinearOpBase<double> > A[3] = { 
       block2x2(F_,Bt_,B_,C_),
       block2x2(F_,Bt_,B_,C_), 
       blkOp };

   // this is what the factory should build
   RCP<const Thyra::LinearOpBase<double> > invA[3] = { 
       block2x2(invF_,zero(Bt_->range(),Bt_->domain()),zero(B_->range(),B_->domain()),invC_),
       block2x2(invF_,zero(Bt_->range(),Bt_->domain()),zero(B_->range(),B_->domain()),invC_),
       invBlkOp };

   // test both constructors
   for(int i=0;i<3;i++) {
      RCP<const Thyra::LinearOpBase<double> > op;

      RCP<Thyra::PreconditionerFactoryBase<double> > fact = fact_array[i];
      RCP<Thyra::PreconditionerBase<double> > prec = fact->createPrec();

      // initialize the preconditioner
      fact->initializePrec(Thyra::defaultLinearOpSource(A[i]), &*prec);

      op = prec->getRightPrecOp();
      TEST_EQUALITY(op,Teuchos::null,
            std::endl << "   tBlockJacobiPreconditionerFactory::test_initializePrec "
            << "using \"" << constrType[i] << "\" constructor " << toString(status) 
            << ": Preconditioner \"getRightPrecOp\" is not null (it should be!)");

      op = prec->getLeftPrecOp();
      TEST_EQUALITY(op,Teuchos::null,
            std::endl << "   tBlockJacobiPreconditionerFactory::test_initializePrec "
            << "using \"" << constrType[i] << "\" constructor " << toString(status) 
            << ": Preconditioner \"getLeftPrecOp\" is not null (it should be!)");

      op = prec->getUnspecifiedPrecOp();
      TEST_NOT_EQUAL(op,Teuchos::null,
            std::endl << "   tBlockJacobiPreconditionerFactory::test_initializePrec "
            << "using \"" << constrType[i] << "\" constructor " << toString(status) 
            << ": Preconditioner \"getUnspecifiedPrecOp\" is null (it should not be!)");

     LinearOpTester<double> tester;
     tester.show_all_tests(true);
     std::stringstream ss;
     Teuchos::FancyOStream fos(rcpFromRef(ss),"      |||");
     const bool result = tester.compare( *invA[i], *op, &fos );
     TEST_ASSERT(result,
            std::endl << "   tBlockJacobiPreconditionerFactory::test_initializePrec "
            << ": Comparing factory generated operator to correct operator");
     if(not result || verbosity>=10) 
        os << ss.str(); 
   }

   return allPassed;
}

int exampleImplicitlyComposedLinearOperators(
  const int n0,
  const int n1,
  const int n2,
  Teuchos::FancyOStream &out,
  const Teuchos::EVerbosityLevel verbLevel,
  typename Teuchos::ScalarTraits<Scalar>::magnitudeType errorTol,
  const bool testAdjoint
  )
{

  // Using and other declarations
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Teuchos::as;
  using Teuchos::RCP;
  using Teuchos::OSTab;
  using Thyra::VectorSpaceBase;
  using Thyra::VectorBase;
  using Thyra::MultiVectorBase;
  using Thyra::LinearOpBase;
  using Thyra::defaultSpmdVectorSpace;
  using Thyra::randomize;
  using Thyra::identity;
  using Thyra::diagonal;
  using Thyra::multiply;
  using Thyra::add;
  using Thyra::subtract;
  using Thyra::scale;
  using Thyra::adjoint;
  using Thyra::block1x2;
  using Thyra::block2x2;
  using Thyra::block2x2;

  out << "\n***"
      << "\n*** Demonstrating building linear operators for scalar type "
      << ST::name()
      << "\n***\n";

  OSTab tab(out);

  //
  // A) Set up the basic objects and other inputs to build the implicitly
  // composed linear operators.
  //
  
  // Create serial vector spaces in this case
  const RCP<const VectorSpaceBase<Scalar> >
    space0 = defaultSpmdVectorSpace<Scalar>(n0),
    space1 = defaultSpmdVectorSpace<Scalar>(n1),
    space2 = defaultSpmdVectorSpace<Scalar>(n2);

  // Create the component linear operators first as multi-vectors
  const RCP<MultiVectorBase<Scalar> >
    mvA = createMembers(space2, n0, "A"),
    mvB = createMembers(space0, n2, "B"),
    mvC = createMembers(space0, n0, "C"),
    mvE = createMembers(space0, n1, "E"),
    mvF = createMembers(space0, n1, "F"),
    mvJ = createMembers(space2, n1, "J"),
    mvK = createMembers(space1, n2, "K"),
    mvL = createMembers(space2, n1, "L"),
    mvN = createMembers(space0, n1, "N"),
    mvP = createMembers(space2, n1, "P"),
    mvQ = createMembers(space0, n2, "Q");

  // Create the vector diagonal for D
  const RCP<VectorBase<Scalar> > d = createMember(space2);

  // Get the constants
  const Scalar
    one = 1.0,
    beta = 2.0,
    gamma = 3.0,
    eta = 4.0;

  // Randomize the values in the Multi-Vector
  randomize( -one, +one, mvA.ptr() );
  randomize( -one, +one, mvB.ptr() );
  randomize( -one, +one, mvC.ptr() );
  randomize( -one, +one, d.ptr() );
  randomize( -one, +one, mvE.ptr() );
  randomize( -one, +one, mvF.ptr() );
  randomize( -one, +one, mvJ.ptr() );
  randomize( -one, +one, mvK.ptr() );
  randomize( -one, +one, mvL.ptr() );
  randomize( -one, +one, mvN.ptr() );
  randomize( -one, +one, mvP.ptr() );
  randomize( -one, +one, mvQ.ptr() );

  // Get the linear operator forms of the basic component linear operators
  const RCP<const LinearOpBase<Scalar> >
    A = mvA,
    B = mvB,
    C = mvC,
    E = mvE,
    F = mvF,
    J = mvJ,
    K = mvK,
    L = mvL,
    N = mvN,
//.........这里部分代码省略.........

void computeCubicSplineCoeff(
    const typename DataStore<Scalar>::DataStoreVector_t & data,
    const Ptr<CubicSplineCoeff<Scalar> > & coeffPtr
    )
{
  using Teuchos::outArg;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  using Thyra::createMember;
  TEUCHOS_TEST_FOR_EXCEPTION( 
      (data.size() < 2), std::logic_error,
      "Error!  A minimum of two data points is required for this cubic spline."
      );
  // time data in the DataStoreVector should be unique and sorted 
  Array<Scalar> t;
  Array<Teuchos::RCP<const Thyra::VectorBase<Scalar> > > x_vec, xdot_vec;
  dataStoreVectorToVector<Scalar>( data, &t, &x_vec, &xdot_vec, NULL );
#ifdef HAVE_RYTHMOS_DEBUG
  assertTimePointsAreSorted<Scalar>( t );
#endif // HAVE_RYTHMOS_DEBUG
  // 11/18/08 tscoffe:  Question:  Should I erase everything in coeffPtr or
  // re-use what I can?  For now, I'll erase and create new each time.
  CubicSplineCoeff<Scalar>& coeff = *coeffPtr;
  // If there are only two points, then we do something special and just create
  // a linear polynomial between the points