C++ Teuchos::parameterList方法代码示例

Teuchos::RCP<SolverManagerBaseType>
makeSolverManagerTmpl (const Teuchos::RCP<Teuchos::ParameterList>& params)
{
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;

  RCP<SolverManagerType> solver = rcp (new SolverManagerType);

  // Some solvers may not like to get a null ParameterList.  If params
  // is null, replace it with an empty parameter list.  The solver
  // will fill in default parameters for that case.  Use the name of
  // the solver's default parameters to name the new empty list.
  RCP<ParameterList> pl;
  if (params.is_null()) {
    pl = parameterList (solver->getValidParameters ()->name ());
  } else {
    pl = params;
  }
  TEUCHOS_TEST_FOR_EXCEPTION(
    pl.is_null(), std::logic_error,
    "Belos::SolverFactory: ParameterList to pass to solver is null.  This "
    "should never happen.  Please report this bug to the Belos developers.");
  solver->setParameters (pl);
  return solver;
}

void DefaultLinearSolverBuilder::readParameters( std::ostream *out )
{
  using Teuchos::parameterList;
  using Teuchos::ptr;
  using Teuchos::updateParametersFromXmlFile;
  using Teuchos::updateParametersFromXmlString;
  using std::endl;
  
  if (!paramList_.get()) {
    paramList_ = parameterList("DefaultLinearSolverBuilder");
  }
  if (paramsXmlFileName().length()) {
    if (out) {
      *out << endl << "Reading parameters from XML file \"" 
           << paramsXmlFileName() << "\" ..." << endl;
    }
    updateParametersFromXmlFile (paramsXmlFileName (), paramList_.ptr());
  }
  if (extraParamsXmlString().length()) {
    if (out) {
      *out << endl << "Appending extra parameters from the XML string \""
           << extraParamsXmlString() << "\" ..." << endl;
    }
    updateParametersFromXmlString (extraParamsXmlString (), paramList_.ptr());
  }
  setParameterList(paramList_);
}

      /// \brief Return a valid parameter list for verifying Tsqr.
      ///
      /// Call this once to get a valid parameter list with all the
      /// defaults filled in.  This list is valid for all the Scalar
      /// types which TsqrVerifierCaller::run tests.
      Teuchos::RCP<const Teuchos::ParameterList>
      getValidParameterList () const
      {
        using Teuchos::ParameterList;
        using Teuchos::parameterList;
        using Teuchos::RCP;

        RCP<ParameterList> plist = parameterList ("FullTsqrVerifier");

        const size_t cacheSizeHint = 0;
        const int numCores = 1;
        const ordinal_type numRowsLocal = 100;
        const ordinal_type numCols = 10;
        const bool contiguousCacheBlocks = false;
        const bool testFactorExplicit = true;
        const bool testRankRevealing = true;
        const bool printFieldNames = true;
        const bool printResults = true;
        const bool failIfInaccurate = true;
        const bool debug = false;

        // Parameters for configuring Tsqr itself.
        plist->set ("cacheSizeHint", cacheSizeHint,
                    "Cache size hint in bytes.  "
                    "Zero means TSQR picks a reasonable default.");
        plist->set ("numCores", numCores,
                    "Number of partition(s) to use for TbbTsqr (if "
                    "applicable).  Must be a positive integer.");

        // Parameters for testing Tsqr.
        plist->set ("numRowsLocal", numRowsLocal,
                    "Number of rows per (MPI) process in the test matrix.  "
                    "Must be >= the number of columns.");
        plist->set ("numCols", numCols,
                    "Number of columns in the test matrix.");
        plist->set ("contiguousCacheBlocks", contiguousCacheBlocks,
                    "Whether to test the factorization with contiguously "
                    "stored cache blocks.");
        plist->set ("testFactorExplicit", testFactorExplicit,
                    "Whether to test TSQR's factorExplicit() (a hopefully "
                    "faster path than calling factor() and explicit_Q() in "
                    "sequence).");
        plist->set ("testRankRevealing", testRankRevealing,
                    "Whether to test TSQR's rank-revealing capability.");
        plist->set ("printFieldNames", printFieldNames,
                    "Whether to print field names (this is only done once, "
                    "for all Scalar types tested).");
        plist->set ("printResults", printResults,
                    "Whether to print test results.");
        plist->set ("failIfInaccurate", failIfInaccurate,
                    "Whether to fail the test if the factorization "
                    "is not sufficiently accurate.");
        plist->set ("debug", debug,
                    "Whether to print debugging output.");
        return plist;
      }

    void
    setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
    {
        using Teuchos::ParameterList;
        using Teuchos::parameterList;
        using Teuchos::RCP;
        using Teuchos::sublist;

        RCP<ParameterList> params = plist.is_null() ?
                                    parameterList (*getValidParameters ()) : plist;
        nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
        distTsqr_->setParameterList (sublist (params, "DistTsqr"));

        this->setMyParamList (params);
    }

    Teuchos::RCP<const Teuchos::ParameterList>
    getValidParameters () const
    {
        using Teuchos::RCP;
        using Teuchos::rcp;
        using Teuchos::ParameterList;
        using Teuchos::parameterList;

        if (defaultParams_.is_null()) {
            RCP<ParameterList> params = parameterList ("TSQR implementation");
            params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
            params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
            defaultParams_ = params;
        }
        return defaultParams_;
    }

  Teuchos::RCP<const Teuchos::ParameterList>
  MinresSolMgr<ScalarType, MV, OP>::defaultParameters()
  {
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::rcpFromRef;
    using Teuchos::EnhancedNumberValidator;
    typedef MagnitudeType MT;
    typedef Teuchos::ScalarTraits<MT> MST;

    // List of parameters accepted by MINRES, and their default values.
    RCP<ParameterList> pl = parameterList ("MINRES");

    pl->set ("Convergence Tolerance", MST::squareroot (MST::eps()),
	     "Relative residual tolerance that needs to be achieved by "
	     "the iterative solver, in order for the linear system to be "
	     "declared converged.",
	     rcp (new EnhancedNumberValidator<MT> (MST::zero(), MST::rmax())));
    pl->set ("Maximum Iterations", static_cast<int>(1000),
	     "Maximum number of iterations allowed for each right-hand "
	     "side solved.",
	     rcp (new EnhancedNumberValidator<int> (0, INT_MAX)));
    pl->set ("Block Size", static_cast<int>(1),
	     "Number of vectors in each block.  WARNING: The current "
	     "implementation of MINRES only accepts a block size of 1, "
	     "since it can only solve for 1 right-hand side at a time.",
	     rcp (new EnhancedNumberValidator<int> (1, 1)));
    pl->set ("Verbosity", (int) Belos::Errors,
	     "The type(s) of solver information that should "
	     "be written to the output stream.");
    pl->set ("Output Style", (int) Belos::General,
	     "What style is used for the solver information written "
	     "to the output stream.");
    pl->set ("Output Frequency", static_cast<int>(-1),
	     "How often (in terms of number of iterations) intermediate "
	     "convergence information should be written to the output stream."
	     "  -1 means never.");
    pl->set ("Output Stream", rcpFromRef(std::cout),
	     "A reference-counted pointer to the output stream where all "
	     "solver output is sent.  The output stream defaults to stdout.");
    pl->set ("Timer Label", std::string("Belos"),
	     "The string to use as a prefix for the timer labels.");
    return pl;
  }

 RCP<typename Kokkos::DefaultKernels<float,int,Node>::SparseOps::template bind_scalar<float>::other_type>
 gen_prob(RCP<Node> node, int N, size_t &totalNNZ)
 {
   typedef typename Kokkos::DefaultKernels<float,int,Node>::SparseOps   DSM;
   typedef typename DSM::template bind_scalar<float>::other_type       fDSM;
   typedef typename fDSM::template graph<int,Node>::graph_type                  GRPH;
   typedef typename fDSM::template matrix<float,int,Node>::matrix_type           MAT;
   // generate symmetric tridiagonal matrix
   RCP<GRPH> G = rcp(new GRPH(N,N,node,null));
   RCP<MAT>  A= rcp(new MAT(G,null));
   // allocate buffers for offsets, indices and values
   totalNNZ = 3*N - 2;
   ArrayRCP<size_t> offsets(N+1);
   ArrayRCP<int>    inds(totalNNZ);
   ArrayRCP<float>  vals(totalNNZ);
   {
     size_t NNZsofar = 0;
     offsets[0] = NNZsofar;
     inds[NNZsofar] = 0; inds[NNZsofar+1] =  1;
     vals[NNZsofar] = 2; vals[NNZsofar+1] = -1;
     NNZsofar += 2;
     for (int i=1; i != N-1; ++i) {
       offsets[i] = NNZsofar;
       inds[NNZsofar] = i-1; inds[NNZsofar+1] = i; inds[NNZsofar+2] = i+1;
       vals[NNZsofar] =  -1; vals[NNZsofar+1] = 2; vals[NNZsofar+2] =  -1;
       NNZsofar += 3;
     }
     offsets[N-1] = NNZsofar;
     inds[NNZsofar] = N-2; inds[NNZsofar+1] = N-1;
     vals[NNZsofar] =  -1; vals[NNZsofar+1] = 2;
     NNZsofar += 2;
     offsets[N]   = NNZsofar;
     TEUCHOS_TEST_FOR_EXCEPT(NNZsofar != totalNNZ);
   }
   G->setStructure(offsets, inds);
   offsets = Teuchos::null;
   inds    = Teuchos::null;
   A->setValues(vals);
   vals    = Teuchos::null;
   fDSM::finalizeGraphAndMatrix(Teuchos::UNDEF_TRI,Teuchos::NON_UNIT_DIAG,*G,*A,parameterList());
   RCP<fDSM> dsm = rcp(new fDSM(node));
   dsm->setGraphAndMatrix(G,A);
   return dsm;
 }

  void run(Teuchos::ParameterList &myMachPL, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, const Teuchos::RCP<Node> &node) 
  {
    using std::pair;
    using std::make_pair;
    using std::plus;
    using std::endl;
    using Teuchos::null;
    using Teuchos::RCP;
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using TpetraExamples::make_pair_op;
    using Tpetra::RTI::reductionGlob;
    using Tpetra::RTI::ZeroOp;
    using Tpetra::RTI::binary_pre_transform_reduce;
    using Tpetra::RTI::binary_transform;

    // Static types
    typedef typename MPStack::type   S;
    typedef int                     LO;
    typedef int                     GO;
    typedef Tpetra::Map<LO,GO,Node>               Map;
    typedef Tpetra::CrsMatrix<S,LO,GO,Node> CrsMatrix;
    typedef Tpetra::Vector<S,LO,GO,Node>       Vector;

    *out << "Running test with Node==" << Teuchos::typeName(*node) << " on rank " << comm->getRank() << "/" << comm->getSize() << std::endl;

    // read the matrix
    RCP<CrsMatrix> A;
    RCP<const Map> rowMap = null;
    RCP<ParameterList> fillParams = parameterList();
    fillParams->set("Preserve Local Graph",true);
    // must preserve the local graph in order to do convert() calls later
    Tpetra::Utils::readHBMatrix(matrixFile,comm,node,A,rowMap,fillParams);
    rowMap = A->getRowMap();

    // init the solver stack
    TpetraExamples::RFPCGInit<S,LO,GO,Node> init(A);
    RCP<ParameterList> db = Tpetra::Ext::initStackDB<MPStack>(*params,init);

    testPassed = true;

    // choose a solution, compute a right-hand-side
    auto x = Tpetra::createVector<S>(rowMap),
         b = Tpetra::createVector<S>(rowMap);
    x->randomize();
    A->apply(*x,*b);
    {
      // init the rhs
      auto bx = db->get<RCP<Vector>>("bx");
      binary_transform( *bx, *b, [](S, S bi) {return bi;}); // bx = b
    }

    // call the solve
    TpetraExamples::recursiveFPCG<MPStack,LO,GO,Node>(out,*db);

    // check that residual is as requested
    {
      auto xhat = db->get<RCP<Vector>>("bx"),
           bhat = Tpetra::createVector<S>(rowMap);
      A->apply(*xhat,*bhat);
      // compute bhat-b, while simultaneously computing |bhat-b|^2 and |b|^2
      auto nrms = binary_pre_transform_reduce(*bhat, *b, 
                                              reductionGlob<ZeroOp<pair<S,S>>>( 
                                                [](S bhati, S bi){ return bi-bhati;}, // bhati = bi-bhat
                                                [](S bhati, S bi){ return make_pair(bhati*bhati, bi*bi); },
                                                make_pair_op<S,S>(plus<S>())) );
      const S enrm = Teuchos::ScalarTraits<S>::squareroot(nrms.first),
              bnrm = Teuchos::ScalarTraits<S>::squareroot(nrms.second);
      // check that residual is as requested
      *out << "|b - A*x|/|b|: " << enrm / bnrm << endl;
      const double tolerance = db->get<double>("tolerance");
      if (MPStack::bottom) {
        // give a little slack
        if (enrm / bnrm > 5*tolerance) testPassed = false;
      }
      else {
        if (enrm / bnrm > tolerance) testPassed = false;
      }
    }

    // 
    // solve again, with the unfused version, just for timings purposes
    if (unfusedTest) 
    {
      // init the rhs
      auto bx = db->get<RCP<Vector>>("bx");
      binary_transform( *bx, *b, [](S, S bi) {return bi;}); // bx = b
      // call the solve
      TpetraExamples::recursiveFPCGUnfused<MPStack,LO,GO,Node>(out,*db);
      //
      // test the result
      auto xhat = db->get<RCP<Vector>>("bx"),
           bhat = Tpetra::createVector<S>(rowMap);
      A->apply(*xhat,*bhat);
      // compute bhat-b, while simultaneously computing |bhat-b|^2 and |b|^2
      auto nrms = binary_pre_transform_reduce(*bhat, *b, 
                                              reductionGlob<ZeroOp<pair<S,S>>>( 
                                                [](S bhati, S bi){ return bi-bhati;}, // bhati = bi-bhat
                                                [](S bhati, S bi){ return make_pair(bhati*bhati, bi*bi); },
                                                make_pair_op<S,S>(plus<S>())) );
//.........这里部分代码省略.........

/// \fn main
/// \brief Benchmark driver for (Mat)OrthoManager subclasses
  int
main (int argc, char *argv[])
{
  using Belos::OrthoManager;
  using Belos::OrthoManagerFactory;
  using Belos::OutputManager;
  using Teuchos::CommandLineProcessor;
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;

  Tpetra::ScopeGuard tpetraScope (&argc, &argv);

  bool success = false;
  bool verbose = false; // Verbosity of output
  try {
    RCP<const Teuchos::Comm<int> > pComm = Tpetra::getDefaultComm();
    // This factory object knows how to make a (Mat)OrthoManager
    // subclass, given a name for the subclass.  The name is not the
    // same as the class' syntactic name: e.g., "TSQR" is the name of
    // TsqrOrthoManager.
    OrthoManagerFactory<scalar_type, MV, OP> factory;

    // The name of the (Mat)OrthoManager subclass to instantiate.
    std::string orthoManName (factory.defaultName());

    // For SimpleOrthoManager: the normalization method to use.  Valid
    // values: "MGS", "CGS".
    std::string normalization ("CGS");

    // Name of the Harwell-Boeing sparse matrix file from which to read
    // the inner product operator matrix.  If name is "" or not provided
    // at the command line, use the standard Euclidean inner product.
    std::string filename;

    bool debug = false;   // Whether to print debugging-level output
    // Whether or not to run the benchmark.  If false, we let this
    // "test" pass trivially.
    bool benchmark = false;

    // Whether to display benchmark results compactly (in a CSV format),
    // or in a human-readable table.
    bool displayResultsCompactly = false;

    // Default _local_ (per MPI process) number of rows.  This will
    // change if a sparse matrix is loaded in as an inner product
    // operator.  Regardless, the number of rows per MPI process must be
    // no less than numCols*numBlocks in order for TSQR to work.  To
    // ensure that the test always passes with default parameters, we
    // scale by the number of processes.  The default value below may be
    // changed by a command-line parameter with a corresponding name.
    int numRowsPerProcess = 100;

    // The OrthoManager is benchmarked with numBlocks multivectors of
    // width numCols each, for numTrials trials.  The values below are
    // defaults and may be changed by the corresponding command-line
    // arguments.
    int numCols = 10;
    int numBlocks = 5;
    int numTrials = 3;

    CommandLineProcessor cmdp (false, true);
    cmdp.setOption ("benchmark", "nobenchmark", &benchmark,
        "Whether to run the benchmark.  If not, this \"test\" "
        "passes trivially.");
    cmdp.setOption ("verbose", "quiet", &verbose,
        "Print messages and results.");
    cmdp.setOption ("debug", "nodebug", &debug,
        "Print debugging information.");
    cmdp.setOption ("compact", "human", &displayResultsCompactly,
        "Whether to display benchmark results compactly (in a "
        "CSV format), or in a human-readable table.");
    cmdp.setOption ("filename", &filename,
        "Filename of a Harwell-Boeing sparse matrix, used as the "
        "inner product operator by the orthogonalization manager."
        "  If not provided, no matrix is read and the Euclidean "
        "inner product is used.");
    {
      std::ostringstream os;
      const int numValid = factory.numOrthoManagers();
      const bool plural = numValid > 1 || numValid == 0;

      os << "OrthoManager subclass to benchmark.  There ";
      os << (plural ? "are " : "is ") << numValid << (plural ? "s: " : ": ");
      factory.printValidNames (os);
      os << ".  If none is provided, the test trivially passes.";
      cmdp.setOption ("ortho", &orthoManName, os.str().c_str());
    }
    cmdp.setOption ("normalization", &normalization,
        "For SimpleOrthoManager (--ortho=Simple): the normalization "
        "method to use.  Valid values: \"MGS\", \"CGS\".");
    cmdp.setOption ("numRowsPerProcess", &numRowsPerProcess,
        "Number of rows per MPI process in the test multivectors.  "
        "If an input matrix is given, this value is ignored, since "
        "the vectors must be commensurate with the dimensions of "
        "the matrix.");
    cmdp.setOption ("numCols", &numCols,
//.........这里部分代码省略.........

SolveStatus<Scalar>
BelosLinearOpWithSolve<Scalar>::solveImpl(
  const EOpTransp M_trans,
  const MultiVectorBase<Scalar> &B,
  const Ptr<MultiVectorBase<Scalar> > &X,
  const Ptr<const SolveCriteria<Scalar> > solveCriteria
  ) const
{

  THYRA_FUNC_TIME_MONITOR("Stratimikos: BelosLOWS");

  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using Teuchos::rcpFromPtr;
  using Teuchos::FancyOStream;
  using Teuchos::OSTab;
  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::describe;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef typename ST::magnitudeType ScalarMag;
  Teuchos::Time totalTimer(""), timer("");
  totalTimer.start(true);

  assertSolveSupports(*this, M_trans, solveCriteria);
  // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
  // solve(...).

  const RCP<FancyOStream> out = this->getOStream();
  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  OSTab tab = this->getOSTab();
  if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)) {
    *out << "\nStarting iterations with Belos:\n";
    OSTab tab2(out);
    *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
    *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
    *out << "With #Eqns="<<B.range()->dim()<<", #RHSs="<<B.domain()->dim()<<" ...\n";
  }

  //
  // Set RHS and LHS
  //

  bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
  TEUCHOS_TEST_FOR_EXCEPTION(
    ret == false, CatastrophicSolveFailure
    ,"Error, the Belos::LinearProblem could not be set for the current solve!"
    );

  //
  // Set the solution criteria
  //

  // Parameter list for the current solve.
  const RCP<ParameterList> tmpPL = Teuchos::parameterList();

  // The solver's valid parameter list.
  RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();

  SolveMeasureType solveMeasureType;
  RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
  if (nonnull(solveCriteria)) {
    solveMeasureType = solveCriteria->solveMeasureType;
    const ScalarMag requestedTol = solveCriteria->requestedTol;
    if (solveMeasureType.useDefault()) {
      tmpPL->set("Convergence Tolerance", defaultTol_);
    }
    else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
      if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
        tmpPL->set("Convergence Tolerance", requestedTol);
      }
      else {
        tmpPL->set("Convergence Tolerance", defaultTol_);
      }
      setResidualScalingType (tmpPL, validPL, "Norm of RHS");
    }
    else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
      if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
        tmpPL->set("Convergence Tolerance", requestedTol);
      }
      else {
        tmpPL->set("Convergence Tolerance", defaultTol_);
      }
      setResidualScalingType (tmpPL, validPL, "Norm of Initial Residual");
    }
    else {
      // Set the most generic (and inefficient) solve criteria
      generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
        *solveCriteria, convergenceTestFrequency_);
      // Set the verbosity level (one level down)
      generalSolveCriteriaBelosStatusTest->setOStream(out);
      generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
      // Set the default convergence tolerance to always converged to allow
      // the above status test to control things.
      tmpPL->set("Convergence Tolerance", 1.0);
    }
    // maximum iterations
    if (nonnull(solveCriteria->extraParameters)) {
      if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,"Maximum Iterations")) {
        tmpPL->set("Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,"Maximum Iterations"));
//.........这里部分代码省略.........

    void
    benchmarkKokkosNodeTsqr (const Teuchos::RCP<NodeType>& node,
                             const int numTrials,
                             const Ordinal numRows,
                             const Ordinal numCols,
                             const int numPartitions,
                             const size_t cacheSizeHint,
                             const bool contiguousCacheBlocks,
                             const bool printFieldNames,
                             const bool humanReadable)
    {
      using Teuchos::ParameterList;
      using Teuchos::parameterList;
      using Teuchos::RCP;
      using Teuchos::TypeNameTraits;
      using std::cerr;
      using std::cout;
      using std::endl;
      typedef TSQR::KokkosNodeTsqr<Ordinal, Scalar, NodeType> node_tsqr_type;
      typedef typename node_tsqr_type::FactorOutput factor_output_type;
      typedef Teuchos::ScalarTraits<Scalar> STS;
      // typedef typename STS::magnitudeType magnitude_type;
      typedef Teuchos::Time timer_type;
      typedef Matrix<Ordinal, Scalar> matrix_type;

      const std::string scalarTypeName = TypeNameTraits<Scalar>::name();

      // Pseudorandom normal(0,1) generator.  Default seed is OK,
      // because this is a benchmark, not an accuracy test.
      TSQR::Random::NormalGenerator<Ordinal, Scalar> gen;

      // Set up TSQR implementation.
      RCP<ParameterList> params = parameterList ("Intranode TSQR");
      params->set ("Cache Size Hint", cacheSizeHint);
      params->set ("Num Tasks", numPartitions);
      node_tsqr_type actor (params);
      actor.setNode (node);

      // Allocate space for test problem.
      matrix_type A (numRows, numCols);
      matrix_type A_copy (numRows, numCols);
      matrix_type Q (numRows, numCols);
      matrix_type R (numCols, numCols);

      // Fill R with zeros, since the factorization may not overwrite
      // the strict lower triangle of R.
      R.fill (STS::zero());

      // Create a test problem
      nodeTestProblem (gen, numRows, numCols, A.get(), A.lda(), false);

      // Copy A into A_copy, since TSQR overwrites the input.  If
      // specified, rearrange the data in A_copy so that the data in
      // each cache block is contiguously stored.
      if (contiguousCacheBlocks) {
        actor.cache_block (numRows, numCols, A_copy.get(), A.get(), A.lda());
      } else {
        deep_copy (A_copy, A);
      }

      // Do a few timing runs and throw away the results, just to warm
      // up any libraries that do autotuning.
      const int numWarmupRuns = 5;
      for (int warmupRun = 0; warmupRun < numWarmupRuns; ++warmupRun) {
        // Factor the matrix in-place in A_copy, and extract the
        // resulting R factor into R.
        factor_output_type factor_output =
          actor.factor (numRows, numCols, A_copy.get(), A_copy.lda(),
                        R.get(), R.lda(), contiguousCacheBlocks);
        // Compute the explicit Q factor (which was stored
        // implicitly in A_copy and factor_output) and store in Q.
        // We don't need to un-cache-block the output, because we
        // aren't verifying it here.
        actor.explicit_Q (numRows, numCols, A_copy.get(), A_copy.lda(),
                          factor_output, numCols, Q.get(), Q.lda(),
                          contiguousCacheBlocks);
      }

      // Benchmark intranode TSQR for numTrials trials.
      //
      // Name of timer doesn't matter here; we only need the timing.
      timer_type timer("KokkosNodeTsqr");
      timer.start();
      for (int trialNum = 0; trialNum < numTrials; ++trialNum) {
        // Factor the matrix in-place in A_copy, and extract the
        // resulting R factor into R.
        factor_output_type factor_output =
          actor.factor (numRows, numCols, A_copy.get(), A_copy.lda(),
                        R.get(), R.lda(), contiguousCacheBlocks);
        // Compute the explicit Q factor (which was stored
        // implicitly in A_copy and factor_output) and store in Q.
        // We don't need to un-cache-block the output, because we
        // aren't verifying it here.
        actor.explicit_Q (numRows, numCols, A_copy.get(), A_copy.lda(),
                          factor_output, numCols, Q.get(), Q.lda(),
                          contiguousCacheBlocks);
      }
      const double timing = timer.stop();

      // Print the results
//.........这里部分代码省略.........

    void
    verifyKokkosNodeTsqr (const Teuchos::RCP<NodeType>& node,
                          TSQR::Random::NormalGenerator<Ordinal, Scalar>& gen,
                          const Ordinal numRows,
                          const Ordinal numCols,
                          const int numPartitions,
                          const size_t cacheSizeHint,
                          const bool contiguousCacheBlocks,
                          const bool printFieldNames,
                          const bool humanReadable,
                          const bool debug)
    {
      using Teuchos::ParameterList;
      using Teuchos::parameterList;
      using Teuchos::RCP;
      using Teuchos::TypeNameTraits;
      using std::cerr;
      using std::cout;
      using std::endl;
      typedef TSQR::KokkosNodeTsqr<Ordinal, Scalar, NodeType> node_tsqr_type;
      typedef typename node_tsqr_type::FactorOutput factor_output_type;
      typedef Teuchos::ScalarTraits<Scalar> STS;
      typedef typename STS::magnitudeType magnitude_type;
      // typedef Teuchos::Time timer_type;
      typedef Matrix<Ordinal, Scalar> matrix_type;
      typedef MatView<Ordinal, Scalar> mat_view_type;

      const std::string scalarTypeName = TypeNameTraits<Scalar>::name();

      // Set up TSQR implementation.
      RCP<ParameterList> params = parameterList ("Intranode TSQR");
      params->set ("Cache Size Hint", cacheSizeHint);
      params->set ("Num Tasks", numPartitions);
      node_tsqr_type actor (params);
      actor.setNode (node);
      if (debug)
        {
          cerr << actor.description() << endl;
          if (contiguousCacheBlocks)
            cerr << "-- Test with contiguous cache blocks" << endl;
        }

      // Allocate space for test problem.
      matrix_type A (numRows, numCols);
      matrix_type A_copy (numRows, numCols);
      matrix_type Q (numRows, numCols);
      matrix_type R (numCols, numCols);
      if (std::numeric_limits<Scalar>::has_quiet_NaN)
        {
          A.fill (std::numeric_limits<Scalar>::quiet_NaN());
          A_copy.fill (std::numeric_limits<Scalar>::quiet_NaN());
          Q.fill (std::numeric_limits<Scalar>::quiet_NaN());
          R.fill (std::numeric_limits<Scalar>::quiet_NaN());
        }
      else
        {
          A.fill (STS::zero());
          A_copy.fill (STS::zero());
          Q.fill (STS::zero());
          R.fill (STS::zero());
        }
      const Ordinal lda = numRows;
      const Ordinal ldq = numRows;
      const Ordinal ldr = numCols;

      // Create a test problem
      nodeTestProblem (gen, numRows, numCols, A.get(), A.lda(), true);

      if (debug)
        {
          cerr << "-- Generated test problem" << endl;
          // Don't print the matrix if it's too big.
          if (A.nrows() <= 30)
            {
              cerr << "A = " << endl;
              print_local_matrix (cerr, A.nrows(), A.ncols(),
                                  A.get(), A.lda());
              cerr << endl << endl;
            }
        }

      // Copy A into A_copy, since TSQR overwrites the input.  If
      // specified, rearrange the data in A_copy so that the data in
      // each cache block is contiguously stored.
      if (! contiguousCacheBlocks) {
        deep_copy (A_copy, A);
        if (debug) {
          cerr << "-- Copied test problem from A into A_copy" << endl;
          // Don't print the matrix if it's too big.
          if (A_copy.nrows() <= 30) {
            cerr << "A_copy = " << endl;
            print_local_matrix (cerr, A_copy.nrows(), A_copy.ncols(),
                                A_copy.get(), A_copy.lda());
            cerr << endl << endl;
          }
        }
      }
      else {
        actor.cache_block (numRows, numCols, A_copy.get(), A.get(), A.lda());
        if (debug) {
//.........这里部分代码省略.........

  void run(Teuchos::ParameterList &myMachPL, const Teuchos::RCP<const Teuchos::Comm<int> > &comm, const Teuchos::RCP<Node> &node) 
  {
    using std::plus;
    using std::endl;
    using Teuchos::null;
    using Teuchos::RCP;
    using Teuchos::ParameterList;
    using Teuchos::parameterList;
    using Tpetra::RTI::ZeroOp;

    // Static types
    typedef int                     LO;
    typedef int                     GO;
    typedef Tpetra::Map<LO,GO,Node>               Map;
    typedef Tpetra::CrsMatrix<S,LO,GO,Node> CrsMatrix;
    typedef Tpetra::Vector<S,LO,GO,Node>       Vector;
    typedef Teuchos::ScalarTraits<S> ST;

    IRTRdetails::fpu_fix<S> ff; ff.fix();

    *out << "Running test with Node==" << Teuchos::typeName(*node) << " on rank " << comm->getRank() << "/" << comm->getSize() << std::endl;

    // read the matrix
    RCP<CrsMatrix> A;
    RCP<const Map> rowMap = null;
    RCP<ParameterList> fillParams = parameterList();
    // must preserve the local graph in order to do convert() calls later
    if (Teuchos::TypeTraits::is_same<S,Sinner>::value) {
      fillParams->set("Preserve Local Graph",false);
    }
    else {
      fillParams->set("Preserve Local Graph",true);
    }
    Tpetra::Utils::readHBMatrix(matrixFile,comm,node,A,rowMap,fillParams);
    rowMap = A->getRowMap();

    testPassed = true;

    // compute an inital vector
    auto x = Tpetra::createVector<S>(rowMap);
    x->randomize();

    // call the solve
    S lambda = TpetraExamples::IRTR<Sinner,S,LO,GO,Node>(out,*params,A,x);

    // check that residual is as requested
    {
      auto r = Tpetra::createVector<S>(rowMap);
      A->apply(*x,*r);
      // compute A*x - x*lambda, while simultaneously computing |A*x - x*lambda|
      const S r_r = XFORM_REDUCE(r, x,                          // fused: 
                                 r - x*lambda,                  //      : r = r - x*lambda = A*x - x*lambda
                                 r*r, ZeroOp<S>, plus<S>() );   //      : sum r'*r
      const S rnrm = Teuchos::ScalarTraits<S>::squareroot(r_r);
      // check that residual is as requested
      *out << "|A*x - x*lambda|/|lambda|: " << rnrm / ST::magnitude(lambda) << endl;
      const double tolerance = params->get<double>("tolerance");
      if (rnrm / lambda > tolerance) testPassed = false;
    }

    ff.unfix();

    // print timings
    Teuchos::TimeMonitor::summarize( *out );
  }

void
cloneAndSolveWithBelos (
  bool& converged,
  int& numItersPerformed,
  const Teuchos::ScalarTraits<ST>::magnitudeType& tol,
  const int maxNumIters,
  const int num_steps,
  const Teuchos::RCP<CloneNode>& clone_node,
  const Teuchos::RCP<multivector_type>& X,
  const Teuchos::RCP<const sparse_matrix_type>& A,
  const Teuchos::RCP<const multivector_type>& B,
  const std::string& prec_type,
  const Teuchos::RCP<const operator_type>& M_left=Teuchos::null,
  const Teuchos::RCP<const operator_type>& M_right=Teuchos::null) {

  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;
  using Teuchos::rcp_dynamic_cast;

  typedef Tpetra::CrsMatrix<ST, LO, GO, CloneNode>    clone_sparse_matrix_type;
  typedef Tpetra::Operator<ST, LO, GO, CloneNode>     clone_operator_type;
  typedef Tpetra::MultiVector<ST, LO, GO, CloneNode>  clone_multi_vector_type;
#ifdef HAVE_TRILINOSCOUPLINGS_MUELU
  typedef typename KokkosClassic::DefaultKernels<ST,LO,CloneNode>::SparseOps clone_sparse_ops;
#endif // HAVE_TRILINOSCOUPLINGS_MUELU
  typedef clone_multi_vector_type MV;
  typedef clone_operator_type OP;

  // Clone Matrix, RHS, LHS
  RCP<ParameterList> plClone = parameterList();
  RCP<clone_sparse_matrix_type> A_clone;
  RCP<clone_multi_vector_type> B_clone, X_clone;
  {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Clone System", clone_system);
    A_clone = A->clone(clone_node, plClone);
    B_clone = B->clone(clone_node);
    X_clone = X->clone(clone_node);
  }

  // Clone preconditioner(s)
  RCP<const clone_operator_type> M_left_clone, M_right_clone;
  {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Clone Preconditioner", clone_prec);

#ifdef HAVE_TRILINOSCOUPLINGS_MUELU
    if (M_left != Teuchos::null && prec_type == "MueLu") {
      RCP< const MueLu::TpetraOperator<ST,LO,GO,Node> > M_muelu =
        rcp_dynamic_cast<const MueLu::TpetraOperator<ST,LO,GO,Node> >(M_left);
      M_left_clone = M_muelu->clone<CloneNode, clone_sparse_ops>(clone_node);
    }
    if (M_right != Teuchos::null && prec_type == "MueLu") {
      RCP< const MueLu::TpetraOperator<ST,LO,GO,Node> > M_muelu =
        rcp_dynamic_cast<const MueLu::TpetraOperator<ST,LO,GO,Node> >(M_right);
      M_right_clone = M_muelu->clone<CloneNode, clone_sparse_ops>(clone_node);
    }
#else
    TEUCHOS_TEST_FOR_EXCEPTION(
      prec_type == "MueLu", std::runtime_error, "Tpetra scaling example: "
      "In order to precondition with MueLu, you must have built Trilinos "
      "with the MueLu package enabled.");
#endif // HAVE_TRILINOSCOUPLINGS_MUELU
  }

  // Solve
  {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Clone Solve", clone_solve);
    IntrepidPoissonExample::solveWithBelos<ST,MV,OP> (
      converged, numItersPerformed, tol, maxNumIters, num_steps,
      X_clone, A_clone, B_clone, M_left_clone, M_right_clone);
  }

  // Copy X_clone back into X
  {
    TEUCHOS_FUNC_TIME_MONITOR_DIFF("Clone Solution", clone_sol);
    RCP<multivector_type> X_host = X_clone->clone(X->getMap()->getNode());
    X->update(1.0, *X_host, 0.0);
  }
}

Piro::RythmosSolver<Scalar>::RythmosSolver(Teuchos::RCP<Teuchos::ParameterList> in_appParams,
                          Teuchos::RCP< Thyra::ModelEvaluatorDefaultBase<Scalar> > in_model,
                          Teuchos::RCP<Rythmos::IntegrationObserverBase<Scalar> > in_observer) :
  appParams(in_appParams),
  model(in_model),
  observer(in_observer)
{
  // For dumping default parameters from Rythmos
  {
    //Rythmos::IntegratorBuilder<double> b;
    //std::cout << *(b.getValidParameters()) << std::endl;
    //Teuchos::writeParameterListToXmlFile(*b.getValidParameters(), "sample.xml");
  }

  using Teuchos::ParameterList;
  using Teuchos::parameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;

  out = Teuchos::VerboseObjectBase::getDefaultOStream();

  num_p = model->createInArgs().Np();
  num_g = model->createOutArgs().Ng();

//   TEUCHOS_TEST_FOR_EXCEPTION(num_p > 1, Teuchos::Exceptions::InvalidParameter,
//                      std::endl << "Error in Piro::RythmosSolver " <<
//                      "Not Implemented for Np>1 : " << num_p << std::endl);
//   TEUCHOS_TEST_FOR_EXCEPTION(num_g > 1, Teuchos::Exceptions::InvalidParameter,
//                      std::endl << "Error in Piro::RythmosSolver " <<
//                      "Not Implemented for Ng>1 : " << num_g << std::endl);

  //
  *out << "\nA) Get the base parameter list ...\n";
  //

  RCP<Teuchos::ParameterList> rythmosPL = sublist(appParams, "Rythmos", true);
  rythmosPL->validateParameters(*getValidRythmosParameters(),0);

  {
    const std::string verbosity = rythmosPL->get("Verbosity Level", "VERB_DEFAULT");
    solnVerbLevel = Teuchos::VERB_DEFAULT;
    if      (verbosity == "VERB_NONE")    solnVerbLevel = Teuchos::VERB_NONE;
    else if (verbosity == "VERB_LOW")     solnVerbLevel = Teuchos::VERB_LOW;
    else if (verbosity == "VERB_MEDIUM")  solnVerbLevel = Teuchos::VERB_MEDIUM;
    else if (verbosity == "VERB_HIGH")    solnVerbLevel = Teuchos::VERB_HIGH;
    else if (verbosity == "VERB_EXTREME") solnVerbLevel = Teuchos::VERB_EXTREME;
  }

  t_final = rythmosPL->get("Final Time", 0.1);
  
  const std::string stepperType = rythmosPL->get("Stepper Type", "Backward Euler");
  
  //
  *out << "\nC) Create and initalize the forward model ...\n";
  //
      
  *out << "\nD) Create the stepper and integrator for the forward problem ...\n";
  //
  
  if (rythmosPL->get<std::string>("Nonlinear Solver Type") == "Rythmos") {
    Teuchos::RCP<Rythmos::TimeStepNonlinearSolver<double> > rythmosTimeStepSolver = 
      Rythmos::timeStepNonlinearSolver<double>();
    if (rythmosPL->getEntryPtr("NonLinear Solver")) {
      RCP<Teuchos::ParameterList> nonlinePL =
	sublist(rythmosPL, "NonLinear Solver", true);
      rythmosTimeStepSolver->setParameterList(nonlinePL);
    }
    fwdTimeStepSolver = rythmosTimeStepSolver;
  }
  else if (rythmosPL->get<std::string>("Nonlinear Solver Type") == "NOX") {
#ifdef Piro_ENABLE_NOX
    Teuchos::RCP<Thyra::NOXNonlinearSolver> nox_solver =  Teuchos::rcp(new Thyra::NOXNonlinearSolver);
    Teuchos::RCP<Teuchos::ParameterList> nox_params = Teuchos::rcp(new Teuchos::ParameterList);
    *nox_params = appParams->sublist("NOX");
    nox_solver->setParameterList(nox_params);
    fwdTimeStepSolver = nox_solver;
#else
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,"Requested NOX solver for a Rythmos Transient solve, Trilinos was not built with NOX enabled.  Please rebuild Trilinos or use the native Rythmos nonlinear solver.");
#endif
    
  }
  
  if (stepperType == "Backward Euler") {
    fwdStateStepper = Rythmos::backwardEulerStepper<Scalar> (model, fwdTimeStepSolver);
    fwdStateStepper->setParameterList(sublist(rythmosPL, "Rythmos Stepper", true));
  }
  else if (stepperType == "Explicit RK") {
    fwdStateStepper = Rythmos::explicitRKStepper<Scalar>(model);
    fwdStateStepper->setParameterList(sublist(rythmosPL, "Rythmos Stepper", true));
  }
  else if (stepperType == "BDF") {
    Teuchos::RCP<Teuchos::ParameterList> BDFparams = 
      Teuchos::sublist(rythmosPL, "Rythmos Stepper", true);
    Teuchos::RCP<Teuchos::ParameterList> BDFStepControlPL =
      Teuchos::sublist(BDFparams,"Step Control Settings");
    
    fwdStateStepper = Teuchos::rcp( new Rythmos::ImplicitBDFStepper<Scalar>(model,fwdTimeStepSolver,BDFparams) );
    fwdStateStepper->setInitialCondition(model->getNominalValues());
    
  }
//.........这里部分代码省略.........