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

Aeras::ShallowWaterResponseL2Error<EvalT, Traits>::
ShallowWaterResponseL2Error(Teuchos::ParameterList& p,
		      const Teuchos::RCP<Albany::Layouts>& dl) :
  sphere_coord("Lat-Long", dl->qp_gradient),
  weighted_measure("Weights", dl->qp_scalar),
  flow_state_field("Flow State", dl->node_vector), 
  BF("BF",dl->node_qp_scalar)
{
  Teuchos::RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());
  // get and validate Response parameter list
  Teuchos::ParameterList* plist = 
    p.get<Teuchos::ParameterList*>("Parameter List");
  Teuchos::RCP<const Teuchos::ParameterList> reflist = 
    this->getValidResponseParameters();
  plist->validateParameters(*reflist,0);

  std::string fieldName = "Flow Field"; //field to integral is the flow field 

  // coordinate dimensions
  std::vector<PHX::DataLayout::size_type> coord_dims;
  dl->qp_vector->dimensions(coord_dims);
  numQPs = coord_dims[1]; //# quad points
  numDims = coord_dims[2]; //# spatial dimensions
  std::vector<PHX::DataLayout::size_type> dims;
  flow_state_field.fieldTag().dataLayout().dimensions(dims);
  vecDim = dims[2]; //# dofs per node
  numNodes =  dims[1]; //# nodes per element

 
  // User-specified parameters
  refSolName = plist->get<std::string>("Reference Solution Name"); //no reference solution by default.
  *out << "Reference Solution Name for Aeras::ShallowWaterResponseL2Error response: " << refSolName << std::endl; 
  inputData = plist->get<double>("Reference Solution Data", 0.0);
   
  if (refSolName == "Zero")
    ref_sol_name = ZERO;
  else if (refSolName == "TC2")
    ref_sol_name  = TC2;
  //Add other test case reference solutions here...
  else if (refSolName == "TC4"){
    ref_sol_name = TC4;
    
    myPi = Aeras::ShallowWaterConstants::self().pi;
    earthRadius = Aeras::ShallowWaterConstants::self().earthRadius;
    gravity = Aeras::ShallowWaterConstants::self().gravity;
    
    Omega = 2.0*myPi/(24.*3600.); //this should be sitting in SW Constants class
    
    rlon0 = 0.;
    rlat0 = myPi/4.;
    npwr = 14.;
    
    su0 = 20.;
    phi0 = 1.0e5;
    alfa = -0.03*(phi0/(2.*Omega*sin(myPi/4.)));
    sigma = (2.*earthRadius/1.0e6)*(2.*earthRadius/1.0e6);
    
  }
  else { 
    TEUCHOS_TEST_FOR_EXCEPTION(
      true, Teuchos::Exceptions::InvalidParameter,
      std::endl << "Error!  Unknown reference solution name " << ref_sol_name <<
      "!" << std::endl;);
  }

Teuchos::ParameterList Ifpack_GetValidParameters()
{
  Teuchos::ParameterList List; // empty list

  // ============================================================ //
  // Parameters are reported from each used file in IFPACK. Files //
  // are listed in alphabetical order, first all *.cpp, then *.h. //
  // Some options not very tested or documented anywhere          //
  // are not reported here.                                       //
  // ============================================================ //
  
  // Ifpack_Amesos.cpp
  List.set("amesos: solver type", "Amesos_Klu");

  // Ifpack_IC.cpp
  List.set("fact: level-of-fill", (int)1);
  List.set("fact: absolute threshold", (double)0.0);
  List.set("fact: relative threshold", (double)0.0);
  List.set("fact: drop tolerance", (double)0.0);

  // Ifpack_ICT.cpp
  List.set("fact: ict level-of-fill", (double)1.0);
  List.set("fact: absolute threshold", (double)0.0);
  List.set("fact: relative threshold", (double)1.0);
  List.set("fact: relax value", (double)0.0);
  List.set("fact: drop tolerance", (double)0.0);

  // Ifpack_ILU.cpp
  List.set("fact: level-of-fill", (int)0);
  List.set("fact: absolute threshold", (double)0.0);
  List.set("fact: relative threshold", (double)1.0);
  List.set("fact: relax value", (double)0.0);

  // Ifpack_ILUT.cpp
  List.set("fact: ilut level-of-fill", (double)1.0);
  List.set("fact: absolute threshold", (double)0.0);
  List.set("fact: relative threshold", (double)1.0);
  List.set("fact: relax value", (double)0.0);

#ifdef HAVE_IFPACK_SUPERLU
  // Ifpack_SILU.cpp
  List.set("fact: drop tolerance",1e-4);
  List.set("fact: zero pivot threshold",1e-2);
  List.set("fact: maximum fill factor",10.0);
  List.set("fact: silu drop rule",9);
#endif

  // Ifpack_METISPartitioner.cpp
  List.set("partitioner: local parts", (int)1);
  List.set("partitioner: overlap", (int)0);
  List.set("partitioner: print level", (int)0);

  // Ifpack_PointRelaxation.cpp
  List.set("relaxation: type", "Jacobi");
  List.set("relaxation: sweeps", (int)1);
  List.set("relaxation: damping factor", (double)1.0);
  List.set("relaxation: min diagonal value", (double)1.0);
  List.set("relaxation: zero starting solution", true);
  List.set("relaxation: backward mode",false);
  List.set("relaxation: use l1",false);
  List.set("relaxation: l1 eta",(double)1.5);

  // Ifpack_SPARSKIT.cpp
  List.set("fact: sparskit: lfil", (int)0);
  List.set("fact: sparskit: tol", (double)0.0);
  List.set("fact: sparskit: droptol", (double)0.0);
  List.set("fact: sparskit: permtol", (double)0.1);
  List.set("fact: sparskit: alph", (double)0.0);
  List.set("fact: sparskit: mbloc", (int)(-1));
  List.set("fact: sparskit: type", ("ILUT"));

  // Additive Schwarz preconditioner
  List.set("schwarz: compute condest", true);
  List.set("schwarz: combine mode", "Zero"); // use std::string mode for this
  List.set("schwarz: reordering type", "none");
  List.set("schwarz: filter singletons", false);

  // Ifpack_BlockRelaxation.h
  // List.set("relaxation: type", "Jacobi"); // already set
  // List.set("relaxation: sweeps", 1); // already set
  // List.get("relaxation: damping factor", 1.0); // already set
  // List.get("relaxation: zero starting solution", true); // already set
  List.set("partitioner: type", "greedy");
  List.set("partitioner: local parts", (int)1);
  List.set("partitioner: overlap", (int)0);

  // Ifpack_METISPartitioner.h
  List.set("partitioner: use symmetric graph", true);

  // Krylov smoother
  List.set("krylov: iterations",(int)5);
  List.set("krylov: tolerance",(double)0.001);
  List.set("krylov: solver",(int)1);
  List.set("krylov: preconditioner",(int)0);
  List.set("krylov: number of sweeps",(int)1);
  List.set("krylov: block size",(int)1);
  List.set("krylov: damping parameter",(double)1.0);
  List.set("krylov: zero starting solution",true);

  return(List);
//.........这里部分代码省略.........

//.........这里部分代码省略.........

    // ------------------------------------------------------------- //
    // complete the interfaces
    // ------------------------------------------------------------- //
    for (int i=0; i<ninter; ++i)
      if (!interfaces[i]->Complete())
      {
         cout << "Interface " << i << " completion returned false\n";
         exit(EXIT_FAILURE);
      }

    // ------------------------------------------------------------- //
    // create an empty MOERTEL::Manager for 3D problems
    // It organizes everything from integration to solution
    // ------------------------------------------------------------- //
    MOERTEL::Manager manager(Comm,MOERTEL::Manager::manager_3D,printlevel);
    
    // ------------------------------------------------------------- //
    // Add the interfaces to the manager
    // ------------------------------------------------------------- //
    for (int i=0; i<ninter; ++i)
      manager.AddInterface(*(interfaces[i]));

    // ------------------------------------------------------------- //
    // for mortar integration, the mortar manager needs to know about
    // the rowmap of the original (uncoupled) problem because it will
    // create coupling matrices D and M matching that rowmap
    // ------------------------------------------------------------- //
    manager.SetProblemMap(&Grid.RowMap());

    // ============================================================= //
    // choose integration parameters
    // ============================================================= //
    Teuchos::ParameterList& moertelparams = manager.Default_Parameters();
    // this does affect this 3D case only
    moertelparams.set("exact values at gauss points",true);
    // 1D interface possible values are 1,2,3,4,5,6,7,8,10 (2 recommended with linear shape functions)
    moertelparams.set("number gaussian points 1D",2);
    // 2D interface possible values are 3,6,12,13,16,19,27 (12 recommended with linear functions)
    moertelparams.set("number gaussian points 2D",12);

    // ============================================================= //
    // Here we are done with the construction phase of the interface
    // so we can integrate the mortar integrals
    // (Note we have not yet evaluated the PDE at all!)
    // ============================================================= //
    manager.Mortar_Integrate();
    
    // print interface information
    // (Manager, Interface, Segment, Node implement the << operator)
    if (printlevel) cout << manager;
        
    // ======================================================== //
    // Prepares the linear system. This requires the definition //
    // of a quadrature formula compatible with the grid, a      //
    // variational formulation, and a problem object which take //
    // care of filling matrix and right-hand side.              //
    // NOTE:
    // we are doing this AFTER we did all the mortar stuff to
    // show that the mortar integration is actually PDE-independent
    // ======================================================== //
    Epetra_CrsMatrix A(Copy, Grid.RowMap(), 0);
    Epetra_Vector    LHS(Grid.RowMap(),true);
    Epetra_Vector    RHS(Grid.RowMap());

    int NumQuadratureNodes = 8;

  void ParameterListInterpreter<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps>::UpdateFactoryManager(Teuchos::ParameterList& paramList,
        const Teuchos::ParameterList& defaultList, FactoryManager& manager) const {
    // NOTE: Factory::SetParameterList must be called prior to Factory::SetFactory, as
    // SetParameterList sets default values for non mentioned parameters, including factories

    // === Smoothing ===
    bool isCustomSmoother =
        paramList.isParameter("smoother: pre or post") ||
        paramList.isParameter("smoother: type")    || paramList.isParameter("smoother: pre type")    || paramList.isParameter("smoother: post type")   ||
        paramList.isSublist  ("smoother: params")  || paramList.isSublist  ("smoother: pre params")  || paramList.isSublist  ("smoother: post params") ||
        paramList.isParameter("smoother: sweeps")  || paramList.isParameter("smoother: pre sweeps")  || paramList.isParameter("smoother: post sweeps") ||
        paramList.isParameter("smoother: overlap") || paramList.isParameter("smoother: pre overlap") || paramList.isParameter("smoother: post overlap");;
    MUELU_READ_2LIST_PARAM(paramList, defaultList, "smoother: pre or post", std::string, "both", PreOrPost);
    if (PreOrPost == "none") {
      manager.SetFactory("Smoother", Teuchos::null);

    } else if (isCustomSmoother) {
      // FIXME: get default values from the factory
      // NOTE: none of the smoothers at the moment use parameter validation framework, so we
      // cannot get the default values from it.
#define TEST_MUTUALLY_EXCLUSIVE(arg1,arg2) \
      TEUCHOS_TEST_FOR_EXCEPTION(paramList.isParameter(#arg1) && paramList.isParameter(#arg2), \
                                 Exceptions::InvalidArgument, "You cannot specify both \""#arg1"\" and \""#arg2"\"");
#define TEST_MUTUALLY_EXCLUSIVE_S(arg1,arg2) \
      TEUCHOS_TEST_FOR_EXCEPTION(paramList.isSublist(#arg1) && paramList.isSublist(#arg2), \
                                 Exceptions::InvalidArgument, "You cannot specify both \""#arg1"\" and \""#arg2"\"");

      TEST_MUTUALLY_EXCLUSIVE  ("smoother: type",    "smoother: pre type");
      TEST_MUTUALLY_EXCLUSIVE  ("smoother: type",    "smoother: post type");
      TEST_MUTUALLY_EXCLUSIVE  ("smoother: sweeps",  "smoother: pre sweeps");
      TEST_MUTUALLY_EXCLUSIVE  ("smoother: sweeps",  "smoother: post sweeps");
      TEST_MUTUALLY_EXCLUSIVE  ("smoother: overlap", "smoother: pre overlap");
      TEST_MUTUALLY_EXCLUSIVE  ("smoother: overlap", "smoother: post overlap");
      TEST_MUTUALLY_EXCLUSIVE_S("smoother: params",  "smoother: pre params");
      TEST_MUTUALLY_EXCLUSIVE_S("smoother: params",  "smoother: post params");
      TEUCHOS_TEST_FOR_EXCEPTION(PreOrPost == "both" && (paramList.isParameter("smoother: pre type") != paramList.isParameter("smoother: post type")),
                                 Exceptions::InvalidArgument, "You must specify both \"smoother: pre type\" and \"smoother: post type\"");

      // Default values
      int overlap = 0;
      ParameterList defaultSmootherParams;
      defaultSmootherParams.set("relaxation: type",           "Symmetric Gauss-Seidel");
      defaultSmootherParams.set("relaxation: sweeps",         Teuchos::OrdinalTraits<LO>::one());
      defaultSmootherParams.set("relaxation: damping factor", Teuchos::ScalarTraits<Scalar>::one());

      RCP<SmootherPrototype> preSmoother = Teuchos::null, postSmoother = Teuchos::null;
      std::string            preSmootherType,             postSmootherType;
      ParameterList          preSmootherParams,           postSmootherParams;

      if (paramList.isParameter("smoother: overlap"))
        overlap = paramList.get<int>("smoother: overlap");

      if (PreOrPost == "pre" || PreOrPost == "both") {
        if (paramList.isParameter("smoother: pre type")) {
          preSmootherType = paramList.get<std::string>("smoother: pre type");
        } else {
          MUELU_READ_2LIST_PARAM(paramList, defaultList, "smoother: type", std::string, "RELAXATION", preSmootherTypeTmp);
          preSmootherType = preSmootherTypeTmp;
        }
        if (paramList.isParameter("smoother: pre overlap"))
          overlap = paramList.get<int>("smoother: pre overlap");

        if (paramList.isSublist("smoother: pre params"))
          preSmootherParams = paramList.sublist("smoother: pre params");
        else if (paramList.isSublist("smoother: params"))
          preSmootherParams = paramList.sublist("smoother: params");
        else if (defaultList.isSublist("smoother: params"))
          preSmootherParams = defaultList.sublist("smoother: params");
        else if (preSmootherType == "RELAXATION")
          preSmootherParams = defaultSmootherParams;

        preSmoother = rcp(new TrilinosSmoother(preSmootherType, preSmootherParams, overlap));
      }

      if (PreOrPost == "post" || PreOrPost == "both") {
        if (paramList.isParameter("smoother: post type"))
          postSmootherType = paramList.get<std::string>("smoother: post type");
        else {
          MUELU_READ_2LIST_PARAM(paramList, defaultList, "smoother: type", std::string, "RELAXATION", postSmootherTypeTmp);
          postSmootherType = postSmootherTypeTmp;
        }

        if (paramList.isSublist("smoother: post params"))
          postSmootherParams = paramList.sublist("smoother: post params");
        else if (paramList.isSublist("smoother: params"))
          postSmootherParams = paramList.sublist("smoother: params");
        else if (defaultList.isSublist("smoother: params"))
          postSmootherParams = defaultList.sublist("smoother: params");
        else if (postSmootherType == "RELAXATION")
          postSmootherParams = defaultSmootherParams;
        if (paramList.isParameter("smoother: post overlap"))
          overlap = paramList.get<int>("smoother: post overlap");

        if (postSmootherType == preSmootherType && areSame(preSmootherParams, postSmootherParams))
          postSmoother = preSmoother;
        else
          postSmoother = rcp(new TrilinosSmoother(postSmootherType, postSmootherParams, overlap));
      }

      manager.SetFactory("Smoother", rcp(new SmootherFactory(preSmoother, postSmoother)));
//.........这里部分代码省略.........

 StochasticStefanBoltzmannMeshManager(Teuchos::ParameterList &parlist)
   : MeshManager_Rectangle<Real>(parlist) {
   nx_ = parlist.sublist("Geometry").get("NX", 3);
   ny_ = parlist.sublist("Geometry").get("NY", 3);
   computeSideSets();
 }

void getValidParameters(Teuchos::ParameterList& params)
{
  //params.clear();
  Teuchos::ParameterList empty;
  params = empty;

  // ============================================================ //
  // Parameters are reported from each used file in IFPACK2. Files //
  // are listed in alphabetical order, first all *.cpp, then *.hpp. //
  // Some options not very tested or documented anywhere          //
  // are not reported here.                                       //
  // ============================================================ //
  
  // Ifpack2_IlukGraph.hpp
  params.set("fact: iluk level-of-fill", (int)1);
  params.set("fact: iluk level-of-overlap", (int)0);

  // Ifpack2_Amesos.cpp
  params.set("amesos: solver type", "Amesos_Klu");

  // Ifpack2_IC.cpp
  params.set("fact: level-of-fill", (int)1);
  params.set("fact: absolute threshold", (double)0.0);
  params.set("fact: relative threshold", (double)0.0);
  params.set("fact: drop tolerance", (double)0.0);

  // Ifpack2_ICT.cpp
  params.set("fact: ict level-of-fill", (double)1.0);
  params.set("fact: absolute threshold", (double)0.0);
  params.set("fact: relative threshold", (double)1.0);
  params.set("fact: relax value", (double)0.0);
  params.set("fact: drop tolerance", (double)0.0);

  // Ifpack2_ILU.cpp
  params.set("fact: level-of-fill", (int)0);
  params.set("fact: absolute threshold", (double)0.0);
  params.set("fact: relative threshold", (double)1.0);
  params.set("fact: relax value", (double)0.0);

  // Ifpack2_ILUT.cpp
  params.set("fact: ilut level-of-fill", (double)1.0);
  params.set("fact: absolute threshold", (double)0.0);
  params.set("fact: relative threshold", (double)1.0);
  params.set("fact: relax value", (double)0.0);

  // Ifpack2_METISPartitioner.cpp
  params.set("partitioner: local parts", (int)1);
  params.set("partitioner: overlap", (int)0);
  params.set("partitioner: print level", (int)0);

  // Ifpack2_Relaxation.cpp
  params.set("relaxation: type", "Jacobi");
  params.set("relaxation: sweeps", (int)1);
  params.set("relaxation: damping factor", (double)1.0);
  params.set("relaxation: min diagonal value", (double)1.0);
  params.set("relaxation: zero starting solution", true);

  // Ifpack2_SPARSKIT.cpp
  params.set("fact: sparskit: lfil", (int)0);
  params.set("fact: sparskit: tol", (double)0.0);
  params.set("fact: sparskit: droptol", (double)0.0);
  params.set("fact: sparskit: permtol", (double)0.1);
  params.set("fact: sparskit: alph", (double)0.0);
  params.set("fact: sparskit: mbloc", (int)(-1));
  params.set("fact: sparskit: type", ("ILUT"));

  // Additive Schwarz preconditioner
  params.set("schwarz: compute condest", true);
  params.set("schwarz: combine mode", "Zero"); // use string mode for this
  params.set("schwarz: reordering type", "none");
  params.set("schwarz: filter singletons", false);

  // Ifpack2_BlockRelaxation.hpp
  // params.set("relaxation: type", "Jacobi"); // already set
  // params.set("relaxation: sweeps", 1); // already set
  // params.get("relaxation: damping factor", 1.0); // already set
  // params.get("relaxation: zero starting solution", true); // already set
  params.set("partitioner: type", "greedy");
  params.set("partitioner: local parts", (int)1);
  params.set("partitioner: overlap", (int)0);

  // Ifpack2_METISPartitioner.hpp
  params.set("partitioner: use symmetric graph", true);
}

 // Constructor
 Brents( Teuchos::ParameterList &parlist ) : LineSearch<Real>(parlist) {
   tol_ = parlist.get("Bracketing Tolerance",1.e-8);
   btls_ = Teuchos::rcp(new BackTracking<Real>(parlist));
 }

//.........这里部分代码省略.........
    Teuchos::RCP<ROL::SampleGenerator<RealT> > sampler
      = Teuchos::rcp(new ROL::MonteCarloGenerator<RealT>(
          nSamp,bounds,bman,false,false,100));
    /*************************************************************************/
    /************* INITIALIZE REDUCED OBJECTIVE FUNCTION *********************/
    /*************************************************************************/
    bool storage = true, fdhess = false;
    Teuchos::RCP<ROL::Objective<RealT> > robj
      = Teuchos::rcp(new ROL::Reduced_Objective_SimOpt<RealT>(
          pobj,pcon,up,zp,lp,gup,gzp,cp,storage,fdhess));
    /*************************************************************************/
    /************* INITIALIZE BOUND CONSTRAINTS ******************************/
    /*************************************************************************/
    std::vector<RealT> Zlo(nx+2,0.0), Zhi(nx+2,10.0);
    for (int i = 0; i < nx+2; i++) {
      if ( i < (int)((nx+2)/3) ) {
        Zlo[i] = -1.0;
        Zhi[i] = 1.0;
      }
      if ( i >= (int)((nx+2)/3) && i < (int)(2*(nx+2)/3) ) {
        Zlo[i] = 1.0;
        Zhi[i] = 5.0;
      }
      if ( i >= (int)(2*(nx+2)/3) ) {
        Zlo[i] = 5.0;
        Zhi[i] = 10.0;
      }
    }
    Teuchos::RCP<ROL::BoundConstraint<RealT> > Zbnd
      = Teuchos::rcp(new L2BoundConstraint<RealT>(Zlo,Zhi,fem));
    /*************************************************************************/
    /************* INITIALIZE OPTIMIZATION PROBLEM ***************************/
    /*************************************************************************/
    Teuchos::ParameterList SOLlist;
    SOLlist.sublist("SOL").set("Stochastic Optimization Type","Risk Averse");
    SOLlist.sublist("SOL").set("Store Sampled Value and Gradient",storage);
    SOLlist.sublist("SOL").sublist("Risk Measure").set("Name","KL Divergence");
    SOLlist.sublist("SOL").sublist("Risk Measure").sublist("KL Divergence").set("Threshold",1.e-2);
    ROL::StochasticProblem<RealT> optProb(SOLlist,robj,sampler,zp,Zbnd);
    /*************************************************************************/
    /************* CHECK DERIVATIVES AND CONSISTENCY *************************/
    /*************************************************************************/
    // CHECK OBJECTIVE DERIVATIVES
    bool derivcheck = false;
    if (derivcheck) {
      int nranks = sampler->numBatches();
      for (int pid = 0; pid < nranks; pid++) {
        if ( pid == sampler->batchID() ) {
          for (int i = sampler->start(); i < sampler->numMySamples(); i++) {
            *outStream << "Sample " << i << "  Rank " << sampler->batchID() << "\n";
            *outStream << "(" << sampler->getMyPoint(i)[0] << ", "
                              << sampler->getMyPoint(i)[1] << ", "
                              << sampler->getMyPoint(i)[2] << ", "
                              << sampler->getMyPoint(i)[3] << ")\n";
            pcon->setParameter(sampler->getMyPoint(i));
            pcon->checkSolve(*up,*zp,*cp,print,*outStream);
            robj->setParameter(sampler->getMyPoint(i));
            *outStream << "\n";
            robj->checkGradient(*zp,*gzp,*yzp,print,*outStream);
            robj->checkHessVec(*zp,*gzp,*yzp,print,*outStream);
            *outStream << "\n\n";
          }
        }
        comm->barrier();
      }
    }

// ***********************************************************
int DG_Prob::Eigenvectors(const double Dt,
                          const Epetra_Map & Map)
{
  printf("Entrou em Eigenvectors\n");
 
#ifdef HAVE_MPI
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  //MPI::COMM_WORLD.Barrier();
  Comm.Barrier();
  Teuchos::RCP<Epetra_FECrsMatrix> M = Teuchos::rcp(new Epetra_FECrsMatrix(Copy, Map,0));//&NNz[0]);
  Teuchos::RCP<Epetra_FEVector> RHS = Teuchos::rcp(new Epetra_FEVector(Map,1));
  
  DG_MatrizVetor_Epetra(Dt,M,RHS);

  Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp(new Epetra_CrsMatrix(Copy, Map,0
                                                                       /* &NNz[0]*/) );
  Epetra_Export Exporter(Map,Map);
  A->PutScalar(0.0);
  A->Export(*(M.ptr()),Exporter,Add);
  A->FillComplete();

  using std::cout;
 // int nx = 5;
  bool boolret;
  int MyPID = Comm.MyPID();
  
  bool verbose = true;
  bool debug = false;
  std::string which("LR");
  
  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
  cmdp.setOption("sort",&which,"Targetted eigenvalues (SM,LM,SR,LR,SI,or LI).");
  
  typedef double ScalarType;
  typedef Teuchos::ScalarTraits<ScalarType>          SCT;
  typedef SCT::magnitudeType               MagnitudeType;
  typedef Epetra_MultiVector                          MV;
  typedef Epetra_Operator                             OP;
  typedef Anasazi::MultiVecTraits<ScalarType,MV>     MVT;
  typedef Anasazi::OperatorTraits<ScalarType,MV,OP>  OPT;
  
 
 // double rho = 2*nx+1;
  
  // Compute coefficients for discrete convection-diffution operator
 // const double one = 1.0;
  
 // int NumEntries, info;
  
  //************************************
  // Start the block Arnoldi iteration
  //***********************************
  //
  //  Variables used for the Block Krylov Schur Method
  //    
  int nev = 10;
  int blockSize = 1;
  int numBlocks = 20;
  int maxRestarts = 500;
  //int stepSize = 5;
  double tol = 1e-8;
  
  // Create a sort manager to pass into the block Krylov-Schur solver manager
  // -->  Make sure the reference-counted pointer is of type Anasazi::SortManager<>
  // -->  The block Krylov-Schur solver manager uses Anasazi::BasicSort<> by default,
  //      so you can also pass in the parameter "Which", instead of a sort manager.
  Teuchos::RCP<Anasazi::SortManager<MagnitudeType> > MySort =     
    Teuchos::rcp( new Anasazi::BasicSort<MagnitudeType>( which ) );
  
  // Set verbosity level
  int verbosity = Anasazi::Errors + Anasazi::Warnings;
  if (verbose) {
    verbosity += Anasazi::FinalSummary + Anasazi::TimingDetails;
  }
  if (debug) {
    verbosity += Anasazi::Debug;
  }
  //
  // Create parameter list to pass into solver manager
  //
  Teuchos::ParameterList MyPL;
  MyPL.set( "Verbosity", verbosity );
  MyPL.set( "Sort Manager", MySort );
  //MyPL.set( "Which", which );  
  MyPL.set( "Block Size", blockSize );
  MyPL.set( "Num Blocks", numBlocks );
  MyPL.set( "Maximum Restarts", maxRestarts );
  //MyPL.set( "Step Size", stepSize );
  MyPL.set( "Convergence Tolerance", tol );
  
  // Create an Epetra_MultiVector for an initial vector to start the solver.
  // Note:  This needs to have the same number of columns as the blocksize.
  Teuchos::RCP<Epetra_MultiVector> ivec = Teuchos::rcp( new Epetra_MultiVector(Map, blockSize) );
//.........这里部分代码省略.........

TEUCHOS_UNIT_TEST(gs_evaluators, gather_constr)
{
    PHX::KokkosDeviceSession session;

    const std::size_t workset_size = 20;
    Teuchos::RCP<panzer::BasisIRLayout> linBasis = buildLinearBasis(workset_size);

    Teuchos::RCP<std::vector<std::string> > fieldNames
        = Teuchos::rcp(new std::vector<std::string>);
    fieldNames->push_back("dog");

    Teuchos::ParameterList pl;
    pl.set("Basis",linBasis);
    pl.set("Field Names",fieldNames);

    Teuchos::RCP<panzer_stk_classic::STK_Interface> mesh = buildMesh(2,2);

    RCP<Epetra_Comm> Comm = Teuchos::rcp(new Epetra_MpiComm(MPI_COMM_WORLD));

    Teuchos::RCP<Teuchos::ParameterList> ipb = Teuchos::parameterList("Physics Blocks");
    std::vector<panzer::BC> bcs;
    testInitialzation(ipb, bcs);

    Teuchos::RCP<panzer::FieldManagerBuilder> fmb =
        Teuchos::rcp(new panzer::FieldManagerBuilder);

    // build physics blocks
    //////////////////////////////////////////////////////////////
    Teuchos::RCP<user_app::MyFactory> eqset_factory = Teuchos::rcp(new user_app::MyFactory);
    user_app::BCFactory bc_factory;
    std::vector<Teuchos::RCP<panzer::PhysicsBlock> > physicsBlocks;

    {
        std::map<std::string,std::string> block_ids_to_physics_ids;
        block_ids_to_physics_ids["eblock-0_0"] = "test physics";
        block_ids_to_physics_ids["eblock-1_0"] = "test physics";

        std::map<std::string,Teuchos::RCP<const shards::CellTopology> > block_ids_to_cell_topo;
        block_ids_to_cell_topo["eblock-0_0"] = mesh->getCellTopology("eblock-0_0");
        block_ids_to_cell_topo["eblock-1_0"] = mesh->getCellTopology("eblock-1_0");

        Teuchos::RCP<panzer::GlobalData> gd = panzer::createGlobalData();

        int default_integration_order = 1;

        panzer::buildPhysicsBlocks(block_ids_to_physics_ids,
                                   block_ids_to_cell_topo,
                                   ipb,
                                   default_integration_order,
                                   workset_size,
                                   eqset_factory,
                                   gd,
                                   false,
                                   physicsBlocks);
    }

    Teuchos::RCP<panzer_stk_classic::WorksetFactory> wkstFactory
        = Teuchos::rcp(new panzer_stk_classic::WorksetFactory(mesh)); // build STK workset factory
    Teuchos::RCP<panzer::WorksetContainer> wkstContainer     // attach it to a workset container (uses lazy evaluation)
        = Teuchos::rcp(new panzer::WorksetContainer(wkstFactory,physicsBlocks,workset_size));

    // build DOF Manager
    /////////////////////////////////////////////////////////////

    // build the connection manager
    const Teuchos::RCP<panzer::ConnManager<int,int> >
    conn_manager = Teuchos::rcp(new panzer_stk_classic::STKConnManager<int>(mesh));

    panzer::DOFManagerFactory<int,int> globalIndexerFactory;
    RCP<panzer::UniqueGlobalIndexer<int,int> > dofManager
        = globalIndexerFactory.buildUniqueGlobalIndexer(Teuchos::opaqueWrapper(MPI_COMM_WORLD),physicsBlocks,conn_manager);

    Teuchos::RCP<panzer::EpetraLinearObjFactory<panzer::Traits,int> > eLinObjFactory
        = Teuchos::rcp(new panzer::EpetraLinearObjFactory<panzer::Traits,int>(Comm.getConst(),dofManager));
    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > linObjFactory = eLinObjFactory;

    // setup field manager build
    /////////////////////////////////////////////////////////////

    // Add in the application specific closure model factory
    user_app::MyModelFactory_TemplateBuilder cm_builder;
    panzer::ClosureModelFactory_TemplateManager<panzer::Traits> cm_factory;
    cm_factory.buildObjects(cm_builder);

    Teuchos::ParameterList closure_models("Closure Models");
    closure_models.sublist("solid").sublist("SOURCE_TEMPERATURE").set<double>("Value",1.0);
    closure_models.sublist("ion solid").sublist("SOURCE_ION_TEMPERATURE").set<double>("Value",1.0);

    Teuchos::ParameterList user_data("User Data");

    fmb->setWorksetContainer(wkstContainer);
    fmb->setupVolumeFieldManagers(physicsBlocks,cm_factory,closure_models,*linObjFactory,user_data);
    fmb->setupBCFieldManagers(bcs,physicsBlocks,*eqset_factory,cm_factory,bc_factory,closure_models,*linObjFactory,user_data);

    fmb->writeVolumeGraphvizDependencyFiles("field_manager",physicsBlocks);

    panzer::AssemblyEngine_TemplateManager<panzer::Traits> ae_tm;
    panzer::AssemblyEngine_TemplateBuilder builder(fmb,linObjFactory);
    ae_tm.buildObjects(builder);

//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
#include <MueLu_UseShortNames.hpp>

  using Teuchos::RCP; // reference count pointers
  using Teuchos::rcp;
  using Teuchos::TimeMonitor;

  // =========================================================================
  // MPI initialization using Teuchos
  // =========================================================================
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
  RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
  int MyPID   = comm->getRank();
  int NumProc = comm->getSize();

  const Teuchos::RCP<Epetra_Comm> epComm = Teuchos::rcp_const_cast<Epetra_Comm>(Xpetra::toEpetra(comm));

  // =========================================================================
  // Convenient definitions
  // =========================================================================
  //SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();

  // Instead of checking each time for rank, create a rank 0 stream
  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
  Teuchos::FancyOStream& fancyout = *fancy;
  fancyout.setOutputToRootOnly(0);



  // =========================================================================
  // Parameters initialization
  // =========================================================================
  Teuchos::CommandLineProcessor clp(false);
  GO nx = 100, ny = 100;
  clp.setOption("nx",                   &nx,              "mesh size in x direction");
  clp.setOption("ny",                   &ny,              "mesh size in y direction");
  std::string xmlFileName = "xml/s3a.xml"; clp.setOption("xml", &xmlFileName,     "read parameters from a file. Otherwise, this example uses by default 'tutorial1a.xml'");
  int mgridSweeps = 1; clp.setOption("mgridSweeps", &mgridSweeps, "number of multigrid sweeps within Multigrid solver.");
  std::string printTimings = "no";   clp.setOption("timings", &printTimings,     "print timings to screen [yes/no]");
  double tol               = 1e-12;  clp.setOption("tol",                   &tol,              "solver convergence tolerance");

  switch (clp.parse(argc,argv)) {
    case Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED:        return EXIT_SUCCESS; break;
    case Teuchos::CommandLineProcessor::PARSE_ERROR:
    case Teuchos::CommandLineProcessor::PARSE_UNRECOGNIZED_OPTION: return EXIT_FAILURE; break;
    case Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL:                               break;
  }

  // =========================================================================
  // Problem construction
  // =========================================================================
  RCP<TimeMonitor> globalTimeMonitor = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time"))), tm;

  comm->barrier();
  tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));

  Teuchos::ParameterList GaleriList;
  GaleriList.set("nx", nx);
  GaleriList.set("ny", ny);
  GaleriList.set("mx", epComm->NumProc());
  GaleriList.set("my", 1);
  GaleriList.set("lx", 1.0); // length of x-axis
  GaleriList.set("ly", 1.0); // length of y-axis
  GaleriList.set("diff", 1e-5);
  GaleriList.set("conv", 1.0);

  // create map
  Teuchos::RCP<Epetra_Map> epMap = Teuchos::rcp(Galeri::CreateMap("Cartesian2D", *epComm, GaleriList));

  // create coordinates
  Teuchos::RCP<Epetra_MultiVector> epCoord = Teuchos::rcp(Galeri::CreateCartesianCoordinates("2D", epMap.get(), GaleriList));

  // create matrix
  Teuchos::RCP<Epetra_CrsMatrix> epA = Teuchos::rcp(Galeri::CreateCrsMatrix("Recirc2D", epMap.get(), GaleriList));

  // Epetra -> Xpetra
  Teuchos::RCP<CrsMatrix> exA = Teuchos::rcp(new Xpetra::EpetraCrsMatrix(epA));
  Teuchos::RCP<CrsMatrixWrap> exAWrap = Teuchos::rcp(new CrsMatrixWrap(exA));

  RCP<Matrix> A = Teuchos::rcp_dynamic_cast<Matrix>(exAWrap);
  int numPDEs = 1;
  A->SetFixedBlockSize(numPDEs);

  // set rhs and solution vector
  RCP<Epetra_Vector> B = Teuchos::rcp(new Epetra_Vector(*epMap));
  RCP<Epetra_Vector> X = Teuchos::rcp(new Epetra_Vector(*epMap));
  B->PutScalar(1.0);
  X->PutScalar(0.0);

  // Epetra -> Xpetra
  RCP<Vector> xB = Teuchos::rcp(new Xpetra::EpetraVector(B));
  RCP<Vector> xX = Teuchos::rcp(new Xpetra::EpetraVector(X));

  xX->setSeed(100);
  xX->randomize();

  // build null space vector
  RCP<const Map> map = A->getRowMap();
  RCP<MultiVector> nullspace = MultiVectorFactory::Build(map, numPDEs);

//.........这里部分代码省略.........

  Teuchos::RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
  CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > op,
                             const Teuchos::ParameterList& inParamList,
                             Teuchos::RCP<Xpetra::MultiVector<double, LocalOrdinal, GlobalOrdinal, Node> > coords = Teuchos::null,
                             Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > nullspace = Teuchos::null) {
    typedef MueLu::HierarchyManager<Scalar,LocalOrdinal,GlobalOrdinal,Node> HierarchyManager;
    typedef MueLu::HierarchyUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> HierarchyUtils;
    typedef MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> Hierarchy;
    typedef MueLu::MLParameterListInterpreter<Scalar,LocalOrdinal,GlobalOrdinal,Node> MLParameterListInterpreter;
    typedef MueLu::ParameterListInterpreter<Scalar,LocalOrdinal,GlobalOrdinal,Node> ParameterListInterpreter;
    typedef Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVectorFactory;

    std::string timerName = "MueLu setup time";
    RCP<Teuchos::Time> tm = Teuchos::TimeMonitor::getNewTimer(timerName);
    tm->start();

    bool hasParamList = inParamList.numParams();

    RCP<HierarchyManager> mueLuFactory;

    // Rip off non-serializable data before validation
    Teuchos::ParameterList nonSerialList,paramList;
    MueLu::ExtractNonSerializableData(inParamList, paramList, nonSerialList);

    std::string syntaxStr = "parameterlist: syntax";
    if (hasParamList && paramList.isParameter(syntaxStr) && paramList.get<std::string>(syntaxStr) == "ml") {
      paramList.remove(syntaxStr);
      mueLuFactory = rcp(new MLParameterListInterpreter(paramList));
    } else {
      mueLuFactory = rcp(new ParameterListInterpreter(paramList,op->getDomainMap()->getComm()));
    }

    // Create Hierarchy
    RCP<Hierarchy> H = mueLuFactory->CreateHierarchy();
    H->setlib(op->getDomainMap()->lib());

    // Stick the non-serializible data on the hierarchy.
    HierarchyUtils::AddNonSerializableDataToHierarchy(*mueLuFactory,*H, nonSerialList);

    // Set fine level operator
    H->GetLevel(0)->Set("A", op);

    // Set coordinates if available
    if (coords != Teuchos::null) {
      H->GetLevel(0)->Set("Coordinates", coords);
    }

    // Wrap nullspace if available, otherwise use constants
    if (nullspace == Teuchos::null) {
      int nPDE = MueLu::MasterList::getDefault<int>("number of equations");
      if (paramList.isSublist("Matrix")) {
        // Factory style parameter list
        const Teuchos::ParameterList& operatorList = paramList.sublist("Matrix");
        if (operatorList.isParameter("PDE equations"))
          nPDE = operatorList.get<int>("PDE equations");

      } else if (paramList.isParameter("number of equations")) {
        // Easy style parameter list
        nPDE = paramList.get<int>("number of equations");
      }

      nullspace = MultiVectorFactory::Build(op->getDomainMap(), nPDE);
      if (nPDE == 1) {
        nullspace->putScalar(Teuchos::ScalarTraits<Scalar>::one());

      } else {
        for (int i = 0; i < nPDE; i++) {
          Teuchos::ArrayRCP<Scalar> nsData = nullspace->getDataNonConst(i);
          for (int j = 0; j < nsData.size(); j++) {
            GlobalOrdinal GID = op->getDomainMap()->getGlobalElement(j) - op->getDomainMap()->getIndexBase();

            if ((GID-i) % nPDE == 0)
              nsData[j] = Teuchos::ScalarTraits<Scalar>::one();
          }
        }
      }
    }
    H->GetLevel(0)->Set("Nullspace", nullspace);


    mueLuFactory->SetupHierarchy(*H);

    tm->stop();
    tm->incrementNumCalls();

    if (H->GetVerbLevel() & Statistics0) {
      const bool alwaysWriteLocal = true;
      const bool writeGlobalStats = true;
      const bool writeZeroTimers  = false;
      const bool ignoreZeroTimers = true;
      const std::string filter    = timerName;
      Teuchos::TimeMonitor::summarize(op->getRowMap()->getComm().ptr(), std::cout, alwaysWriteLocal, writeGlobalStats,
                                      writeZeroTimers, Teuchos::Union, filter, ignoreZeroTimers);
    }

    tm->reset();

    return H;
  }

  TEUCHOS_UNIT_TEST(SaPFactory_kokkos, EpetraVsTpetra)
  {
#   include "MueLu_UseShortNames.hpp"
    MueLu::VerboseObject::SetDefaultOStream(Teuchos::rcpFromRef(out));

    out << "version: " << MueLu::Version() << std::endl;
    out << "Compare results of Epetra and Tpetra" << std::endl;
    out << "for 3 level AMG solver using smoothed aggregation with" << std::endl;
    out << "one SGS sweep on each multigrid level as pre- and postsmoother" << std::endl;

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

    typedef Teuchos::ScalarTraits<SC> STS;
    SC zero = STS::zero(), one = STS::one();

    Array<STS::magnitudeType> results(2);

    // run test only on 1 proc
    if(comm->getSize() == 1) {
      Xpetra::UnderlyingLib lib = Xpetra::UseEpetra;

      // run Epetra and Tpetra test
      for (int run = 0; run < 2; run++) { //TODO: create a subfunction instead or Tuple of UnderlyingLib
        if (run == 0) lib = Xpetra::UseEpetra;
        else          lib = Xpetra::UseTpetra;

        // generate problem
        LO maxLevels = 3;
        LO its       = 10;
        GO nEle      = 63;
        const RCP<const Map> map = MapFactory::Build(lib, nEle, 0, comm);
        Teuchos::ParameterList matrixParameters;
        matrixParameters.set("nx", nEle);

        RCP<Galeri::Xpetra::Problem<Map,CrsMatrixWrap,MultiVector> > Pr =
          Galeri::Xpetra::BuildProblem<SC,LO,GO,Map,CrsMatrixWrap,MultiVector>("Laplace1D", map, matrixParameters);
        RCP<Matrix> Op = Pr->BuildMatrix();

        // build nullspace
        RCP<MultiVector> nullSpace = MultiVectorFactory::Build(map,1);
        nullSpace->putScalar(one);
        Array<STS::magnitudeType> norms(1);
        nullSpace->norm1(norms);
        if (comm->getRank() == 0)
          out << "||NS|| = " << norms[0] << std::endl;

        // fill hierarchy
        RCP<Hierarchy> H = rcp( new Hierarchy() );
        H->setDefaultVerbLevel(Teuchos::VERB_HIGH);

        RCP<Level> Finest = H->GetLevel();
        Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
        Finest->Set("A",Op);                      // set fine level matrix
        Finest->Set("Nullspace",nullSpace);       // set null space information for finest level

        // define transfer operators
        RCP<CoupledAggregationFactory> CoupledAggFact = rcp(new CoupledAggregationFactory());
        CoupledAggFact->SetMinNodesPerAggregate(3);
        CoupledAggFact->SetMaxNeighAlreadySelected(0);
        CoupledAggFact->SetOrdering("natural");
        CoupledAggFact->SetPhase3AggCreation(0.5);

        RCP<TentativePFactory> Ptentfact = rcp(new TentativePFactory());
        RCP<SaPFactory>        Pfact = rcp( new SaPFactory());
        RCP<Factory>      Rfact = rcp( new TransPFactory() );
        RCP<RAPFactory>        Acfact = rcp( new RAPFactory() );
        H->SetMaxCoarseSize(1);

        // setup smoothers
        Teuchos::ParameterList smootherParamList;
        smootherParamList.set("relaxation: type", "Symmetric Gauss-Seidel");
        smootherParamList.set("relaxation: sweeps", (LO) 1);
        smootherParamList.set("relaxation: damping factor", (SC) 1.0);
        RCP<SmootherPrototype> smooProto = rcp( new TrilinosSmoother("RELAXATION", smootherParamList) );
        RCP<SmootherFactory> SmooFact = rcp( new SmootherFactory(smooProto) );
        Acfact->setVerbLevel(Teuchos::VERB_HIGH);

        RCP<SmootherFactory> coarseSolveFact = rcp(new SmootherFactory(smooProto, Teuchos::null));

        FactoryManager M;
        M.SetFactory("P", Pfact);
        M.SetFactory("R", Rfact);
        M.SetFactory("A", Acfact);
        M.SetFactory("Ptent", Ptentfact);
        M.SetFactory("Aggregates", CoupledAggFact);
        M.SetFactory("Smoother", SmooFact);
        M.SetFactory("CoarseSolver", coarseSolveFact);

        H->Setup(M, 0, maxLevels);

        // test some basic multigrid data
        RCP<Level> coarseLevel = H->GetLevel(1);
        TEST_EQUALITY(coarseLevel->IsRequested("A",MueLu::NoFactory::get()), false);
        TEST_EQUALITY(coarseLevel->IsRequested("P",MueLu::NoFactory::get()), false);
        TEST_EQUALITY(coarseLevel->IsRequested("PreSmoother",MueLu::NoFactory::get()), false);
        TEST_EQUALITY(coarseLevel->IsRequested("PostSmoother",MueLu::NoFactory::get()), false);
        TEST_EQUALITY(coarseLevel->IsRequested("R",MueLu::NoFactory::get()), false);
        TEST_EQUALITY(coarseLevel->IsAvailable("A",MueLu::NoFactory::get()), true);
        TEST_EQUALITY(coarseLevel->IsAvailable("P",MueLu::NoFactory::get()), true);
        TEST_EQUALITY(coarseLevel->IsAvailable("PreSmoother",MueLu::NoFactory::get()), true);
//.........这里部分代码省略.........

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

  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
  Teuchos::RCP<const Teuchos::Comm<int> > commptr =
    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;
  Teuchos::RCP<std::ostream> outStream;
  Teuchos::oblackholestream bhs; // outputs nothing
  if (iprint > 0 && commptr->getRank() == 0)
    outStream = Teuchos::rcp(&std::cout, false);
  else
    outStream = Teuchos::rcp(&bhs, false);

  int errorFlag  = 0;

  try {
    /**********************************************************************************************/
    /************************* CONSTRUCT SOL COMPONENTS *******************************************/
    /**********************************************************************************************/
    // Set random seed
    srand(123456789);
    // Build samplers
    size_t dimension = 1;

    // Initialize distribution
    Teuchos::RCP<ROL::Distribution<RealT> > dist;
    std::vector<Teuchos::RCP<ROL::Distribution<RealT> > > distVec(dimension);
    Teuchos::ParameterList Dlist;
    Dlist.sublist("SOL").sublist("Distribution").set("Name","Beta");
    RealT alpha = 1., beta = 4.;
    // Fill moment vector and initial guess
    for (size_t d = 0; d < dimension; d++) {
      // Build distribution for dimension d
      alpha++; beta++;
      Dlist.sublist("SOL").sublist("Distribution").sublist("Beta").set("Shape 1",alpha);
      Dlist.sublist("SOL").sublist("Distribution").sublist("Beta").set("Shape 2",beta);
      dist = ROL::DistributionFactory<RealT>(Dlist);
      distVec[d] = ROL::DistributionFactory<RealT>(Dlist);
    }

    // Get ROL parameterlist
    std::string filename = "input_04.xml";
    Teuchos::RCP<Teuchos::ParameterList> parlist = Teuchos::rcp( new Teuchos::ParameterList() );
    Teuchos::updateParametersFromXmlFile( filename, parlist.ptr() );

    Teuchos::ParameterList &list = parlist->sublist("SOL").sublist("Sample Generator").sublist("SROM");
    Teuchos::Array<int> moments = Teuchos::getArrayFromStringParameter<int>(list,"Moments");
    size_t numMoments = static_cast<size_t>(moments.size());

    std::clock_t timer = std::clock();
    Teuchos::RCP<ROL::BatchManager<RealT> > bman =
      Teuchos::rcp(new ROL::TeuchosBatchManager<RealT,int>(commptr));
    Teuchos::RCP<ROL::SampleGenerator<RealT> > sampler =
      Teuchos::rcp(new ROL::SROMGenerator<RealT>(*parlist,bman,distVec));
    *outStream << std::endl << "Sample Time: "
               << (std::clock()-timer)/(RealT)CLOCKS_PER_SEC << " seconds"
               << std::endl;

    RealT val = 0., error = 0., data = 0., sum = 0.;
    *outStream << std::endl;
    *outStream << std::scientific << std::setprecision(11);
    *outStream << std::right << std::setw(20) << "Computed Moment"
                             << std::setw(20) << "True Moment"
                             << std::setw(20) << "Relative Error"
                             << std::endl;
    for (size_t m = 0; m < numMoments; m++) {
      for (size_t d = 0; d < dimension; d++) {
        val = 0.; data = distVec[d]->moment(moments[m]);
        for (size_t k = 0; k < (size_t)sampler->numMySamples(); k++) {
          val += sampler->getMyWeight(k)*std::pow((sampler->getMyPoint(k))[d],moments[m]);
        }
        bman->sumAll(&val,&sum,1);
        error = std::abs(sum-data)/std::abs(data);
        if ( error > 1.e-1 ) {
          errorFlag++;
        }
        *outStream << std::right << std::setw(20) << sum
                                 << std::setw(20) << data
                                 << std::setw(20) << error
                                 << std::endl;
      }
    }
    *outStream << std::endl;

//    std::ofstream file;
//    std::stringstream name;
//    name << "samples." << commptr->getRank() << ".txt";
//    file.open(name.str().c_str()); 
//    for (size_t k = 0; k < (size_t)sampler->numMySamples(); k++) {
//      for (size_t d = 0; d < dimension; d++) {
//        file << std::setprecision(std::numeric_limits<RealT>::digits10)
//             << std::scientific
//             << (sampler->getMyPoint(k))[d];
//        file << "  ";
//      }
//      file << std::setprecision(std::numeric_limits<RealT>::digits10)
//           << std::scientific
//           << sampler->getMyWeight(k) << std::endl;
//.........这里部分代码省略.........

bool Options::setOptions(Teuchos::ParameterList& nlParams)
{

  // Set status tests if not already set
  if( Teuchos::is_null(testCombo) )
  {
    // Check for MaxIters option
    int maxIters;
#if  (PETSC_VERSION_MAJOR >= 3) || (PETSC_VERSION_MINOR >= 5)
    PetscBool lflg;
#else
    PetscTruth lflg;  // Needed to permit two ways of specification
#endif
    ierr = PetscOptionsGetInt(PETSC_NULL,"-snes_max_it", &maxIters, &flg);CHKERRQ(ierr);
    ierr = PetscOptionsGetInt(PETSC_NULL,"-nox_conv_maxiters", &maxIters, &lflg);CHKERRQ(ierr);
    if(flg || lflg)
    {
      testMaxIters = Teuchos::rcp( new NOX::StatusTest::MaxIters(maxIters) );
      if( Teuchos::is_null(testCombo) )
        testCombo = Teuchos::rcp( new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, testMaxIters) );
      else
        testCombo->addStatusTest(testMaxIters);
    }

    // Check for (absolute) residual norm (L2-norm) tolerance
    double absResNorm;
    PetscReal petscVal;
    ierr = PetscOptionsGetReal(PETSC_NULL,"-snes_atol", &petscVal, &flg);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(PETSC_NULL,"-nox_conv_abs_res", &petscVal, &lflg);CHKERRQ(ierr);
    if(flg || lflg)
    {
      absResNorm = (double) petscVal;
      testNormF = Teuchos::rcp( new NOX::StatusTest::NormF(absResNorm) );
      if( Teuchos::is_null(testCombo) )
        testCombo = Teuchos::rcp( new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, testNormF) );
      else
        testCombo->addStatusTest(testNormF);
    }

    // Check for update norm (L2-norm) tolerance
    double absUpdateNorm;
    ierr = PetscOptionsGetReal(PETSC_NULL,"-snes_stol", &petscVal, &flg);CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(PETSC_NULL,"-nox_conv_update", &petscVal, &lflg);CHKERRQ(ierr);
    if(flg || lflg)
    {
      absUpdateNorm = (double) petscVal;
      testNormUpdate = Teuchos::rcp( new NOX::StatusTest::NormUpdate(absUpdateNorm) );
      if( Teuchos::is_null(testCombo) )
        testCombo = Teuchos::rcp( new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, testNormUpdate) );
      else
        testCombo->addStatusTest(testNormUpdate);
    }

    // Finally, provide a default test if none specified
    if( Teuchos::is_null(testCombo) ) // No tests specified by the uesr
    {
      assert( Teuchos::is_null(testMaxIters) );
      testMaxIters = Teuchos::rcp( new NOX::StatusTest::MaxIters(20) );
      assert( Teuchos::is_null(testNormF) );
      testNormF = Teuchos::rcp( new NOX::StatusTest::NormF(1.e-12) );
      testCombo = Teuchos::rcp( new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR, testMaxIters, testNormF) );
    }


  } // End of StatusTest construction


  // Allow solution-type to be specified
  ierr = PetscOptionsHasName(PETSC_NULL,"-nox_trustregion_based",&flg);
         CHKERRQ(ierr);
  if(flg)
    nlParams.set("Nonlinear Solver", "Trust Region Based");
  else // default
    // This is done to allow PetscOptions to register that this option was used
    ierr = PetscOptionsHasName(PETSC_NULL,"-nox_linesearch_based",&flg);
           CHKERRQ(ierr);
    nlParams.set("Nonlinear Solver", "Line Search Based");

  // Now allow linesearch type to be specified
  Teuchos::ParameterList& searchParams = nlParams.sublist("Line Search");
  ierr = PetscOptionsGetString(PETSC_NULL,"-nox_linesearch_type",
               optionString, maxStringLength, &flg);CHKERRQ(ierr);
  if(flg)
  {
    if( !strcmp(optionString, "full_step") )
      searchParams.set("Method", "Full Step");
    if( !strcmp(optionString, "polynomial") )
      searchParams.set("Method", "Polynomial");
    if( !strcmp(optionString, "backtrack") )
      searchParams.set("Method", "Backtrack");
    if( !strcmp(optionString, "more_thuente") )
      searchParams.set("Method", "More'-Thuente");
#ifdef WITH_PRERELEASE
    if( !strcmp(optionString, "nonlinearcg") )
      searchParams.set("Method", "NonlinearCG");
#endif
  }
  else // default
    searchParams.set("Method", "Full Step");

//.........这里部分代码省略.........