C++ Teuchos类代码示例

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

  using std::endl;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;
  using Teuchos::CommandLineProcessor;
  
  bool result, success = true;

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

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

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

  try {

    //
    // Read commandline options
    //

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

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

    bool dumpFinalSolutions = false;
    clp.setOption(
      "dump-final-solutions", "no-dump-final-solutions", &dumpFinalSolutions,
      "Determine if the final solutions are dumpped or not." );
    
    CommandLineProcessor::EParseCommandLineReturn parse_return = clp.parse(argc,argv);
    if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL ) return parse_return;
    
    if ( Teuchos::VERB_DEFAULT == verbLevel )
      verbLevel = Teuchos::VERB_LOW;

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

    //
    // Get the base parameter list that all other parameter lists will be read
    // from.
    //
    
    RCP<ParameterList>
      paramList = Teuchos::parameterList();

    //
    // Create the underlying EpetraExt::ModelEvaluator
    //

    RCP<LorenzModel>
      lorenzModel = rcp(new LorenzModel( epetra_comm, *paramList ));

    //
    // Create the Thyra-wrapped ModelEvaluator
    //
    
    RCP<Thyra::ModelEvaluator<double> >
      thyraLorenzModel = Thyra::epetraModelEvaluator(lorenzModel,Teuchos::null);

    //
    // Create the Rythmos GAASP ErrorEstimator
    //

    RCP<Rythmos::GAASPErrorEstimator> gaaspEE = rcp(new Rythmos::GAASPErrorEstimator);
    gaaspEE->setModel(thyraLorenzModel);
    //gaaspEE->setQuantityOfInterest( AVERAGE_ERROR_QTY );  // Not passed through yet.
    Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
    pl->sublist("GAASP Interface Parameters").set<double>("eTime",1.0);
    pl->sublist("GAASP Interface Parameters").set<double>("timeStep",0.1);
    //pl->sublist("GAASP Interface Parameters").set<double>("timeStep",0.5);
    gaaspEE->setParameterList(pl);

    //RCP<const Rythmos::ErrorEstimateBase<double> > error = gaaspEE->getErrorEstimate();

    //double uTOL = 1.0e-8;
    double uTOL = 1.0e-2;
    RCP<const Rythmos::ErrorEstimateBase<double> > error = gaaspEE->controlGlobalError(uTOL);


    double err = error->getTotalError();
    out->precision(15);
    *out << "err = " << err << std::endl;

  }
  TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,success);
//.........这里部分代码省略.........

TEUCHOS_UNIT_TEST(Belos_Hypre, Laplace2D){
  const double tol = 1E-7;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;
  typedef Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator>  LinearProblem;

  //
  // Create Laplace2D
  //
#ifdef HAVE_MPI
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm();
#endif
  Teuchos::ParameterList GaleriList;
  int nx = 10 * Comm.NumProc();
  int ny = 10 * Comm.NumProc();
  GaleriList.set("nx", nx);
  GaleriList.set("ny", ny);
  Epetra_Map Map(nx*ny,0,Comm);
  RCP<Epetra_CrsMatrix> Crs_Matrix   = rcp(Galeri::CreateCrsMatrix("Laplace2D", &Map, GaleriList));
  int NumProc = Crs_Matrix->Comm().NumProc();

  //
  // Create the hypre preconditioner
  //
  RCP<Ifpack_Hypre> preconditioner = rcp(new Ifpack_Hypre(Crs_Matrix.get()));
  TEST_EQUALITY(preconditioner->Initialize(),0);
  TEST_EQUALITY(preconditioner->SetParameter(Preconditioner, ParaSails),0); // Use a Euclid Preconditioner (but not really used)
  TEST_EQUALITY(preconditioner->SetParameter(Preconditioner),0); // Solve the problem
  TEST_EQUALITY(preconditioner->Compute(),0);

  //
  // Create the solution vector and rhs
  //
  int numVec = 1;
  RCP<Epetra_MultiVector> X = rcp(new Epetra_MultiVector(Crs_Matrix->OperatorDomainMap(), numVec));
  RCP<Epetra_MultiVector> KnownX = rcp(new Epetra_MultiVector(Crs_Matrix->OperatorDomainMap(), numVec));
  KnownX->Random();
  RCP<Epetra_MultiVector> B = rcp(new Epetra_MultiVector(Crs_Matrix->OperatorRangeMap(), numVec));
  Crs_Matrix->Apply(*KnownX, *B);

  //
  // Test the EpetraExt wrapper
  // amk November 24, 2015: Should we deprecate this?
  //
//  RCP<ParameterList> pl = rcp(new ParameterList());
//  TEST_EQUALITY(X->PutScalar(0.0),0);
//  HYPRE_IJMatrix hypre_mat = preconditioner->HypreMatrix();
//  RCP<EpetraExt_HypreIJMatrix> Hyp_Matrix = rcp(new EpetraExt_HypreIJMatrix(hypre_mat));
//  TEST_EQUALITY(Hyp_Matrix->SetParameter(Preconditioner, ParaSails),0);
//  TEST_EQUALITY(Hyp_Matrix->SetParameter(Preconditioner),0);
//  TEST_EQUALITY(EquivalentMatrices(*Hyp_Matrix, *Crs_Matrix, tol), true);
//  RCP<LinearProblem> problem1 = rcp(new LinearProblem(Crs_Matrix,X,B));
//  problem1->setLeftPrec(Hyp_Matrix);
//  TEST_EQUALITY(problem1->setProblem(),true);
//  Belos::PseudoBlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator> solMgr1(problem1,pl);
//  Belos::ReturnType rv1 = solMgr1.solve(); // TEST_EQUALITY(solMgr2.solve(),Belos::Converged);
//  TEST_EQUALITY(rv1,Belos::Converged);
//  TEST_EQUALITY(EquivalentVectors(*X, *KnownX, tol*10*pow(10.0,NumProc)), true);

  //
  // Test the Ifpack hypre interface
  //
  RCP<ParameterList> pl2 = rcp(new ParameterList());
  RCP<Epetra_Operator> invOp = rcp(new Epetra_InvOperator(preconditioner.get()));
  TEST_EQUALITY(X->PutScalar(0.0),0);
  RCP<LinearProblem> problem2 = rcp(new LinearProblem(Crs_Matrix,X,B));
  problem2->setLeftPrec(invOp);
  TEST_EQUALITY(problem2->setProblem(),true);
  Belos::PseudoBlockGmresSolMgr<double,Epetra_MultiVector,Epetra_Operator> solMgr2(problem2,pl2);
  Belos::ReturnType rv2 = solMgr2.solve(); // TEST_EQUALITY(solMgr2.solve(),Belos::Converged);
  TEST_EQUALITY(rv2,Belos::Converged);
  TEST_EQUALITY(EquivalentVectors(*X, *KnownX, tol*10*pow(10.0,NumProc)), true);
}

TEUCHOS_UNIT_TEST(tIterativePreconditionerFactory, parameter_list_init)
{
   // build global (or serial communicator)
   #ifdef HAVE_MPI
      Epetra_MpiComm Comm(MPI_COMM_WORLD);
   #else
      Epetra_SerialComm Comm;
   #endif

   Teko::LinearOp  A = build2x2(Comm,1,2,3,4);

   Thyra::LinearOpTester<double> tester;
   tester.show_all_tests(true);

   {
      RCP<Teuchos::ParameterList> pl = buildLibPL(4,"Amesos");
      RCP<Teko::IterativePreconditionerFactory> precFact = rcp(new Teko::IterativePreconditionerFactory());
      RCP<Teko::InverseFactory> invFact = rcp(new Teko::PreconditionerInverseFactory(precFact,Teuchos::null));

      try {
         precFact->initializeFromParameterList(*pl);
         out << "Passed correct parameter list" << std::endl;

         Teko::LinearOp prec = Teko::buildInverse(*invFact,A);
      }
      catch(...) {
         success = false; 
         out << "Failed correct parameter list" << std::endl;
      }
   }

   {
      Teuchos::ParameterList pl;
      pl.set("Preconditioner Type","Amesos");

      RCP<Teko::IterativePreconditionerFactory> precFact = rcp(new Teko::IterativePreconditionerFactory());
      RCP<Teko::InverseFactory> invFact = rcp(new Teko::PreconditionerInverseFactory(precFact,Teuchos::null));

      try {
         precFact->initializeFromParameterList(pl);
         out << "Passed iteration count" << std::endl;

         Teko::LinearOp prec = Teko::buildInverse(*invFact,A);
      }
      catch(...) {
         out << "Failed iteration count" << std::endl;
      }
   }

   {
      Teuchos::ParameterList pl;
      pl.set("Iteration Count",4);
      pl.set("Precondiioner Type","Amesos");

      RCP<Teko::IterativePreconditionerFactory> precFact = rcp(new Teko::IterativePreconditionerFactory());

      try {
         precFact->initializeFromParameterList(pl);
         success = false;
         out << "Failed preconditioner type" << std::endl;

         // these should not be executed
         RCP<Teko::InverseFactory> invFact = rcp(new Teko::PreconditionerInverseFactory(precFact,Teuchos::null));
         Teko::LinearOp prec = Teko::buildInverse(*invFact,A);
      }
      catch(const std::exception & exp) {
         out << "Passed preconditioner type" << std::endl;
      }
   }
}

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

  typedef double Scalar;
  typedef Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
  typedef Tpetra::Map<>::local_ordinal_type LO;
  typedef Tpetra::Map<>::global_ordinal_type GO;
  typedef Tpetra::DefaultPlatform::DefaultPlatformType Platform;

  typedef Tpetra::CrsMatrix<Scalar,LO,GO> MAT;
  typedef Tpetra::MultiVector<Scalar,LO,GO> MV;

  using Tpetra::global_size_t;
  using Teuchos::tuple;
  using Teuchos::RCP;
  using Teuchos::rcp;

  Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
  Teuchos::RCP<const Teuchos::Comm<int> > comm = platform.getComm();
  size_t myRank = comm->getRank();

  RCP<Teuchos::FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));

  *fos << Amesos2::version() << std::endl << std::endl;

  bool printMatrix   = false;
  bool printSolution = false;
  bool printTiming   = false;
  bool printResidual = false;
  bool printLUStats  = false;
  bool verbose       = false;
  std::string solver_name = "SuperLU";
  std::string filedir = "../test/matrices/";
  std::string filename = "arc130.mtx";
  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("filedir",&filedir,"Directory where matrix-market files are located");
  cmdp.setOption("filename",&filename,"Filename for Matrix-Market test matrix.");
  cmdp.setOption("print-matrix","no-print-matrix",&printMatrix,"Print the full matrix after reading it.");
  cmdp.setOption("print-solution","no-print-solution",&printSolution,"Print solution vector after solve.");
  cmdp.setOption("print-timing","no-print-timing",&printTiming,"Print solver timing statistics");
  cmdp.setOption("print-residual","no-print-residual",&printResidual,"Print solution residual");
  cmdp.setOption("print-lu-stats","no-print-lu-stats",&printLUStats,"Print nnz in L and U factors");
  cmdp.setOption("solver", &solver_name, "Which TPL solver library to use.");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }

  // Before we do anything, check that the solver is enabled
  if( !Amesos2::query(solver_name) ){
    std::cerr << solver_name << " not enabled.  Exiting..." << std::endl;
    return EXIT_SUCCESS;        // Otherwise CTest will pick it up as
                                // failure, which it isn't really
  }

  const size_t numVectors = 1;

  // create a Map
  global_size_t nrows = 6;
  RCP<Tpetra::Map<LO,GO> > map
    = rcp( new Tpetra::Map<LO,GO>(nrows,0,comm) );

  std::string mat_pathname = filedir + filename;
  RCP<MAT> A = Tpetra::MatrixMarket::Reader<MAT>::readSparseFile(mat_pathname,comm);

  if( printMatrix ){
    A->describe(*fos, Teuchos::VERB_EXTREME);
  }
  else if( verbose && myRank==0 ){
    *fos << std::endl << A->description() << std::endl << std::endl;
  }

  // get the maps
  RCP<const Tpetra::Map<LO,GO> > dmnmap = A->getDomainMap();
  RCP<const Tpetra::Map<LO,GO> > rngmap = A->getRangeMap();

  // Create random X
  RCP<MV> Xhat = rcp( new MV(dmnmap,numVectors) );
  RCP<MV> X = rcp( new MV(dmnmap,numVectors) );
  X->setObjectLabel("X");
  Xhat->setObjectLabel("Xhat");
  X->randomize();

  RCP<MV> B = rcp(new MV(rngmap,numVectors));
  A->apply(*X, *B);

  // Constructor from Factory
  RCP<Amesos2::Solver<MAT,MV> > solver;
  try{
    solver = Amesos2::create<MAT,MV>(solver_name, A, Xhat, B);
  } catch (std::invalid_argument e){
    *fos << e.what() << std::endl;
    return 0;
  }

  #ifdef SHYLU_NODEBASKER
  if( Amesos2::query("Basker") ) {
    Teuchos::ParameterList amesos2_params("Amesos2");
      amesos2_params.sublist(solver_name).set("num_threads", 1, "Number of threads");
    solver->setParameters( Teuchos::rcpFromRef(amesos2_params) );
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
  //
  int MyPID = 0;
#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
  MyPID = Comm.MyPID();
#else
  Epetra_SerialComm Comm;
#endif
  //
  typedef double                            ST;
  typedef Teuchos::ScalarTraits<ST>        SCT;
  typedef SCT::magnitudeType                MT;
  typedef Epetra_MultiVector                MV;
  typedef Epetra_Operator                   OP;
  typedef Belos::MultiVecTraits<ST,MV>     MVT;
  typedef Belos::OperatorTraits<ST,MV,OP>  OPT;

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

  bool success = false;
  bool verbose = false;
  try {
    bool proc_verbose = false;
    int frequency = -1;        // frequency of status test output.
    int blocksize = 1;         // blocksize
    int numrhs = 1;            // number of right-hand sides to solve for
    int maxrestarts = 15;      // maximum number of restarts allowed
    int maxiters = -1;         // maximum number of iterations allowed per linear system
    int maxsubspace = 25;      // maximum number of blocks the solver can use for the subspace
    std::string filename("orsirr1.hb");
    MT tol = 1.0e-5;           // relative residual tolerance

    // Specify whether to use RHS as initial guess. If false, use zero.
    bool useRHSAsInitialGuess = false;

    Teuchos::CommandLineProcessor cmdp(false,true);
    cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
    cmdp.setOption("use-rhs","use-zero",&useRHSAsInitialGuess,"Use RHS as initial guess.");
    cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
    cmdp.setOption("filename",&filename,"Filename for test matrix.  Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
    cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
    cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
    cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
    cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
    cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
    cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
    if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
      return EXIT_FAILURE;
    }

    if (!verbose)
      frequency = -1;  // reset frequency if test is not verbose

    //
    // *************Get the problem*********************
    //
    RCP<Epetra_Map> Map;
    RCP<Epetra_CrsMatrix> A;
    RCP<Epetra_MultiVector> B, X;
    RCP<Epetra_Vector> vecB, vecX;
    EpetraExt::readEpetraLinearSystem(filename, Comm, &A, &Map, &vecX, &vecB);
    A->OptimizeStorage();
    proc_verbose = verbose && (MyPID==0);  /* Only print on the zero processor */

    // Check to see if the number of right-hand sides is the same as requested.
    if (numrhs>1) {
      X = rcp( new Epetra_MultiVector( *Map, numrhs ) );
      B = rcp( new Epetra_MultiVector( *Map, numrhs ) );
      X->Random();
      OPT::Apply( *A, *X, *B );
      X->PutScalar( 0.0 );
    }
    else {
      X = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecX);
      B = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecB);
    }

    // If requested, use a copy of B as initial guess
    if (useRHSAsInitialGuess)
    {
      X->Update(1.0, *B, 0.0);
    }

    //
    // ************Construct preconditioner*************
    //
    ParameterList ifpackList;

    // allocates an IFPACK factory. No data is associated
    // to this object (only method Create()).
    Ifpack Factory;

    // create the preconditioner. For valid PrecType values,
    // please check the documentation
    std::string PrecType = "ILU"; // incomplete LU
//.........这里部分代码省略.........

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

  typedef Tpetra::Vector<SC,LO,GO,NO>                  TVEC;
  typedef Tpetra::MultiVector<SC,LO,GO,NO>             TMV;
  typedef Tpetra::CrsMatrix<SC,LO,GO,NO,LMO>           TCRS;
  typedef Xpetra::CrsMatrix<SC,LO,GO,NO,LMO>           XCRS;
  typedef Xpetra::TpetraCrsMatrix<SC,LO,GO,NO,LMO>     XTCRS;
  typedef Xpetra::Matrix<SC,LO,GO,NO,LMO>              XMAT;
  typedef Xpetra::CrsMatrixWrap<SC,LO,GO,NO,LMO>       XWRAP;

  typedef Belos::OperatorT<TMV>                        TOP;
  typedef Belos::OperatorTraits<SC,TMV,TOP>            TOPT;
  typedef Belos::MultiVecTraits<SC,TMV>                TMVT;
  typedef Belos::LinearProblem<SC,TMV,TOP>             TProblem;
  typedef Belos::SolverManager<SC,TMV,TOP>             TBelosSolver;
  typedef Belos::BlockGmresSolMgr<SC,TMV,TOP>          TBelosGMRES;

  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::TimeMonitor;

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

  GO nx,ny,nz;
  nx=100; ny=100; nz=100;
  double stretchx, stretchy, stretchz, h, delta;
  stretchx=1.0; stretchy=1.0; stretchz=1.0;
  h=0.01; delta=2.0;
  int PMLXL, PMLXR, PMLYL, PMLYR, PMLZL, PMLZR;
  PMLXL=10; PMLXR=10; PMLYL=10; PMLYR=10; PMLZL=10; PMLZR=10;
  double omega, shift;
  omega=20.0*M_PI;
  shift=0.5;

  Galeri::Xpetra::Parameters<GO> matrixParameters(clp, nx, ny, nz, "Helmholtz1D", 0, stretchx, stretchy, stretchz,
						  h, delta, PMLXL, PMLXR, PMLYL, PMLYR, PMLZL, PMLZR, omega, shift);
  Xpetra::Parameters             xpetraParameters(clp);

  RCP<TimeMonitor> globalTimeMonitor = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time")));
  RCP<TimeMonitor> tm = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));

  Teuchos::ParameterList pl = matrixParameters.GetParameterList();
  RCP<MultiVector> coordinates;
  Teuchos::ParameterList galeriList;
  galeriList.set("nx", pl.get("nx", nx));
  galeriList.set("ny", pl.get("ny", ny));
  galeriList.set("nz", pl.get("nz", nz));
  RCP<const Map> map;

  if (matrixParameters.GetMatrixType() == "Helmholtz1D") {
    map = MapFactory::Build(xpetraParameters.GetLib(), matrixParameters.GetNumGlobalElements(), 0, comm);
    coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC, LO, GO, Map, MultiVector>("1D", map, matrixParameters.GetParameterList());
  }
  else if (matrixParameters.GetMatrixType() == "Helmholtz2D") {
    map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(), "Cartesian2D", comm, galeriList);
    coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC, LO, GO, Map, MultiVector>("2D", map, matrixParameters.GetParameterList());
  }
  else if (matrixParameters.GetMatrixType() == "Helmholtz3D") {
    map = Galeri::Xpetra::CreateMap<LO, GO, Node>(xpetraParameters.GetLib(), "Cartesian3D", comm, galeriList);
    coordinates = Galeri::Xpetra::Utils::CreateCartesianCoordinates<SC, LO, GO, Map, MultiVector>("3D", map, matrixParameters.GetParameterList());
  }

  RCP<const Tpetra::Map<LO, GO, NO> > tmap = Xpetra::toTpetra(map);

  Teuchos::ParameterList matrixParams = matrixParameters.GetParameterList();

  // Build problem
  RCP<Galeri::Xpetra::Problem_Helmholtz<Map,CrsMatrixWrap,MultiVector> > Pr =
      Galeri::Xpetra::BuildProblem_Helmholtz<SC,LO,GO,Map,CrsMatrixWrap,MultiVector>(matrixParameters.GetMatrixType(), map, matrixParams);
  RCP<Matrix> A = Pr->BuildMatrix();

  RCP<MultiVector> nullspace = MultiVectorFactory::Build(map,1);
  nullspace->putScalar( (SC) 1.0);
 
  comm->barrier();

  tm = Teuchos::null;

  // Construct a multigrid preconditioner
  tm = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 2 - MueLu Setup")));

  // Multigrid Hierarchy
  RCP<Hierarchy> H = rcp(new Hierarchy(A));
  H->GetLevel(0)->Set("Nullspace",nullspace);
  FactoryManager Manager;
  H->Setup(Manager, 0, 5);
  //H->Write(-1,-1);

  tm = Teuchos::null;

  // Solve Ax = b
  tm = rcp (new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 3 - LHS and RHS initialization")));
  RCP<TVEC> X = Tpetra::createVector<SC,LO,GO,NO>(tmap);
  RCP<TVEC> B = Tpetra::createVector<SC,LO,GO,NO>(tmap);  
  X->putScalar((SC) 0.0);
  B->putScalar((SC) 0.0);
  if(comm->getRank()==0) {
//.........这里部分代码省略.........

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

  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ArrayRCP;
  using Teuchos::TimeMonitor;
  using Teuchos::ParameterList;

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

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

  // =========================================================================
  // Parameters initialization
  // =========================================================================
  Teuchos::CommandLineProcessor clp(false);

  GO nx = 100, ny = 100, nz = 100;
  Galeri::Xpetra::Parameters<GO> galeriParameters(clp, nx, ny, nz, "Laplace2D"); // manage parameters of the test case
  Xpetra::Parameters             xpetraParameters(clp);                          // manage parameters of Xpetra

  std::string xmlFileName       = "scalingTest.xml"; clp.setOption("xml",                   &xmlFileName,      "read parameters from a file [default = 'scalingTest.xml']");
  bool        printTimings      = true;              clp.setOption("timings", "notimings",  &printTimings,     "print timings to screen");
  int         writeMatricesOPT  = -2;                clp.setOption("write",                 &writeMatricesOPT, "write matrices to file (-1 means all; i>=0 means level i)");
  std::string dsolveType        = "cg", solveType;   clp.setOption("solver",                &dsolveType,       "solve type: (none | cg | gmres | standalone)");
  double      dtol              = 1e-12, tol;        clp.setOption("tol",                   &dtol,             "solver convergence tolerance");

  std::string mapFile;                               clp.setOption("map",                   &mapFile,          "map data file");
  std::string matrixFile;                            clp.setOption("matrix",                &matrixFile,       "matrix data file");
  std::string coordFile;                             clp.setOption("coords",                &coordFile,        "coordinates data file");
  int         numRebuilds       = 0;                 clp.setOption("rebuild",               &numRebuilds,      "#times to rebuild hierarchy");
  int         maxIts            = 200;               clp.setOption("its",                   &maxIts,           "maximum number of solver iterations");
  bool        scaleResidualHistory = true;              clp.setOption("scale", "noscale",  &scaleResidualHistory, "scaled Krylov residual history");

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

  Xpetra::UnderlyingLib lib = xpetraParameters.GetLib();

  ParameterList paramList;
  Teuchos::updateParametersFromXmlFileAndBroadcast(xmlFileName, Teuchos::Ptr<ParameterList>(&paramList), *comm);
  bool isDriver = paramList.isSublist("Run1");
  if (isDriver) {
    // update galeriParameters with the values from the XML file
    ParameterList& realParams = galeriParameters.GetParameterList();

    for (ParameterList::ConstIterator it = realParams.begin(); it != realParams.end(); it++) {
      const std::string& name = realParams.name(it);
      if (paramList.isParameter(name))
        realParams.setEntry(name, paramList.getEntry(name));
    }
  }

  // Retrieve matrix parameters (they may have been changed on the command line)
  // [for instance, if we changed matrix type from 2D to 3D we need to update nz]
  ParameterList galeriList = galeriParameters.GetParameterList();

  // =========================================================================
  // Problem construction
  // =========================================================================
  std::ostringstream galeriStream;
  comm->barrier();
  RCP<TimeMonitor> globalTimeMonitor = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: S - Global Time")));
  RCP<TimeMonitor> tm                = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1 - Matrix Build")));

  RCP<Matrix>      A;
  RCP<const Map>   map;
  RCP<MultiVector> coordinates;
  RCP<MultiVector> nullspace;
  if (matrixFile.empty()) {
    galeriStream << "========================================================\n" << xpetraParameters << galeriParameters;

    // Galeri will attempt to create a square-as-possible distribution of subdomains di, e.g.,
    //                                 d1  d2  d3
    //                                 d4  d5  d6
    //                                 d7  d8  d9
    //                                 d10 d11 d12
    // A perfect distribution is only possible when the #processors is a perfect square.
    // This *will* result in "strip" distribution if the #processors is a prime number or if the factors are very different in
    // size. For example, np=14 will give a 7-by-2 distribution.
    // If you don't want Galeri to do this, specify mx or my on the galeriList.
    std::string matrixType = galeriParameters.GetMatrixType();

    // Create map and coordinates
    // In the future, we hope to be able to first create a Galeri problem, and then request map and coordinates from it
    // At the moment, however, things are fragile as we hope that the Problem uses same map and coordinates inside
    if (matrixType == "Laplace1D") {
//.........这里部分代码省略.........

void
SupportGraph<MatrixType>::findSupport ()
{
  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
                            global_ordinal_type, node_type> crs_matrix_type;
  typedef Tpetra::Vector<scalar_type, local_ordinal_type,
                         global_ordinal_type, node_type> vec_type;
  
typedef std::pair<int, int> E;
  typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS, 
                                boost::no_property, 
                                boost::property<boost::edge_weight_t, 
                                                magnitude_type> > graph_type;
  typedef typename boost::graph_traits<graph_type>::edge_descriptor edge_type;
  typedef typename boost::graph_traits<graph_type>::vertex_descriptor 
    vertex_type;

  const scalar_type zero = STS::zero();
  const scalar_type one = STS::one();

  //Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cout));

  size_t num_verts = A_local_->getNodeNumRows();
  size_t num_edges  
    = (A_local_->getNodeNumEntries() - A_local_->getNodeNumDiags())/2;

  // Create data structures for the BGL code
  // and temp data structures for extraction
  E *edge_array = new E[num_edges];
  magnitude_type *weights = new magnitude_type[num_edges];

  size_t num_entries;
  size_t max_num_entries = A_local_->getNodeMaxNumRowEntries();

  std::vector<scalar_type> valuestemp (max_num_entries);
  std::vector<local_ordinal_type> indicestemp (max_num_entries);
  
  std::vector<magnitude_type> diagonal (num_verts);

  Tpetra::ArrayView<scalar_type> values (valuestemp);
  Tpetra::ArrayView<local_ordinal_type> indices (indicestemp);

  // Extract from the tpetra matrix keeping only one edge per pair 
  // (assume symmetric)
  size_t offDiagCount = 0;
  for (size_t row = 0; row < num_verts; ++row) {
    A_local_->getLocalRowCopy (row, indices, values, num_entries);
    for (size_t colIndex = 0; colIndex < num_entries; ++colIndex) {
      if(row == Teuchos::as<size_t>(indices[colIndex])) {
        diagonal[row] = values[colIndex];
      }

      if((row < Teuchos::as<size_t>(indices[colIndex])) 
         && (values[colIndex] < zero)) {
        edge_array[offDiagCount] = E(row, indices[colIndex]);
        weights[offDiagCount] = values[colIndex];
        if (Randomize_) {
          // Add small random pertubation.
          weights[offDiagCount] *= one + 
            STS::magnitude(STS::rmin() * STS::random());
        }

        offDiagCount++;
      }
    }
  }

  // Create BGL graph
  graph_type g(edge_array, edge_array + num_edges, weights, num_verts);
  typedef typename boost::property_map 
    <graph_type, boost::edge_weight_t>::type type;
  type weight = get (boost::edge_weight, g);
  std::vector<edge_type> spanning_tree;

  // Run Kruskal, actually maximal weight ST since edges are negative
  boost::kruskal_minimum_spanning_tree(g, std::back_inserter (spanning_tree));

  // Create array to store the exact number of non-zeros per row
  Teuchos::ArrayRCP<size_t> NumNz (num_verts, 1);

  typedef typename std::vector<edge_type>::iterator edge_iterator_type;

  // Find the degree of all the vertices
  for (edge_iterator_type ei = spanning_tree.begin(); ei != spanning_tree.end();
       ++ei) {
    local_ordinal_type localsource = source(*ei,g);
    local_ordinal_type localtarget = target(*ei,g);

    // We only want upper triangular entries, might need to swap
    if (localsource > localtarget) {
      localsource = target(*ei, g);
      localtarget = source(*ei, g);
    }

    NumNz[localsource] += 1;
  }


  // Create an stl vector of stl vectors to hold indices and values
  std::vector<std::vector<local_ordinal_type> > Indices (num_verts); 
//.........这里部分代码省略.........

void
SupportGraph<MatrixType>::
apply (const Tpetra::MultiVector<scalar_type,
                                 local_ordinal_type,
                                 global_ordinal_type,
                                 node_type>& X,
       Tpetra::MultiVector<scalar_type,
                           local_ordinal_type,
                           global_ordinal_type,
                           node_type>& Y,
       Teuchos::ETransp mode,
       scalar_type alpha,
       scalar_type beta) const
{
  using Teuchos::FancyOStream;
  using Teuchos::getFancyOStream;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::rcpFromRef;
  using Teuchos::Time;
  using Teuchos::TimeMonitor;
  typedef scalar_type DomainScalar;
  typedef scalar_type RangeScalar;
  typedef Tpetra::MultiVector<DomainScalar, local_ordinal_type,
    global_ordinal_type, node_type> MV;

  RCP<FancyOStream> out = getFancyOStream(rcpFromRef(std::cout));

  // Create a timer for this method, if it doesn't exist already.
  // TimeMonitor::getNewCounter registers the timer, so that
  // TimeMonitor's class methods like summarize() will report the
  // total time spent in successful calls to this method.
  const std::string timerName ("Ifpack2::SupportGraph::apply");
  RCP<Time> timer = TimeMonitor::lookupCounter(timerName);
  if (timer.is_null()) {
    timer = TimeMonitor::getNewCounter(timerName);
  }

  { // Start timing here.
    Teuchos::TimeMonitor timeMon (*timer);

    TEUCHOS_TEST_FOR_EXCEPTION(
      ! isComputed(), std::runtime_error,
      "Ifpack2::SupportGraph::apply: You must call compute() to compute the "
      "incomplete factorization, before calling apply().");

    TEUCHOS_TEST_FOR_EXCEPTION(
      X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
      "Ifpack2::SupportGraph::apply: X and Y must have the same number of "
      "columns.  X has " << X.getNumVectors() << " columns, but Y has "
      << Y.getNumVectors() << " columns.");

    TEUCHOS_TEST_FOR_EXCEPTION(
      beta != STS::zero(), std::logic_error,
      "Ifpack2::SupportGraph::apply: This method does not currently work when "
      "beta != 0.");

    // If X and Y are pointing to the same memory location,
    // we need to create an auxiliary vector, Xcopy
    RCP<const MV> Xcopy;
    if (X.getLocalMV().getValues() == Y.getLocalMV().getValues()) {
      Xcopy = rcp (new MV(X));
    }
    else {
      Xcopy = rcpFromRef(X);
    }

    if (alpha != STS::one()) {
      Y.scale(alpha);
    }

    RCP<MV> Ycopy = rcpFromRef(Y);

    solver_->setB(Xcopy);
    solver_->setX(Ycopy);

    solver_->solve ();
  } // Stop timing here.

  ++NumApply_;

  // timer->totalElapsedTime() returns the total time over all timer
  // calls.  Thus, we use = instead of +=.
  ApplyTime_ = timer->totalElapsedTime();
}

int main(int argc, char *argv[]) {
  //
#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
#endif
  //
  using Teuchos::ParameterList;
  using Teuchos::RCP;
  using Teuchos::rcp;

  bool success = false;
  bool verbose = false;
  try {
    //
    // Get test parameters from command-line processor
    //
    bool proc_verbose = false;
    int frequency = -1;                  // frequency of status test output.
    std::string filename("gr_30_30.hb"); // default input filename
    double tol = 1.0e-10;                // relative residual tolerance
    int numBlocks = 30;                  // maximum number of blocks the solver can use for the Krylov subspace
    int recycleBlocks = 3;               // maximum number of blocks the solver can use for the recycle space
    int numrhs = 1;                      // number of right-hand sides to solve for
    int maxiters = -1;                   // maximum number of iterations allowed per linear system

    Teuchos::CommandLineProcessor cmdp(false,true);
    cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
    cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
    cmdp.setOption("filename",&filename,"Filename for test matrix.");
    cmdp.setOption("tol",&tol,"Relative residual tolerance used by the RCG solver.");
    cmdp.setOption("max-subspace",&numBlocks,"Maximum number of vectors in search space (not including recycle space).");
    cmdp.setOption("recycle",&recycleBlocks,"Number of vectors in recycle space.");
    cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
    cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
    if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
      return -1;
    }
    if (!verbose)
      frequency = -1;  // reset frequency if test is not verbose
    //
    // Get the problem
    //
    int MyPID;
    RCP<Epetra_CrsMatrix> A;
    RCP<Epetra_MultiVector> B, X;
    int return_val =Belos::createEpetraProblem(filename,NULL,&A,&B,&X,&MyPID);
    if(return_val != 0) return return_val;
    proc_verbose = ( verbose && (MyPID==0) );
    //
    // Solve using Belos
    //
    typedef double                          ST;
    typedef Epetra_Operator                 OP;
    typedef Epetra_MultiVector              MV;
    typedef Belos::OperatorTraits<ST,MV,OP> OPT;
    typedef Belos::MultiVecTraits<ST,MV>    MVT;
    //
    // *****Construct initial guess and right-hand sides *****
    //
    if (numrhs != 1) {
      X = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
      X->Random();
      B = rcp( new Epetra_MultiVector( A->Map(), numrhs ) );
      OPT::Apply( *A, *X, *B );
      MVT::MvInit( *X, 0.0 );
    }
    else { // initialize exact solution to be vector of ones
      MVT::MvInit( *X, 1.0 );
      OPT::Apply( *A, *X, *B );
      MVT::MvInit( *X, 0.0 );
    }
    //
    // ********Other information used by block solver***********
    // *****************(can be user specified)******************
    //
    const int NumGlobalElements = B->GlobalLength();
    if (maxiters == -1)
      maxiters = NumGlobalElements - 1; // maximum number of iterations to run
    //
    ParameterList belosList;
    belosList.set( "Maximum Iterations", maxiters );       // Maximum number of iterations allowed
    belosList.set( "Num Blocks", numBlocks);               // Maximum number of blocks in Krylov space
    belosList.set( "Num Recycled Blocks", recycleBlocks ); // Number of vectors in recycle space
    belosList.set( "Convergence Tolerance", tol );         // Relative convergence tolerance requested
    if (verbose) {
      belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
          Belos::TimingDetails + Belos::FinalSummary + Belos::StatusTestDetails );
      if (frequency > 0)
        belosList.set( "Output Frequency", frequency );
    }
    else
      belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
    //
    // Construct an unpreconditioned linear problem instance.
    //
    Belos::LinearProblem<double,MV,OP> problem( A, X, B );
    bool set = problem.setProblem();
    if (set == false) {
      if (proc_verbose)
//.........这里部分代码省略.........

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

  //Check number of arguments
  if (argc < 4) {
    std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
    std::cout <<"Usage:\n\n";
    std::cout <<"  ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ verbose\n\n";
    std::cout <<" where \n";
    std::cout <<"   int deg             - polynomial degree to be used (assumed >= 1) \n";
    std::cout <<"   int NX              - num intervals in x direction (assumed box domain, 0,1) \n";
    std::cout <<"   int NY              - num intervals in y direction (assumed box domain, 0,1) \n";
    std::cout <<"   int NZ              - num intervals in y direction (assumed box domain, 0,1) \n";
    std::cout <<"   verbose (optional)  - any character, indicates verbose output \n\n";
    exit(1);
  }
  
  // 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 > 2)
    outStream = Teuchos::rcp(&std::cout, false);
  else
    outStream = Teuchos::rcp(&bhs, false);
  
  // Save the format state of the original std::cout.
  Teuchos::oblackholestream oldFormatState;
  oldFormatState.copyfmt(std::cout);
  
  *outStream								\
    << "===============================================================================\n" \
    << "|                                                                             |\n" \
    << "|  Example: Apply Stiffness Matrix for                                        |\n" \
    << "|                   Poisson Equation on Hexahedral Mesh                       |\n" \
    << "|                                                                             |\n" \
    << "|  Questions? Contact  Pavel Bochev  ([email protected]),                    |\n" \
    << "|                      Denis Ridzal  ([email protected]),                    |\n" \
    << "|                      Kara Peterson ([email protected]).                    |\n" \
    << "|                                                                             |\n" \
    << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
    << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
    << "|                                                                             |\n" \
    << "===============================================================================\n";

  
  // ************************************ GET INPUTS **************************************
  
  int deg          = atoi(argv[1]);  // polynomial degree to use
  int NX           = atoi(argv[2]);  // num intervals in x direction (assumed box domain, 0,1)
  int NY           = atoi(argv[3]);  // num intervals in y direction (assumed box domain, 0,1)
  int NZ           = atoi(argv[4]);  // num intervals in y direction (assumed box domain, 0,1)
  

  // *********************************** CELL TOPOLOGY **********************************
  
  // Get cell topology for base hexahedron
  typedef shards::CellTopology    CellTopology;
  CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() );
  
  // Get dimensions 
  int numNodesPerElem = hex_8.getNodeCount();
  int spaceDim = hex_8.getDimension();
  
  // *********************************** GENERATE MESH ************************************
  
  *outStream << "Generating mesh ... \n\n";
  
  *outStream << "   NX" << "   NY" << "   NZ\n";
  *outStream << std::setw(5) << NX <<
    std::setw(5) << NY << std::setw(5) << NZ << "\n\n";
  
  // Print mesh information
  int numElems = NX*NY*NZ;
  int numNodes = (NX+1)*(NY+1)*(NZ+1);
  *outStream << " Number of Elements: " << numElems << " \n";
  *outStream << "    Number of Nodes: " << numNodes << " \n\n";
  
  // Cube
  double leftX = 0.0, rightX = 1.0;
  double leftY = 0.0, rightY = 1.0;
  double leftZ = 0.0, rightZ = 1.0;

  // Mesh spacing
  double hx = (rightX-leftX)/((double)NX);
  double hy = (rightY-leftY)/((double)NY);
  double hz = (rightZ-leftZ)/((double)NZ);

  // Get nodal coordinates
  FieldContainer<double> nodeCoord(numNodes, spaceDim);
  FieldContainer<int> nodeOnBoundary(numNodes);
  int inode = 0;
  for (int k=0; k<NZ+1; k++) 
    {
      for (int j=0; j<NY+1; j++) 
	{
	  for (int i=0; i<NX+1; i++) 
	    {
	      nodeCoord(inode,0) = leftX + (double)i*hx;
	      nodeCoord(inode,1) = leftY + (double)j*hy;
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
  typedef double                            ST;
  typedef Teuchos::ScalarTraits<ST>        SCT;
  typedef SCT::magnitudeType                MT;
  typedef Tpetra::MultiVector<>             MV;
  typedef Tpetra::Operator<>                OP;
  typedef Belos::MultiVecTraits<ST,MV>     MVT;
  typedef Belos::OperatorTraits<ST,MV,OP>  OPT;
  typedef Tpetra::CrsMatrix<>              CrsMatrix;
  typedef Ifpack2::Preconditioner<>        Prec;

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

  // ************************* Initialize MPI **************************
  Teuchos::oblackholestream blackhole;
  Teuchos::GlobalMPISession mpiSession (&argc, &argv, &blackhole);

  // ************** Get the default communicator and node **************
  RCP<const Teuchos::Comm<int> > comm = Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
  const int myRank = comm->getRank();

  bool verbose = true;
  bool success = true;
  bool proc_verbose = false;
  bool leftprec = true;      // left preconditioning or right.
  int frequency = -1;        // frequency of status test output.
  int numrhs = 1;            // number of right-hand sides to solve for
  int maxiters = -1;         // maximum number of iterations allowed per linear system
  std::string filename("cage4.mtx");
  MT tol = 1.0e-5;           // relative residual tolerance

  // ***************** Read the command line arguments *****************
  Teuchos::CommandLineProcessor cmdp(false,false);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("left-prec","right-prec",&leftprec,"Left preconditioning or right.");
  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
  cmdp.setOption("filename",&filename,"Filename for test matrix.  Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
  cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
    return -1;
  }
  if (!verbose)
    frequency = -1;  // reset frequency if test is not verbose
  proc_verbose = verbose && (myRank==0); /* Only print on the zero processor */

  // ************************* Get the problem *************************
  RCP<CrsMatrix> A = Tpetra::MatrixMarket::Reader<CrsMatrix>::readSparseFile(filename,comm);
  RCP<MV> B = rcp(new MV(A->getRowMap(),numrhs,false));
  RCP<MV> X = rcp(new MV(A->getRowMap(),numrhs,false));
  OPT::Apply(*A, *X, *B);
  MVT::MvInit(*X);
  MVT::MvInit(*B,1);

  // ******************** Construct preconditioner *********************
  Ifpack2::Factory factory;
  RCP<Prec> M = factory.create("RELAXATION", A.getConst());
  ParameterList ifpackParams;
  ifpackParams.set("relaxation: type","Jacobi");
  M->setParameters(ifpackParams);
  M->initialize();
  M->compute();

  // ******* Create parameter list for the Belos solver manager ********
  const int NumGlobalElements = MVT::GetGlobalLength(*B);
  if (maxiters == -1)
    maxiters = NumGlobalElements - 1; // maximum number of iterations to run
  //
  ParameterList belosList;
  belosList.set( "Maximum Iterations", maxiters );       // Maximum number of iterations allowed
  belosList.set( "Convergence Tolerance", tol );         // Relative convergence tolerance requested
  if (numrhs > 1) {
    belosList.set( "Show Maximum Residual Norm Only", true );  // Show only the maximum residual norm
  }
  if (verbose) {
    belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
		   Belos::TimingDetails + Belos::StatusTestDetails );
    if (frequency > 0)
      belosList.set( "Output Frequency", frequency );
  }
  else
    belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + Belos::FinalSummary );

  // ************ Construct a preconditioned linear problem ************
  RCP<Belos::LinearProblem<double,MV,OP> > problem
    = rcp( new Belos::LinearProblem<double,MV,OP>( A, X, B ) );
  if (leftprec) {
    problem->setLeftPrec( M );
  }
  else {
    problem->setRightPrec( M );
  }
  bool set = problem->setProblem();
  if (set == false) {
    if (proc_verbose)
      std::cout << std::endl << "ERROR:  Belos::LinearProblem failed to set up correctly!" << std::endl;
    return -1;
//.........这里部分代码省略.........

void Piro::RythmosSolver<Scalar>::initialize(
#endif
    const Teuchos::RCP<Teuchos::ParameterList> &appParams,
    const Teuchos::RCP< Thyra::ModelEvaluator<Scalar> > &in_model,
    const Teuchos::RCP<Rythmos::IntegrationObserverBase<Scalar> > &observer)
{

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

    // set some internals
    model = in_model;
    num_p = in_model->Np();
    num_g = in_model->Ng();

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


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

        {
            const std::string verbosity = rythmosPL->get("Verbosity Level", "VERB_DEFAULT");
            if      (verbosity == "VERB_NONE")    solnVerbLevel = Teuchos::VERB_NONE;
            else if (verbosity == "VERB_DEFAULT") solnVerbLevel = Teuchos::VERB_DEFAULT;
            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;
            else TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,"Unknown verbosity option specified in Piro_RythmosSolver.");
        }

        t_initial = rythmosPL->get("Initial Time", 0.0);
        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<Scalar> > rythmosTimeStepSolver =
                Rythmos::timeStepNonlinearSolver<Scalar>();
            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 HAVE_PIRO_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 == "Forward Euler") {
            fwdStateStepper = Rythmos::forwardEulerStepper<Scalar> (model);
            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());

        }
        else {
            // first (before failing) check to see if the user has added stepper factory
            typename std::map<std::string,Teuchos::RCP<Piro::RythmosStepperFactory<Scalar> > >::const_iterator
            stepFactItr = stepperFactories.find(stepperType);
            if(stepFactItr!=stepperFactories.end()) {
                // the user has added it, hot dog lets build a new stepper!
//.........这里部分代码省略.........

Teuchos::RCP< std::vector< Teuchos::RCP<PHX::Evaluator<panzer::Traits> > > > 
user_app::MyModelFactory<EvalT>::
buildClosureModels(const std::string& model_id,
		   const Teuchos::ParameterList& models, 
		   const panzer::FieldLayoutLibrary& fl,
		   const Teuchos::RCP<panzer::IntegrationRule>& ir,
		   const Teuchos::ParameterList& default_params,
		   const Teuchos::ParameterList& user_data,
		   const Teuchos::RCP<panzer::GlobalData>& global_data,
		   PHX::FieldManager<panzer::Traits>& fm) const
{

  using std::string;
  using std::vector;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::ParameterList;
  using PHX::Evaluator;

  RCP< vector< RCP<Evaluator<panzer::Traits> > > > evaluators = 
    rcp(new vector< RCP<Evaluator<panzer::Traits> > > );

  if (!models.isSublist(model_id)) {
    models.print(std::cout);
    std::stringstream msg;
    msg << "Falied to find requested model, \"" << model_id 
	<< "\", for equation set:\n" << std::endl;
    TEUCHOS_TEST_FOR_EXCEPTION(!models.isSublist(model_id), std::logic_error, msg.str());
  }

  std::vector<Teuchos::RCP<const panzer::PureBasis> > bases;
  fl.uniqueBases(bases);

  const ParameterList& my_models = models.sublist(model_id);

  for (ParameterList::ConstIterator model_it = my_models.begin(); 
       model_it != my_models.end(); ++model_it) {
    
    bool found = false;
    
    const std::string key = model_it->first;
    ParameterList input;
    const Teuchos::ParameterEntry& entry = model_it->second;
    const ParameterList& plist = Teuchos::getValue<Teuchos::ParameterList>(entry);

    #ifdef HAVE_STOKHOS
    if (plist.isType<double>("Value") && plist.isType<double>("UQ") 
                           && plist.isParameter("Expansion")
                           && (typeid(EvalT)==typeid(panzer::Traits::SGResidual) || 
                               typeid(EvalT)==typeid(panzer::Traits::SGJacobian)) ) {
      { // at IP
	input.set("Name", key);
	input.set("Value", plist.get<double>("Value"));
	input.set("UQ", plist.get<double>("UQ"));
	input.set("Expansion", plist.get<Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > >("Expansion"));
	input.set("Data Layout", ir->dl_scalar);
	RCP< Evaluator<panzer::Traits> > e = 
	  rcp(new user_app::ConstantModel<EvalT,panzer::Traits>(input));
	evaluators->push_back(e);
      }
      
      for (std::vector<Teuchos::RCP<const panzer::PureBasis> >::const_iterator basis_itr = bases.begin();
	   basis_itr != bases.end(); ++basis_itr) { // at BASIS
	input.set("Name", key);
	input.set("Value", plist.get<double>("Value"));
	input.set("UQ", plist.get<double>("UQ"));
	input.set("Expansion", plist.get<Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > >("Expansion"));
	Teuchos::RCP<const panzer::BasisIRLayout> basis = basisIRLayout(*basis_itr,*ir);
	input.set("Data Layout", basis->functional);
	RCP< Evaluator<panzer::Traits> > e = 
	  rcp(new user_app::ConstantModel<EvalT,panzer::Traits>(input));
	evaluators->push_back(e);
      }
      found = true;
    }
    else 
    #endif
    if (plist.isType<std::string>("Type")) {
      
      if (plist.get<std::string>("Type") == "Parameter") {
	{ // at IP
	  RCP< Evaluator<panzer::Traits> > e = 
	    rcp(new panzer::Parameter<EvalT,panzer::Traits>(key,ir->dl_scalar,plist.get<double>("Value"),*global_data->pl));
	  evaluators->push_back(e);
	}
	
	for (std::vector<Teuchos::RCP<const panzer::PureBasis> >::const_iterator basis_itr = bases.begin();
	   basis_itr != bases.end(); ++basis_itr) { // at BASIS
	  Teuchos::RCP<const panzer::BasisIRLayout> basis = basisIRLayout(*basis_itr,*ir);
	  RCP< Evaluator<panzer::Traits> > e = 
	    rcp(new panzer::Parameter<EvalT,panzer::Traits>(key,basis->functional,plist.get<double>("Value"),*global_data->pl));
	  evaluators->push_back(e);
	}
	
	found = true;
      }
  
    }
    else if (plist.isType<double>("Value")) {
      { // at IP
//.........这里部分代码省略.........

  void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
  initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
    using Teuchos::rcp_dynamic_cast;

    // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
    typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>                     XpMap;
    typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node>      XpOp;
    typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>       XpThyUtils;
    typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>        XpCrsMat;
    typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
    typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>           XpMat;
    typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>      XpMultVec;
    typedef Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node>      XpMultVecDouble;
    typedef Thyra::LinearOpBase<Scalar>                                      ThyLinOpBase;
#ifdef HAVE_MUELU_TPETRA
    typedef MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> MueTpOp;
    typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>      TpOp;
    typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
#endif

    // Check precondition
    TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
    TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
    TEUCHOS_ASSERT(prec);

    // Create a copy, as we may remove some things from the list
    ParameterList paramList = *paramList_;

    // Retrieve wrapped concrete Xpetra matrix from FwdOp
    const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
    TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));

    // Check whether it is Epetra/Tpetra
    bool bIsEpetra  = XpThyUtils::isEpetra(fwdOp);
    bool bIsTpetra  = XpThyUtils::isTpetra(fwdOp);
    bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
    TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true  && bIsTpetra == true));
    TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
    TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);

    RCP<XpMat> A = Teuchos::null;
    if(bIsBlocked) {
      Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
          Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));

      TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);

      Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));

      RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));

      // MueLu needs a non-const object as input
      RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));

      // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
      RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));

      RCP<const XpMap> rowmap00 = A00->getRowMap();
      RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();

      // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
      RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));

      // save blocked matrix
      A = bMat;
    } else {
      RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));

      // MueLu needs a non-const object as input
      RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
      TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));

      // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
      A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
    }
    TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));

    // Retrieve concrete preconditioner object
    const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
    TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));

    // extract preconditioner operator
    RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
    thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);

    // Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy
    RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = Teuchos::null;

    // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
    // rebuild preconditioner if startingOver == true
    // reuse preconditioner if startingOver == false
    const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");

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