C++ teuchos::RCP类代码示例

TEUCHOS_UNIT_TEST(tTpetra_GlbEvalData, filtered_dofs)
{
  using Teuchos::RCP;
  using Teuchos::rcp;

  typedef Tpetra::Vector<double,int,panzer::Ordinal64> Tpetra_Vector;
  typedef Tpetra::Map<int,panzer::Ordinal64> Tpetra_Map;
  typedef Tpetra::Import<int,panzer::Ordinal64> Tpetra_Import;


  Teuchos::RCP<Teuchos::MpiComm<int> > comm = Teuchos::rcp(new Teuchos::MpiComm<int>(MPI_COMM_WORLD));

  // This is required
  TEST_ASSERT(comm->getSize()==2);

  std::vector<panzer::Ordinal64> ghosted(4);
  std::vector<panzer::Ordinal64> owned(2);

  // This is a line with 6 notes (numbered 0-5). The boundaries
  // are removed at 0 and 5.
  if(comm->getRank()==0) {
    owned[0] = 1;
    owned[1] = 2;

    ghosted[0] = 0;
    ghosted[1] = 1;
    ghosted[2] = 2;
    ghosted[3] = 3;
  }
  else {
    owned[0] = 3;
    owned[1] = 4;

    ghosted[0] = 2;
    ghosted[1] = 3;
    ghosted[2] = 4;
    ghosted[3] = 5;
  }

  RCP<const Tpetra_Map> ownedMap = rcp(new Tpetra_Map(-1,owned,0,comm));
  RCP<const Tpetra_Map> ghostedMap = rcp(new Tpetra_Map(-1,ghosted,0,comm));
  RCP<const Tpetra_Import> importer = rcp(new Tpetra_Import(ownedMap,ghostedMap));

  TpetraVector_ReadOnly_GlobalEvaluationData<double,int,panzer::Ordinal64> ged;
 
  std::vector<panzer::Ordinal64> constIndex(1);
 
  // setup filtered values
  constIndex[0] = 0;
  ged.useConstantValues(constIndex,2.0);

  constIndex[0] = 5;
  ged.useConstantValues(constIndex,3.0);

  ged.initialize(importer,ghostedMap,ownedMap);

  TEST_THROW(ged.useConstantValues(constIndex,4.0),std::logic_error);

  {
    RCP<Tpetra_Vector> ownedVecTp = rcp(new Tpetra_Vector(ownedMap));
    auto uv_2d = ownedVecTp->getLocalView<Kokkos::HostSpace> ();
    auto ownedVec = Kokkos::subview (uv_2d, Kokkos::ALL (), 0);

    if(comm->getRank()==0) {
      ownedVec(0) = 3.14;
      ownedVec(1) = 1.82;
    }
    else {
      ownedVec(0) = 2.72;
      ownedVec(1) = 6.23;
    }

    // set the owned vector, assure that const can be used
    ged.setOwnedVector_Tpetra(ownedVecTp.getConst());
  }

  // actually do something...
  ged.initializeData();
  ged.globalToGhost(0);

  // check values making sure that the constants are there
  {
    const Tpetra_Vector & ghostedVecTp = *ged.getGhostedVector_Tpetra();
    auto uv_2d = ghostedVecTp.getLocalView<Kokkos::HostSpace> ();
    auto ghostedVecKv = Kokkos::subview (uv_2d, Kokkos::ALL (), 0);

    if(comm->getRank()==0) {
      TEST_EQUALITY(ghostedVecKv(0),2.0);   // <= Replaced constant value
      TEST_EQUALITY(ghostedVecKv(1),3.14);
      TEST_EQUALITY(ghostedVecKv(2),1.82);
      TEST_EQUALITY(ghostedVecKv(3),2.72);
    }
    else {
      TEST_EQUALITY(ghostedVecKv(0),1.82);
      TEST_EQUALITY(ghostedVecKv(1),2.72);
      TEST_EQUALITY(ghostedVecKv(2),6.23);
      TEST_EQUALITY(ghostedVecKv(3),3.0);   // <= Replaced constant value
    }
  }
}

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

  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
Kokkos::initialize();
  // This little trick lets us print to std::cout only if
  // a (dummy) command-line argument is provided.
  int iprint     = argc - 1;
  Teuchos::RCP<std::ostream> outStream;
  Teuchos::oblackholestream bhs; // outputs nothing
  if (iprint > 0)
    outStream = Teuchos::rcp(&std::cout, false);
  else
    outStream = Teuchos::rcp(&bhs, false);

  // Save the format state of the original std::cout.
  Teuchos::oblackholestream oldFormatState;
  oldFormatState.copyfmt(std::cout);

  *outStream \
    << "===============================================================================\n" \
    << "|                                                                             |\n" \
    << "|                    Unit Test (Basis_HGRAD_TET_C1_FEM)                       |\n" \
    << "|                                                                             |\n" \
    << "|     1) Patch test involving mass and stiffness matrices,                    |\n" \
    << "|        for the Neumann problem on a tetrahedral patch                       |\n" \
    << "|        Omega with boundary Gamma.                                           |\n" \
    << "|                                                                             |\n" \
    << "|        - div (grad u) + u = f  in Omega,  (grad u) . n = g  on Gamma        |\n" \
    << "|                                                                             |\n" \
    << "|  Questions? Contact  Pavel Bochev  ([email protected]),                    |\n" \
    << "|                      Denis Ridzal  ([email protected]),                    |\n" \
    << "|                      Kara Peterson ([email protected]).                    |\n" \
    << "|                                                                             |\n" \
    << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
    << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
    << "|                                                                             |\n" \
    << "===============================================================================\n"\
    << "| TEST 1: Patch test                                                          |\n"\
    << "===============================================================================\n";

  
  int errorFlag = 0;

  outStream -> precision(16);


  try {

    int max_order = 1;                                                                  // max total order of polynomial solution
    DefaultCubatureFactory<double>  cubFactory;                                         // create factory
    shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >());  // create parent cell topology
    shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >());     // create relevant subcell (side) topology
    int cellDim = cell.getDimension();
    int sideDim = side.getDimension();

    // Define array containing points at which the solution is evaluated, on the reference tet.
    int numIntervals = 10;
    int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals + 3))/6;
    FieldContainer<double> interp_points_ref(numInterpPoints, 3);
    int counter = 0;
    for (int k=0; k<=numIntervals; k++) {
      for (int j=0; j<=numIntervals; j++) {
        for (int i=0; i<=numIntervals; i++) {
          if (i+j+k <= numIntervals) {
            interp_points_ref(counter,0) = i*(1.0/numIntervals);
            interp_points_ref(counter,1) = j*(1.0/numIntervals);
            interp_points_ref(counter,2) = k*(1.0/numIntervals);
            counter++;
          }
        }
      }
    }

    /* Definition of parent cell. */
    FieldContainer<double> cell_nodes(1, 4, cellDim);
    // funky tet
    cell_nodes(0, 0, 0) = -1.0;
    cell_nodes(0, 0, 1) = -2.0;
    cell_nodes(0, 0, 2) = 0.0;
    cell_nodes(0, 1, 0) = 6.0;
    cell_nodes(0, 1, 1) = 2.0;
    cell_nodes(0, 1, 2) = 0.0;
    cell_nodes(0, 2, 0) = -5.0;
    cell_nodes(0, 2, 1) = 1.0;
    cell_nodes(0, 2, 2) = 0.0;
    cell_nodes(0, 3, 0) = -4.0;
    cell_nodes(0, 3, 1) = -1.0;
    cell_nodes(0, 3, 2) = 3.0;
    // perturbed reference tet
    /*cell_nodes(0, 0, 0) = 0.1;
    cell_nodes(0, 0, 1) = -0.1;
    cell_nodes(0, 0, 2) = 0.2;
    cell_nodes(0, 1, 0) = 1.2;
    cell_nodes(0, 1, 1) = -0.1;
    cell_nodes(0, 1, 2) = 0.05;
    cell_nodes(0, 2, 0) = 0.0;
    cell_nodes(0, 2, 1) = 0.9;
    cell_nodes(0, 2, 2) = 0.1;
    cell_nodes(0, 3, 0) = 0.1;
    cell_nodes(0, 3, 1) = -0.1;
//.........这里部分代码省略.........

TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(dof_pointfield,value,EvalType)
{
  PHX::KokkosDeviceSession session;

  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcp_dynamic_cast;

  // build a dummy workset
  //////////////////////////////////////////////////////////
  int numCells = 2;

  // build basis values
  Teuchos::RCP<shards::CellTopology> topo
    = Teuchos::rcp(new shards::CellTopology(shards::getCellTopologyData< shards::Quadrilateral<4> >()));

  Teuchos::RCP<panzer::PureBasis> sourceBasis = Teuchos::rcp(new panzer::PureBasis("HGrad",1,numCells,topo));
  Teuchos::RCP<panzer::PureBasis> targetBasis = Teuchos::rcp(new panzer::PureBasis("HGrad",2,numCells,topo));

  Teuchos::RCP<PHX::FieldManager<panzer::Traits> > fm
     = Teuchos::rcp(new PHX::FieldManager<panzer::Traits>); 

  // add in some evaluators
  ///////////////////////////////////////////////////
  {
     // fill the basis with some dummy values
     Teuchos::ParameterList p;
     p.set("Name","Pressure");
     p.set("Basis",sourceBasis);
     RCP<panzer::DummyFieldEvaluator<EvalType,panzer::Traits> > e = 
       rcp(new panzer::DummyFieldEvaluator<EvalType,panzer::Traits>(p));

     fm->registerEvaluator<EvalType>(e);
  }

  // add evaluator under test
  ///////////////////////////////////////////////////

  {
    RCP<panzer::DOF_BasisToBasis<EvalType,panzer::Traits> > e = 
      rcp(new panzer::DOF_BasisToBasis<EvalType,panzer::Traits>("Pressure",*sourceBasis,*targetBasis));
    
    fm->registerEvaluator<EvalType>(e);
    
    Teuchos::RCP<PHX::FieldTag> ft = e->evaluatedFields()[0];
    fm->requireField<EvalType>(*ft);
  }

  Teuchos::RCP<panzer::Workset> workset = Teuchos::rcp(new panzer::Workset);
  workset->num_cells = numCells;

  panzer::Traits::SetupData setupData;
  setupData.worksets_ = rcp(new std::vector<panzer::Workset>);
  setupData.worksets_->push_back(*workset);

  std::vector<PHX::index_size_type> derivative_dimensions;
  derivative_dimensions.push_back(4);
  fm->setKokkosExtendedDataTypeDimensions<panzer::Traits::Jacobian>(derivative_dimensions);
  fm->postRegistrationSetup(setupData);

  //fm->writeGraphvizFile();

  panzer::Traits::PreEvalData preEvalData;

  fm->preEvaluate<EvalType>(preEvalData);
  fm->evaluateFields<EvalType>(*workset);
  fm->postEvaluate<EvalType>(0);

  typedef typename EvalType::ScalarT ScalarT;
  typedef typename PHX::MDFieldTypeTraits<PHX::MDField<ScalarT,Cell,BASIS> >::return_type  ScalarView;

  typename PHX::MDField<ScalarT,Cell,BASIS> s("Pressure",sourceBasis->functional);
  typename PHX::MDField<ScalarT,Cell,BASIS> t("Pressure",targetBasis->functional);

  fm->getFieldData<ScalarT,EvalType>(s);
  fm->getFieldData<ScalarT,EvalType>(t);

  typename Teuchos::ScalarTraits<ScalarT>::magnitudeType tol =
    100.0 * Teuchos::ScalarTraits<ScalarT>::eps();

  TEST_FLOATING_EQUALITY(ScalarT(t(0,0)),ScalarT(s(0,0)),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,1)),ScalarT(s(0,1)),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,2)),ScalarT(s(0,2)),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,3)),ScalarT(s(0,3)),tol);

  TEST_FLOATING_EQUALITY(ScalarT(t(1,0)),ScalarT(s(1,0)),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(1,1)),ScalarT(s(1,1)),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(1,2)),ScalarT(s(1,2)),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(1,3)),ScalarT(s(1,3)),tol);

  TEST_FLOATING_EQUALITY(ScalarT(t(0,4)),ScalarT(1.5),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,5)),ScalarT(2.0),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,6)),ScalarT(1.5),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,7)),ScalarT(1.0),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(0,8)),ScalarT(1.5),tol);

  TEST_FLOATING_EQUALITY(ScalarT(t(1,4)),ScalarT(2.5),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(1,5)),ScalarT(3.0),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(1,6)),ScalarT(2.5),tol);
  TEST_FLOATING_EQUALITY(ScalarT(t(1,7)),ScalarT(2.0),tol);
//.........这里部分代码省略.........

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

  // Initialize MPI
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Create a communicator for Epetra objects
#ifdef HAVE_MPI
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  Epetra_SerialComm Comm;
#endif

  Teuchos::CommandLineProcessor clp( false );

  // Default run-time options that can be changed from the command line
  CouplingSolveMethod method          = MATRIX_FREE   ;
  bool                verbose         = true          ;
  int                 NumGlobalNodes  = 20            ;
  int                 probSizeRatio   = 1             ;
  bool                useMatlab       = false         ;
  bool                doOffBlocks     = false         ;
  std::string         solvType        = "seidel"      ;
  // Coupling parameters
  double              alpha           = 0.50          ;
  double              beta            = 0.40          ;
  // Physical parameters
  double              radiation       = 5.67          ;
  double              initVal         = 0.995         ;
  string              outputDir       = "."           ;
  string              goldDir         = "."           ;


  clp.setOption<CouplingSolveMethod>( "solvemethod", &method, 4, SolveMethodValues, SolveMethodNames, "Selects the coupling method to use");
  clp.setOption( "verbose", "no-verbose", &verbose, "Verbosity on or off." );
  clp.setOption( "n", &NumGlobalNodes, "Number of elements" );
  clp.setOption( "nratio", &probSizeRatio, "Ratio of size of problem 2 to problem 1" );
  clp.setOption( "offblocks", "no-offblocks", &doOffBlocks, "Include off-diagonal blocks in preconditioning matrix" );
  clp.setOption( "matlab", "no-matlab", &useMatlab, "Use Matlab debugging engine" );
  clp.setOption( "solvType", &solvType, "Solve Type.  Valid choices are: jacobi, seidel" );
  clp.setOption( "alpha", &alpha, "Interfacial coupling coefficient, alpha" );
  clp.setOption( "beta", &beta, "Interfacial coupling coefficient, beta" );
  clp.setOption( "radiation", &radiation, "Radiation source term coefficient, R" );
  clp.setOption( "initialval", &initVal, "Initial guess for solution values" );
  clp.setOption( "outputdir", &outputDir, "Directory to output mesh and results into. Default is \"./\"" );
  clp.setOption( "golddir", &goldDir, "Directory to read gold test from. Default is \"./\"" );

  Teuchos::CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);

  if( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) 
    return parse_return;

  outputDir += "/";
  goldDir   += "/";

  // Create and reset the Timer
  Epetra_Time myTimer(Comm);
  double startWallTime = myTimer.WallTime();

  // Get the process ID and the total number of processors
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  NumGlobalNodes++; // convert #elements to #nodes

  // The number of unknowns must be at least equal to the number of processors.
  if (NumGlobalNodes < NumProc) 
  {
    cout << "numGlobalNodes = " << NumGlobalNodes 
	 << " cannot be < number of processors = " << NumProc << endl;
    exit(1);
  }

  // Begin Nonlinear Solver ************************************

  // NOTE: For now these parameters apply to all problems handled by
  // Problem_Manager.  Each problem could be made to have its own
  // parameter list as wwell as its own convergence test(s).

  // Create the top level parameter list
  Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr =
    Teuchos::rcp(new Teuchos::ParameterList);
  Teuchos::ParameterList& nlParams = *(nlParamsPtr.get());

  // Set the nonlinear solver method
  nlParams.set("Nonlinear Solver", "Line Search Based");

  // Set the printing parameters in the "Printing" sublist
  Teuchos::ParameterList& printParams = nlParams.sublist("Printing");
  printParams.set("MyPID", MyPID); 
  printParams.set("Output Precision", 3);
  printParams.set("Output Processor", 0);
  printParams.set("Output Information", 
			NOX::Utils::Warning                  +
			NOX::Utils::OuterIteration           + 
			NOX::Utils::InnerIteration           +
			NOX::Utils::Parameters               + 
			NOX::Utils::Details                  + 
//.........这里部分代码省略.........

 Teuchos::RCP<Vector<Real> > clone() const {
   Real stat = 0.0;
   Teuchos::RCP<Vector<Real> > vec = Teuchos::rcp_dynamic_cast<Vector<Real> >(
     Teuchos::rcp_const_cast<Vector<Real> >(vec_->clone()));
   return Teuchos::rcp( new RiskVector( vec, augmented_, stat ) );  
 }

int main(int argc, char **argv)
{
  try {

    // Basis of dimension 3, order 5
    const int d = 3;
    const int p = 5;
    Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d); 
    for (int i=0; i<d; i++) {
      bases[i] = Teuchos::rcp(new Stokhos::PecosOneDOrthogPolyBasis<int,double>(Teuchos::rcp(new Pecos::HermiteOrthogPolynomial), "Hermite", p));
    }
    Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis = 
      Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));

    // Quadrature method
    Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad = 
        Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));

    // Triple product tensor
    Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
      basis->computeTripleProductTensor(basis->size());

    // Expansion method
    Stokhos::QuadOrthogPolyExpansion<int,double> expn(basis, Cijk, quad);

    // Polynomial expansions
    Stokhos::OrthogPolyApprox<int,double> u(basis), v(basis), w(basis);
    u.term(0,0) = 1.0;
    for (int i=0; i<d; i++) {
      u.term(i,1) = 0.4 / d;
      u.term(i,2) = 0.06 / d;
      u.term(i,3) = 0.002 / d;
    }

    // Compute expansion
    expn.log(v,u);
    expn.times(w,v,v);
    expn.plusEqual(w,1.0);
    expn.divide(v,1.0,w);
    //expn.times(v,u,u);

    // Print u and v
    std::cout << "v = 1.0 / (log(u)^2 + 1):" << std::endl;
    std::cout << "\tu = ";
    u.print(std::cout);
    std::cout << "\tv = ";
    v.print(std::cout);

    // Compute moments
    double mean = v.mean();
    double std_dev = v.standard_deviation();

    // Evaluate PCE and function at a point = 0.25 in each dimension
    Teuchos::Array<double> pt(d); 
    for (int i=0; i<d; i++) 
      pt[i] = 0.25;
    double up = u.evaluate(pt);
    double vp = 1.0/(std::log(up)*std::log(up)+1.0);
    double vp2 = v.evaluate(pt);
    
    // Print results
    std::cout << "\tv mean         = " << mean << std::endl;
    std::cout << "\tv std. dev.    = " << std_dev << std::endl;
    std::cout << "\tv(0.25) (true) = " << vp << std::endl;
    std::cout << "\tv(0.25) (pce)  = " << vp2 << std::endl;
    
    // Check the answer
    if (std::abs(vp - vp2) < 1e-2)
      std::cout << "\nExample Passed!" << std::endl;
  }
  catch (std::exception& e) {
    std::cout << e.what() << std::endl;
  }
}

void CEigenLSS::solve()
{
#ifdef CF3_HAVE_TRILINOS
  Timer timer;
  const Uint nb_rows = size();
  cf3_assert(nb_rows == m_system_matrix.outerSize());

  Epetra_SerialComm comm;
  Epetra_Map map(nb_rows, 0, comm);

  Epetra_Vector ep_rhs(View, map, m_rhs.data());
  Epetra_Vector ep_sol(View, map, m_solution.data());

  // Count non-zeros
  std::vector<int> nnz(nb_rows, 0);
  for(int row=0; row < nb_rows; ++row)
  {
    for(MatrixT::InnerIterator it(m_system_matrix, row); it; ++it)
    {
      ++nnz[row];
    }
    cf3_assert(nnz[row]);
  }

  Epetra_CrsMatrix ep_A(Copy, map, &nnz[0]);
  time_matrix_construction = timer.elapsed(); timer.restart();

  // Fill the matrix
  for(int row=0; row < nb_rows; ++row)
  {
    std::vector<int> indices; indices.reserve(nnz[row]);
    std::vector<Real> values; values.reserve(nnz[row]);
    for(MatrixT::InnerIterator it(m_system_matrix, row); it; ++it)
    {
      indices.push_back(it.col());
      values.push_back(it.value());
    }
    ep_A.InsertGlobalValues(row, nnz[row], &values[0], &indices[0]);
  }

  ep_A.FillComplete();

  time_matrix_fill = timer.elapsed(); timer.restart();

///////////////////////////////////////////////////////////////////////////////////////////////
//BEGIN////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////

  Teuchos::RCP<Epetra_CrsMatrix> epetra_A=Teuchos::rcpFromRef(ep_A);
  Teuchos::RCP<Epetra_Vector>    epetra_x=Teuchos::rcpFromRef(ep_sol);
  Teuchos::RCP<Epetra_Vector>    epetra_b=Teuchos::rcpFromRef(ep_rhs);

  const URI config_uri = options().option("config_file").value<URI>();
  const std::string config_path = config_uri.path();

  Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder(config_path); // the most important in general setup

  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream(); // TODO: decouple from fancyostream to ostream or to C stdout when possible
  typedef Teuchos::ParameterList::PrintOptions PLPrintOptions;
  Teuchos::CommandLineProcessor  clp(false); // false: don't throw exceptions

  Teuchos::RCP<const Thyra::LinearOpBase<double> > A = Thyra::epetraLinearOp( epetra_A );
  Teuchos::RCP<Thyra::VectorBase<double> >         x = Thyra::create_Vector( epetra_x, A->domain() );
  Teuchos::RCP<const Thyra::VectorBase<double> >   b = Thyra::create_Vector( epetra_b, A->range() );

  // r = b - A*x, initial L2 norm
  double nrm_r=0.;
  Real systemResidual=-1.;
  {
    Epetra_Vector epetra_r(*epetra_b);
    Epetra_Vector epetra_A_x(epetra_A->OperatorRangeMap());
    epetra_A->Apply(*epetra_x,epetra_A_x);
    epetra_r.Update(-1.0,epetra_A_x,1.0);
    epetra_r.Norm2(&nrm_r);
  }

  // Reading in the solver parameters from the parameters file and/or from
  // the command line.  This was setup by the command-line options
  // set by the setupCLP(...) function above.
  linearSolverBuilder.readParameters(0); // out.get() if want confirmation about the xml file within trilinos
  Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> > lowsFactory = linearSolverBuilder.createLinearSolveStrategy(""); // create linear solver strategy
  lowsFactory->setVerbLevel(Teuchos::VERB_NONE); // set verbosity

//  // print back default and current settings
//  if (opts->trilinos.dumpDefault!=0) {
//    fflush(stdout); cout << flush;
//    _MMESSAGE_(0,1,"Dumping Trilinos/Stratimikos solver defaults to files: 'trilinos_default.txt' and 'trilinos_default.xml'...\n");
//    fflush(stdout); cout << flush;
//    std::ofstream ofs("./trilinos_default.txt");
//    linearSolverBuilder.getValidParameters()->print(ofs,PLPrintOptions().indent(2).showTypes(true).showDoc(true)); // the last true-false is the deciding about whether printing documentation to option or not
//    ofs.flush();ofs.close();
//    ofs.open("./trilinos_default.xml");
//    Teuchos::writeParameterListToXmlOStream(*linearSolverBuilder.getValidParameters(),ofs);
//    ofs.flush();ofs.close();
//  }
//  if (opts->trilinos.dumpCurrXML!=0) {
//    fflush(stdout); cout << flush;
//    _MMESSAGE_(0,1,"Dumping Trilinos/Stratimikos current settings to file: 'trilinos_current.xml'...\n");
//    fflush(stdout); cout << flush;
//    linearSolverBuilder.writeParamsFile(*lowsFactory,"./trilinos_current.xml");
//.........这里部分代码省略.........

Teuchos::RCP<ROL::Vector<Real> > clone() const {
  return Teuchos::rcp( new ConDualStdVector( Teuchos::rcp(new std::vector<Element>(std_vec_->size())) ) );
}

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

  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
  Teuchos::RCP<const Teuchos::Comm<int> > comm
    = Teuchos::DefaultComm<int>::getComm();

  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
  int iprint = argc - 1;
  bool print = (iprint>0);
  Teuchos::RCP<std::ostream> outStream;
  Teuchos::oblackholestream bhs; // outputs nothing
  if (print)
    outStream = Teuchos::rcp(&std::cout, false);
  else
    outStream = Teuchos::rcp(&bhs, false);
    
  bool print0 = print && !comm->getRank();
  Teuchos::RCP<std::ostream> outStream0;
  if (print0)
    outStream0 = Teuchos::rcp(&std::cout, false); 
  else
    outStream0 = Teuchos::rcp(&bhs, false);

  int errorFlag  = 0;

  // *** Example body.

  try {
    /*************************************************************************/
    /************* INITIALIZE BURGERS FEM CLASS ******************************/
    /*************************************************************************/
    int nx      = 256;   // Set spatial discretization.
    RealT alpha = 1.e-3; // Set penalty parameter.
    RealT nl    = 1.0;   // Nonlinearity parameter (1 = Burgers, 0 = linear).
    RealT cH1   = 1.0;   // Scale for derivative term in H1 norm.
    RealT cL2   = 0.0;   // Scale for mass term in H1 norm.
    Teuchos::RCP<BurgersFEM<RealT> > fem
      = Teuchos::rcp(new BurgersFEM<RealT>(nx,nl,cH1,cL2));
    fem->test_inverse_mass(*outStream0);
    fem->test_inverse_H1(*outStream0);
    /*************************************************************************/
    /************* INITIALIZE SIMOPT OBJECTIVE FUNCTION **********************/
    /*************************************************************************/
    Teuchos::RCP<std::vector<RealT> > ud_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx, 1.0) );
    Teuchos::RCP<ROL::Vector<RealT> > ud
      = Teuchos::rcp(new L2VectorPrimal<RealT>(ud_rcp,fem));
    Teuchos::RCP<ROL::Objective_SimOpt<RealT> > pobj
      = Teuchos::rcp(new Objective_BurgersControl<RealT>(fem,ud,alpha));
    /*************************************************************************/
    /************* INITIALIZE SIMOPT EQUALITY CONSTRAINT *********************/
    /*************************************************************************/
    bool hess = true;
    Teuchos::RCP<ROL::EqualityConstraint_SimOpt<RealT> > pcon
      = Teuchos::rcp(new EqualityConstraint_BurgersControl<RealT>(fem,hess));
    /*************************************************************************/
    /************* INITIALIZE VECTOR STORAGE *********************************/
    /*************************************************************************/
    // INITIALIZE CONTROL VECTORS
    Teuchos::RCP<std::vector<RealT> > z_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx+2, 1.0) );
    Teuchos::RCP<std::vector<RealT> > gz_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx+2, 1.0) );
    Teuchos::RCP<std::vector<RealT> > yz_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx+2, 1.0) );
    for (int i=0; i<nx+2; i++) {
      (*yz_rcp)[i] = 2.0*random<RealT>(comm)-1.0;
    }
    Teuchos::RCP<ROL::Vector<RealT> > zp
      = Teuchos::rcp(new PrimalControlVector(z_rcp,fem));
    Teuchos::RCP<ROL::Vector<RealT> > gzp
      = Teuchos::rcp(new DualControlVector(gz_rcp,fem));
    Teuchos::RCP<ROL::Vector<RealT> > yzp
      = Teuchos::rcp(new PrimalControlVector(yz_rcp,fem));
    std::vector<RealT> zvar(1,0.0*random<RealT>(comm));
    std::vector<RealT> gvar(1,random<RealT>(comm));
    std::vector<RealT> yvar(1,random<RealT>(comm));
    ROL::RiskVector<RealT> z(zp,zvar,true), g(gzp,gvar,true), y(yzp,yvar,true);
    // INITIALIZE STATE VECTORS
    Teuchos::RCP<std::vector<RealT> > u_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx, 1.0) );
    Teuchos::RCP<std::vector<RealT> > gu_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx, 1.0) );
    Teuchos::RCP<ROL::Vector<RealT> > up
      = Teuchos::rcp(new PrimalStateVector(u_rcp,fem));
    Teuchos::RCP<ROL::Vector<RealT> > gup
      = Teuchos::rcp(new DualStateVector(gu_rcp,fem));
    // INITIALIZE CONSTRAINT VECTORS
    Teuchos::RCP<std::vector<RealT> > c_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx, 1.0) );
    Teuchos::RCP<std::vector<RealT> > l_rcp
      = Teuchos::rcp( new std::vector<RealT> (nx, 1.0) );
    for (int i=0; i<nx; i++) {
      (*l_rcp)[i] = random<RealT>(comm);
    }
    Teuchos::RCP<ROL::Vector<RealT> > cp
      = Teuchos::rcp(new PrimalConstraintVector(c_rcp,fem));
    Teuchos::RCP<ROL::Vector<RealT> > lp
      = Teuchos::rcp(new DualConstraintVector(l_rcp,fem));
    /*************************************************************************/
//.........这里部分代码省略.........

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

#ifdef HAVE_MPI
  // Initialize MPI
  //
  MPI_Init(&argc,&argv);
#endif

  bool success = false;
  try {
    // Create an Epetra communicator
    //
#ifdef HAVE_MPI
    Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
    Epetra_SerialComm Comm;
#endif

    // Create an Anasazi output manager
    //
    BasicOutputManager<double> printer;
    printer.stream(Errors) << Anasazi_Version() << std::endl << std::endl;

    // Get the sorting std::string from the command line
    //
    std::string which("LM");
    Teuchos::CommandLineProcessor cmdp(false,true);
    cmdp.setOption("sort",&which,"Targetted eigenvalues (SM or LM).");
    if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
      throw -1;
    }


    typedef Epetra_MultiVector MV;
    typedef Epetra_Operator OP;
    typedef MultiVecTraits<double, Epetra_MultiVector> MVT;

    // Number of dimension of the domain
    const int space_dim = 2;

    // Size of each of the dimensions of the domain
    std::vector<double> brick_dim( space_dim );
    brick_dim[0] = 1.0;
    brick_dim[1] = 1.0;

    // Number of elements in each of the dimensions of the domain
    std::vector<int> elements( space_dim );
    elements[0] = 10;
    elements[1] = 10;

    // Create problem
    Teuchos::RCP<ModalProblem> testCase =
      Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );

    // Get the stiffness and mass matrices
    Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
    Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );

    // Eigensolver parameters
    int nev = 10;
    int blockSize = 5;
    int maxIters = 500;
    double tol = 1.0e-8;

    Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(K->OperatorDomainMap(), blockSize) );
    ivec->Random();

    // Create the eigenproblem.
    Teuchos::RCP<BasicEigenproblem<double, MV, OP> > MyProblem =
      Teuchos::rcp( new BasicEigenproblem<double, MV, OP>(K, M, ivec) );

    // Inform the eigenproblem that the operator A is symmetric
    MyProblem->setHermitian(true);

    // Set the number of eigenvalues requested
    MyProblem->setNEV( nev );

    // Inform the eigenproblem that you are finishing passing it information
    bool boolret = MyProblem->setProblem();
    if (boolret != true) {
      printer.print(Errors,"Anasazi::BasicEigenproblem::setProblem() returned an error.\n");
      throw -1;
    }

    // Create parameter list to pass into the solver manager
    //
    Teuchos::ParameterList MyPL;
    MyPL.set( "Which", which );
    MyPL.set( "Block Size", blockSize );
    MyPL.set( "Maximum Iterations", maxIters );
    MyPL.set( "Convergence Tolerance", tol );
    MyPL.set( "Full Ortho", true );
    MyPL.set( "Use Locking", true );
    //
    // Create the solver manager
    LOBPCGSolMgr<double, MV, OP> MySolverMan(MyProblem, MyPL);

    // Solve the problem
    //
    ReturnType returnCode = MySolverMan.solve();
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  Teuchos::GlobalMPISession mpiSession(&argc, &argv);

  bool success = false;
  bool verbose = false;
  try {
    // Parse the command line
    using Teuchos::CommandLineProcessor;
    CommandLineProcessor  clp;
    clp.throwExceptions(false);
    clp.addOutputSetupOptions(true);
    clp.setOption( "v", "disable-verbosity", &verbose, "Enable verbosity" );

    CommandLineProcessor::EParseCommandLineReturn
      parse_return = clp.parse(argc,argv,&std::cerr);

    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
      return parse_return;

    if (verbose)
      std::cout << "Verbosity Activated" << std::endl;
    else
      std::cout << "Verbosity Disabled" << std::endl;

    // Create a communicator for Epetra objects
#ifdef HAVE_MPI
    Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
    Epetra_SerialComm Comm;
#endif

    const int num_elements = 400;

    // Check we have only one processor since this problem doesn't work
    // for more than one proc
    if (Comm.NumProc() > num_elements) {
      std::cerr << "Error! Number of elements must be greate than number of processors!"
        << std::endl;
      return EXIT_FAILURE;
    }

    // Create the model evaluator object
    double paramC = 0.99;
    Teuchos::RCP<ModelEvaluatorHeq<double> > model =
      modelEvaluatorHeq<double>(Teuchos::rcp(&Comm,false),num_elements,paramC);

    ::Stratimikos::DefaultLinearSolverBuilder builder;

    Teuchos::RCP<Teuchos::ParameterList> p =
      Teuchos::rcp(new Teuchos::ParameterList);
    p->set("Linear Solver Type", "AztecOO");
    p->sublist("Linear Solver Types").sublist("AztecOO").sublist("Forward Solve").sublist("AztecOO Settings").set("Output Frequency",20);
    p->set("Preconditioner Type", "Ifpack");
    builder.setParameterList(p);

    Teuchos::RCP< ::Thyra::LinearOpWithSolveFactoryBase<double> >
      lowsFactory = builder.createLinearSolveStrategy("");

    model->set_W_factory(lowsFactory);

    // Create the initial guess
    Teuchos::RCP< ::Thyra::VectorBase<double> >
      initial_guess = model->getNominalValues().get_x()->clone_v();

    Thyra::V_S(initial_guess.ptr(),Teuchos::ScalarTraits<double>::one());

    Teuchos::RCP<NOX::Thyra::Group> nox_group =
      Teuchos::rcp(new NOX::Thyra::Group(*initial_guess, model));
    //Teuchos::rcp(new NOX::Thyra::Group(*initial_guess, model, model->create_W_op(), lowsFactory, Teuchos::null, Teuchos::null));

    nox_group->computeF();

    // Create the NOX status tests and the solver
    // Create the convergence tests
    Teuchos::RCP<NOX::StatusTest::NormF> absresid =
      Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-8));
    Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
      Teuchos::rcp(new NOX::StatusTest::NormWRMS(1.0e-2, 1.0e-8));
    Teuchos::RCP<NOX::StatusTest::Combo> converged =
      Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
    converged->addStatusTest(absresid);
    converged->addStatusTest(wrms);
    Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
      Teuchos::rcp(new NOX::StatusTest::MaxIters(20));
    Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
      Teuchos::rcp(new NOX::StatusTest::FiniteValue);
    Teuchos::RCP<NOX::StatusTest::Combo> combo =
      Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
    combo->addStatusTest(fv);
    combo->addStatusTest(converged);
    combo->addStatusTest(maxiters);

    // Create nox parameter list
    Teuchos::RCP<Teuchos::ParameterList> nl_params =
      Teuchos::rcp(new Teuchos::ParameterList);
    nl_params->set("Nonlinear Solver", "Anderson Accelerated Fixed-Point");
    nl_params->sublist("Anderson Parameters").set("Storage Depth", 5);
    nl_params->sublist("Anderson Parameters").set("Mixing Parameter", 1.0);
    nl_params->sublist("Anderson Parameters").set("Acceleration Start Iteration", 5);
//.........这里部分代码省略.........

  // implementation of "toThyra" for full specialization on SC=double, LO=GO=int and NO=EpetraNode
  // We need the specialization in the cpp file due to a circle dependency in the .hpp files for BlockedCrsMatrix
  Teuchos::RCP<Thyra::LinearOpBase<double> >
  ThyraUtils<double, int, int, EpetraNode>::toThyra(const Teuchos::RCP<Xpetra::BlockedCrsMatrix<double, int, int, EpetraNode> >& mat) {

    int nRows = mat->Rows();
    int nCols = mat->Cols();

    Teuchos::RCP<Xpetra::Matrix<double, int, int, EpetraNode> > Ablock = mat->getMatrix(0,0);
    Teuchos::RCP<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode> > Ablock_wrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<double, int, int, EpetraNode>>(Ablock);
    TEUCHOS_TEST_FOR_EXCEPT(Ablock_wrap.is_null() == true);

    bool bTpetra = false;
    bool bEpetra = false;
#ifdef HAVE_XPETRA_TPETRA
    // Note: Epetra is enabled
#if ((defined(EPETRA_HAVE_OMP)  && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)) || \
     (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT) && defined(HAVE_TPETRA_INST_DOUBLE)))
    Teuchos::RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpetraMat = Teuchos::rcp_dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Ablock_wrap->getCrsMatrix());
    if(tpetraMat!=Teuchos::null) bTpetra = true;
#else
    bTpetra = false;
#endif
#endif

#ifdef HAVE_XPETRA_EPETRA
    Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > epetraMat = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(Ablock_wrap->getCrsMatrix());
    if(epetraMat!=Teuchos::null) bEpetra = true;
#endif

    TEUCHOS_TEST_FOR_EXCEPT(bTpetra == bEpetra); // we only allow Epetra OR Tpetra

    // create new Thyra blocked operator
    Teuchos::RCP<Thyra::PhysicallyBlockedLinearOpBase<Scalar> > blockMat =
        Thyra::defaultBlockedLinearOp<Scalar>();

    blockMat->beginBlockFill(nRows,nCols);

    for (int r=0; r<nRows; ++r) {
      for (int c=0; c<nCols; ++c) {
        Teuchos::RCP<Matrix> xpmat = mat->getMatrix(r,c);
        Teuchos::RCP<CrsMatrixWrap> xpwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(xpmat);
        TEUCHOS_TEST_FOR_EXCEPT(xpwrap.is_null() == true);

        Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thBlock =
            Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpwrap->getCrsMatrix());
        std::stringstream label; label << "A" << r << c;
        blockMat->setBlock(r,c,thBlock);
      }
    }

    blockMat->endBlockFill();

    return blockMat;
  }

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

  // Initialize MPI
#ifdef HAVE_MPI
  MPI_Init(&argc,&argv);
#endif

  // Create a communicator for Epetra objects
#ifdef HAVE_MPI
  Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
  Epetra_SerialComm Comm;
#endif

  // Get the process ID and the total number of processors
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  // Check verbosity level
  bool verbose = false;
  if (argc > 1)
    if (argv[1][0]=='-' && argv[1][1]=='v')
      verbose = true;

  // Get the number of elements from the command line
  int NumGlobalElements = 0;
  if ((argc > 2) && (verbose))
    NumGlobalElements = atoi(argv[2]) + 1;
  else if ((argc > 1) && (!verbose))
    NumGlobalElements = atoi(argv[1]) + 1;
  else
    NumGlobalElements = 101;

  // The number of unknowns must be at least equal to the
  // number of processors.
  if (NumGlobalElements < NumProc) {
    std::cout << "numGlobalBlocks = " << NumGlobalElements
     << " cannot be < number of processors = " << NumProc << std::endl;
    std::cout << "Test failed!" << std::endl;
    throw "NOX Error";
  }

  // Create the interface between NOX and the application
  // This object is derived from NOX::Epetra::Interface
  Teuchos::RCP<Interface> interface =
    Teuchos::rcp(new Interface(NumGlobalElements, Comm));

  // Set the PDE factor (for nonlinear forcing term).  This could be specified
  // via user input.
  interface->setPDEfactor(1000.0);

  // Use a scaled vector space.  The scaling must also be registered
  // with the linear solver so the linear system is consistent!
  Teuchos::RCP<Epetra_Vector> scaleVec =
    Teuchos::rcp(new Epetra_Vector( *(interface->getSolution())));
  scaleVec->PutScalar(2.0);
  Teuchos::RCP<NOX::Epetra::Scaling> scaling =
    Teuchos::rcp(new NOX::Epetra::Scaling);
  scaling->addUserScaling(NOX::Epetra::Scaling::Left, scaleVec);

  // Use a weighted vector space for scaling all norms
  Teuchos::RCP<NOX::Epetra::VectorSpace> weightedVectorSpace =
    Teuchos::rcp(new NOX::Epetra::VectorSpaceScaledL2(scaling));

  // Get the vector from the Problem
  Teuchos::RCP<Epetra_Vector> soln = interface->getSolution();
  Teuchos::RCP<NOX::Epetra::Vector> noxSoln =
    Teuchos::rcp(new NOX::Epetra::Vector(soln,
                     NOX::Epetra::Vector::CreateCopy,
                     NOX::DeepCopy,
                     weightedVectorSpace));

  // Begin Nonlinear Solver ************************************

  // Create the top level parameter list
  Teuchos::RCP<Teuchos::ParameterList> nlParamsPtr =
    Teuchos::rcp(new Teuchos::ParameterList);
  Teuchos::ParameterList& nlParams = *(nlParamsPtr.get());

  // Set the nonlinear solver method
  nlParams.set("Nonlinear Solver", "Line Search Based");

  // Set the printing parameters in the "Printing" sublist
  Teuchos::ParameterList& printParams = nlParams.sublist("Printing");
  printParams.set("MyPID", MyPID);
  printParams.set("Output Precision", 3);
  printParams.set("Output Processor", 0);
  if (verbose)
    printParams.set("Output Information",
                 NOX::Utils::OuterIteration +
                 NOX::Utils::OuterIterationStatusTest +
                 NOX::Utils::InnerIteration +
                 NOX::Utils::LinearSolverDetails +
                 NOX::Utils::Parameters +
                 NOX::Utils::Details +
                 NOX::Utils::Warning +
                             NOX::Utils::Debug +
                 NOX::Utils::TestDetails +
                 NOX::Utils::Error);
//.........这里部分代码省略.........

// hide the original parental method AMET->evalModelImpl():
void
Aeras::HVDecorator::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<ST>& inArgsT,
    const Thyra::ModelEvaluatorBase::OutArgs<ST>& outArgsT) const
{

  std::cout << "DEBUG WHICH HVDecorator: " << __PRETTY_FUNCTION__ << "\n";
	
  Teuchos::TimeMonitor Timer(*timer); //start timer
  //
  // Get the input arguments
  //
  const Teuchos::RCP<const Tpetra_Vector> xT =
    ConverterT::getConstTpetraVector(inArgsT.get_x());

  const Teuchos::RCP<const Tpetra_Vector> x_dotT =
    Teuchos::nonnull(inArgsT.get_x_dot()) ?
    ConverterT::getConstTpetraVector(inArgsT.get_x_dot()) :
    Teuchos::null;

  // AGS: x_dotdot time integrators not imlemented in Thyra ME yet
  //const Teuchos::RCP<const Tpetra_Vector> x_dotdotT =
  //  Teuchos::nonnull(inArgsT.get_x_dotdot()) ?
  //  ConverterT::getConstTpetraVector(inArgsT.get_x_dotdot()) :
  //  Teuchos::null;
  const Teuchos::RCP<const Tpetra_Vector> x_dotdotT = Teuchos::null;


  const double alpha = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_alpha() : 0.0;
  // AGS: x_dotdot time integrators not imlemented in Thyra ME yet
  // const double omega = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_omega() : 0.0;
  const double omega = 0.0;
  const double beta = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_beta() : 1.0;
  const double curr_time = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_t() : 0.0;

  for (int l = 0; l < inArgsT.Np(); ++l) {
    const Teuchos::RCP<const Thyra::VectorBase<ST> > p = inArgsT.get_p(l);
    if (Teuchos::nonnull(p)) {
      const Teuchos::RCP<const Tpetra_Vector> pT = ConverterT::getConstTpetraVector(p);
      const Teuchos::ArrayRCP<const ST> pT_constView = pT->get1dView();

      ParamVec &sacado_param_vector = sacado_param_vec[l];
      for (unsigned int k = 0; k < sacado_param_vector.size(); ++k) {
        sacado_param_vector[k].baseValue = pT_constView[k];
      }
    }
  }

  //
  // Get the output arguments
  //
  const Teuchos::RCP<Tpetra_Vector> fT_out =
    Teuchos::nonnull(outArgsT.get_f()) ?
    ConverterT::getTpetraVector(outArgsT.get_f()) :
    Teuchos::null;

  const Teuchos::RCP<Tpetra_Operator> W_op_outT =
    Teuchos::nonnull(outArgsT.get_W_op()) ?
    ConverterT::getTpetraOperator(outArgsT.get_W_op()) :
    Teuchos::null;

#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
  //IK, 4/24/15: adding object to hold mass matrix to be written to matrix market file
  const Teuchos::RCP<Tpetra_Operator> Mass =
    Teuchos::nonnull(outArgsT.get_W_op()) ?
    ConverterT::getTpetraOperator(outArgsT.get_W_op()) :
    Teuchos::null;
  //IK, 4/24/15: needed for writing mass matrix out to matrix market file
  const Teuchos::RCP<Tpetra_Vector> ftmp =
    Teuchos::nonnull(outArgsT.get_f()) ?
    ConverterT::getTpetraVector(outArgsT.get_f()) :
    Teuchos::null;
#endif

  // Cast W to a CrsMatrix, throw an exception if this fails
  const Teuchos::RCP<Tpetra_CrsMatrix> W_op_out_crsT =
    Teuchos::nonnull(W_op_outT) ?
    Teuchos::rcp_dynamic_cast<Tpetra_CrsMatrix>(W_op_outT, true) :
    Teuchos::null;

#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
  //IK, 4/24/15: adding object to hold mass matrix to be written to matrix market file
  const Teuchos::RCP<Tpetra_CrsMatrix> Mass_crs =
    Teuchos::nonnull(Mass) ?
    Teuchos::rcp_dynamic_cast<Tpetra_CrsMatrix>(Mass, true) :
    Teuchos::null;
#endif

  //
  // Compute the functions
  //
  bool f_already_computed = false;

  // W matrix
  if (Teuchos::nonnull(W_op_out_crsT)) {
    app->computeGlobalJacobianT(
        alpha, beta, omega, curr_time, x_dotT.get(), x_dotdotT.get(),  *xT,
        sacado_param_vec, fT_out.get(), *W_op_out_crsT);
    f_already_computed = true;
//.........这里部分代码省略.........

 RiskVector( const Teuchos::RCP<Vector<Real> > &vec,
             const bool augmented = false,
             const Real stat = 0.)
   : stat_(stat), vec_(vec), augmented_(augmented) {
   dual_vec1_ = vec->dual().clone();
 }