C++ FactoryManager类代码示例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

//.........这里部分代码省略.........
    // build transfer operators
    RCP<TentativePFactory> TentPFact = rcp(new TentativePFactory());

    /*TentPFact->setStridingData(stridingInfo);
    TentPFact->setStridedBlockId(stridedBlockId);
    TentPFact->setDomainMapOffset(offset);*/

    RCP<CoarseMapFactory> coarseMapFact = Teuchos::rcp(new CoarseMapFactory());
    coarseMapFact->setStridingData(stridingInfo);
    coarseMapFact->setStridedBlockId(stridedBlockId);
    coarseMapFact->setDomainMapOffset(offset);

    RCP<SaPFactory> Pfact  = rcp( new SaPFactory() );
    //RCP<PgPFactory> Pfact  = rcp( new PgPFactory() );
    //RCP<TentativePFactory> Pfact  = rcp( new TentativePFactory() );
    RCP<Factory>   Rfact  = rcp( new TransPFactory() );
    //RCP<Factory>   Rfact  = rcp( new GenericRFactory() );

    // RAP Factory
    RCP<RAPFactory> Acfact = rcp( new RAPFactory() );
    Acfact->setVerbLevel(Teuchos::VERB_HIGH);

    // register aggregation export factory in RAPFactory
    //RCP<MueLu::AggregationExportFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node, LocalMatOps> > aggExpFact = rcp(new MueLu::AggregationExportFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node, LocalMatOps>());
    //aggExpFact->SetParameter("Output filename","aggs_level%LEVELID_proc%PROCID.out");
    //Acfact->AddTransferFactory(aggExpFact);

    // build level smoothers
    RCP<SmootherPrototype> smooProto;
    std::string ifpackType;
    Teuchos::ParameterList ifpackList;
    ifpackList.set("relaxation: sweeps", (LO) 1);
    ifpackList.set("relaxation: damping factor", (SC) 1.0); // 0.7
    ifpackType = "RELAXATION";
    ifpackList.set("relaxation: type", "Symmetric Gauss-Seidel");

    smooProto = Teuchos::rcp( new TrilinosSmoother(ifpackType, ifpackList) );
    RCP<SmootherFactory> SmooFact;
    if (maxLevels > 1)
        SmooFact = rcp( new SmootherFactory(smooProto) );

    // create coarsest smoother
    RCP<SmootherPrototype> coarsestSmooProto;
    std::string type = "";
    Teuchos::ParameterList coarsestSmooList;
#if defined(HAVE_AMESOS_SUPERLU)
    coarsestSmooProto = Teuchos::rcp( new DirectSolver("Superlu", coarsestSmooList) );
#else
    coarsestSmooProto = Teuchos::rcp( new DirectSolver("Klu", coarsestSmooList) );
#endif
    RCP<SmootherFactory> coarsestSmooFact;
    coarsestSmooFact = rcp(new SmootherFactory(coarsestSmooProto, Teuchos::null));

    FactoryManager M;
    //M.SetFactory("Graph", dropFact);
    //M.SetFactory("UnAmalgamationInfo", dropFact);
    M.SetFactory("Aggregates", CoupledAggFact);
    M.SetFactory("Ptent", TentPFact);
    M.SetFactory("P", Pfact);
    M.SetFactory("R", Rfact);
    M.SetFactory("A", Acfact);
    M.SetFactory("Smoother", SmooFact);
    M.SetFactory("CoarseSolver", coarsestSmooFact);
    M.SetFactory("CoarseMap", coarseMapFact);

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

    RCP<Vector> xLsg = VectorFactory::Build(map);

    // Use AMG directly as an iterative method
    {
        xLsg->putScalar( (SC) 0.0);

        H->Iterate(*xRhs,10,*xLsg);

        //xLsg->describe(*out,Teuchos::VERB_EXTREME);
    }

    //
    // Solve Ax = b using AMG as a preconditioner in AztecOO
    //
    {
        RCP<Epetra_Vector> X = rcp(new Epetra_Vector(epv->Map()));
        X->PutScalar(0.0);
        Epetra_LinearProblem epetraProblem(epA.get(), X.get(), epv.get());

        AztecOO aztecSolver(epetraProblem);
        aztecSolver.SetAztecOption(AZ_solver, AZ_cg);

        MueLu::EpetraOperator aztecPrec(H);
        aztecSolver.SetPrecOperator(&aztecPrec);

        int maxIts = 50;
        double tol = 1e-8;

        aztecSolver.Iterate(maxIts, tol);
    }

    return xLsg;
}

//.........这里部分代码省略.........
  RCP<Hierarchy> H = rcp ( new Hierarchy() );
  H->setDefaultVerbLevel(Teuchos::VERB_HIGH);
  H->SetMaxCoarseSize(maxCoarseSize);
  H->setlib(Xpetra::UseEpetra);

  // build finest Level
  RCP<MueLu::Level> Finest = H->GetLevel();
  Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
  Finest->setlib(Xpetra::UseEpetra);
  Finest->Set("A",A);
  Finest->Set("Nullspace",nullspace);

  // create factories for transfer operators
  RCP<TentativePFactory> PFact = Teuchos::rcp(new TentativePFactory());
  RCP<TransPFactory>     RFact = Teuchos::rcp(new TransPFactory());
  RFact->SetFactory("P", PFact);

  // build level smoothers
  // use symmetric Gauss-Seidel both for fine and coarse level smoother
  RCP<SmootherPrototype> smooProto;
  std::string ifpackType;
  Teuchos::ParameterList ifpackList;
  ifpackList.set("relaxation: sweeps", (LO) 1);
  ifpackList.set("relaxation: damping factor", (SC) 1.0);
  ifpackType = "RELAXATION";
  ifpackList.set("relaxation: type", "Symmetric Gauss-Seidel");

  smooProto = Teuchos::rcp( new TrilinosSmoother(ifpackType, ifpackList) );
  RCP<SmootherFactory> SmooFact;
  if (maxLevels > 1)
    SmooFact = rcp( new SmootherFactory(smooProto) );

  // design multigrid hierarchy
  FactoryManager M;
  M.SetFactory("P", PFact);
  M.SetFactory("R", RFact);
  M.SetFactory("Nullspace", PFact);
  M.SetFactory("Smoother", SmooFact);
  M.SetFactory("CoarseSolver", SmooFact);

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

  comm->barrier();
  tm = Teuchos::null;

  // =========================================================================
  // System solution (Ax = b)
  // =========================================================================

    //
    // generate exact solution using a direct solver
    //
    RCP<Epetra_Vector> exactLsgVec = rcp(new Epetra_Vector(X->Map()));
    {
      fancyout << "========================================================\nCalculate exact solution." << std::endl;
      tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 3 - direct solve")));
      exactLsgVec->PutScalar(0.0);
      exactLsgVec->Update(1.0,*X,1.0);
      Epetra_LinearProblem epetraProblem(epA.get(), exactLsgVec.get(), B.get());

      Amesos amesosFactory;
      RCP<Amesos_BaseSolver> rcp_directSolver = Teuchos::rcp(amesosFactory.Create("Amesos_Klu", epetraProblem));
      rcp_directSolver->SymbolicFactorization();
      rcp_directSolver->NumericFactorization();
      rcp_directSolver->Solve();

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

  ///////////////////////////////////////// define CoalesceDropFactory and Aggregation for A11
  // set up amalgamation for A11. Note: we're using a default null space factory (Teuchos::null)
  RCP<AmalgamationFactory> amalgFact11 = rcp(new AmalgamationFactory());
  amalgFact11->SetFactory("A", A11Fact);

  amalgFact11->setDefaultVerbLevel(Teuchos::VERB_EXTREME);
  RCP<CoalesceDropFactory> dropFact11 = rcp(new CoalesceDropFactory());
  dropFact11->SetFactory("A", A11Fact);
  dropFact11->SetFactory("UnAmalgamationInfo", amalgFact11);
  dropFact11->setDefaultVerbLevel(Teuchos::VERB_EXTREME);
  RCP<UncoupledAggregationFactory> CoupledAggFact11 = rcp(new UncoupledAggregationFactory());
  CoupledAggFact11->SetFactory("Graph", dropFact11);
  CoupledAggFact11->SetMinNodesPerAggregate(9);
  CoupledAggFact11->SetMaxNeighAlreadySelected(2);
  CoupledAggFact11->SetOrdering(MueLu::AggOptions::NATURAL);
  //CoupledAggFact11->SetPhase3AggCreation(0.5);

  ///////////////////////////////////////// define transfer ops for A11
#if 0
  // use PG-AMG
  RCP<PgPFactory> P11Fact = rcp(new PgPFactory());

  RCP<GenericRFactory> R11Fact = rcp(new GenericRFactory());
  Teuchos::RCP<NullspaceFactory> nspFact11 = Teuchos::rcp(new NullspaceFactory("Nullspace1",P11tentFact));

  Teuchos::RCP<NullspaceFactory> nspFact11 = Teuchos::rcp(new NullspaceFactory("Nullspace1"));

  RCP<CoarseMapFactory> coarseMapFact11 = Teuchos::rcp(new CoarseMapFactory());
  coarseMapFact11->setStridingData(stridingInfo);
  coarseMapFact11->setStridedBlockId(0);

  //////////////////////////////// define factory manager for (1,1) block
  RCP<FactoryManager> M11 = rcp(new FactoryManager());
  M11->SetFactory("A", A11Fact);
  M11->SetFactory("P", P11Fact);
  M11->SetFactory("R", R11Fact);
  M11->SetFactory("Aggregates", CoupledAggFact11);
  M11->SetFactory("UnAmalgamationInfo", amalgFact11);
  M11->SetFactory("Nullspace", nspFact11);
  // M11->SetFactory("Ptent", P11tentFact);
  M11->SetFactory("CoarseMap", coarseMapFact11);
#else
  RCP<TentativePFactory> P11Fact = rcp(new TentativePFactory());

  RCP<TransPFactory> R11Fact = rcp(new TransPFactory());

  Teuchos::RCP<NullspaceFactory> nspFact11 = Teuchos::rcp(new NullspaceFactory("Nullspace1"));
  nspFact11->SetFactory("Nullspace1",P11Fact);

  RCP<CoarseMapFactory> coarseMapFact11 = Teuchos::rcp(new CoarseMapFactory());
  coarseMapFact11->setStridingData(stridingInfo);
  coarseMapFact11->setStridedBlockId(0);

  //////////////////////////////// define factory manager for (1,1) block
  RCP<FactoryManager> M11 = rcp(new FactoryManager());
  M11->SetFactory("A", A11Fact);
  M11->SetFactory("P", P11Fact);
  M11->SetFactory("R", R11Fact);
  M11->SetFactory("Aggregates", CoupledAggFact11);
  M11->SetFactory("UnAmalgamationInfo", amalgFact11);
  M11->SetFactory("Nullspace", nspFact11);
  // M11->SetFactory("Ptent", P11Fact);
  M11->SetFactory("CoarseMap", coarseMapFact11);
#endif
  M11->SetIgnoreUserData(true);               // always use data from factories defined in factory manager

//.........这里部分代码省略.........
            RepartitionFact2->SetParameterList(paramList);
        }
        RepartitionFact2->SetFactory("A", rebA22Fact);
        RepartitionFact2->SetFactory("Partition", repInterface2); // this is not valid

        ////////////////////////////////////////// build non-rebalanced matrix blocks
        // build factories for transfer operator P(1,1) and R(1,1)
        RCP<AmalgamationFactory> amalgFact11 = rcp(new AmalgamationFactory());
        amalgFact11->SetFactory("A", A11Fact);
        amalgFact11->setDefaultVerbLevel(Teuchos::VERB_EXTREME);

        RCP<CoalesceDropFactory> dropFact11 = rcp(new CoalesceDropFactory());
        dropFact11->SetFactory("A", A11Fact);
        dropFact11->SetFactory("UnAmalgamationInfo", amalgFact11);
        dropFact11->setDefaultVerbLevel(Teuchos::VERB_EXTREME);

        RCP<UncoupledAggregationFactory> UncoupledAggFact11 = rcp(new UncoupledAggregationFactory());
        UncoupledAggFact11->SetFactory("Graph", dropFact11);
        UncoupledAggFact11->SetMinNodesPerAggregate(9);
        UncoupledAggFact11->SetMaxNeighAlreadySelected(2);
        UncoupledAggFact11->SetOrdering("natural");

        RCP<CoarseMapFactory> coarseMapFact11 = Teuchos::rcp(new CoarseMapFactory());
        coarseMapFact11->setStridingData(stridingInfo);
        coarseMapFact11->setStridedBlockId(0);

        RCP<TentativePFactory> P11Fact = rcp(new TentativePFactory());
        RCP<TransPFactory> R11Fact = rcp(new TransPFactory());

        Teuchos::RCP<NullspaceFactory> nspFact11 = Teuchos::rcp(new NullspaceFactory("Nullspace1"));
        nspFact11->SetFactory("Nullspace1",P11Fact); // pick "Nullspace1" from Finest level

        //////////////////////////////// define factory manager for (1,1) block
        RCP<FactoryManager> M11 = rcp(new FactoryManager());
        M11->SetFactory("A", A11Fact);  // rebalanced fine-level block operator
        M11->SetFactory("P", P11Fact);  // non-rebalanced transfer operator block P(1,1)
        M11->SetFactory("R", R11Fact);  // non-rebalanced transfer operator block R(1,1)
        M11->SetFactory("Aggregates", UncoupledAggFact11);
        M11->SetFactory("Graph", dropFact11);
        M11->SetFactory("DofsPerNode", dropFact11);
        M11->SetFactory("UnAmalgamationInfo", amalgFact11);
        M11->SetFactory("Nullspace", nspFact11); // TODO check me?
        M11->SetFactory("CoarseMap", coarseMapFact11);
        M11->SetIgnoreUserData(true);               // always use data from factories defined in factory manager

        ////////////////////////////////////////// build non-rebalanced matrix blocks
        // build factories for transfer operator P(2,2) and R(2,2)
        RCP<AmalgamationFactory> amalgFact22 = rcp(new AmalgamationFactory());
        RCP<TentativePFactory> P22Fact = rcp(new TentativePFactory());
        RCP<TransPFactory> R22Fact = rcp(new TransPFactory());

        // connect null space and tentative PFactory
        Teuchos::RCP<NullspaceFactory> nspFact22 = Teuchos::rcp(new NullspaceFactory("Nullspace2"));
        nspFact22->SetFactory("Nullspace2", P22Fact); // define null space generated by P22Fact as null space for coarse level (non-rebalanced)

        RCP<CoarseMapFactory> coarseMapFact22 = Teuchos::rcp(new CoarseMapFactory());
        coarseMapFact22->setStridingData(stridingInfo);
        coarseMapFact22->setStridedBlockId(1);

        //////////////////////////////// define factory manager for (2,2) block
        RCP<FactoryManager> M22 = rcp(new FactoryManager());
        M22->SetFactory("A", A22Fact); // rebalanced fine-level block operator
        M22->SetFactory("P", P22Fact); // non-rebalanced transfer operator P(2,2)
        M22->SetFactory("R", R22Fact); // non-rebalanced transfer operator R(2,2)
        M22->SetFactory("Aggregates", UncoupledAggFact11); // aggregates from block (1,1)
        M22->SetFactory("Nullspace", nspFact22);

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

  using Teuchos::RCP;
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);

  typedef Teuchos::ScalarTraits<SC> ST;

  bool success = false;
  bool verbose = true;
  try {
    RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
    //
    // Parameters
    //

    Teuchos::CommandLineProcessor clp(false);
    Galeri::Xpetra::Parameters<GO> matrixParameters(clp, 8748);
    Xpetra::Parameters xpetraParameters(clp);

    bool optRecycling           = true;  clp.setOption("recycling",             "no-recycling",             &optRecycling,           "Enable recycling of the multigrid preconditioner");

    /* DO NOT WORK YET
       bool optRecyclingRAPpattern = true;  clp.setOption("recycling-rap-pattern", "no-recycling-rap-pattern", &optRecyclingRAPpattern, "Enable recycling of Ac=RAP pattern");
       bool optRecyclingAPpattern  = false; clp.setOption("recycling-ap-pattern",  "no-recycling-ap-pattern",  &optRecyclingAPpattern,  "Enable recycling of AP pattern");
       */
    bool optRecyclingRAPpattern = false;
    bool optRecyclingAPpattern  = false;

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

    // option dependencies
    if (optRecycling == false) {
      optRecyclingRAPpattern = false;
      optRecyclingAPpattern  = false;
    }

    //
    // Construct the problems
    //

    RCP<const Map> map = MapFactory::Build(xpetraParameters.GetLib(), matrixParameters.GetNumGlobalElements(), 0, comm);
    Teuchos::RCP<Galeri::Xpetra::Problem<Map,CrsMatrixWrap,MultiVector> > Pr =
      Galeri::Xpetra::BuildProblem<SC,LO,GO,Map,CrsMatrixWrap,MultiVector>(matrixParameters.GetMatrixType(), map, matrixParameters.GetParameterList());
    RCP<Matrix> A1 = Pr->BuildMatrix();

    RCP<Matrix> A2 = Pr->BuildMatrix(); // TODO: generate another problem would be more meaningful (ex: scale A1)

    //
    // First solve
    //

    FactoryManager M;

    Hierarchy H(A1);

    if (optRecycling) {
      // Configuring "Keep" options before running the first setup.

      // Note: "Keep" flags should not be set on the default factories provided by the FactoryManager because in the current
      // implementation of FactoryManager, those factories are freed and reallocated between Hierarchy::Setup() calls (see FactoryManager::Clean() and Hierarchy::Setup()).
      // So we define our own factories here.

      // AGGREGATES:
      // Note: aggregates are only used to build Ptent, so it is not useful to keep them. Keeping Ptent is enough.
      // RCP<Factory> AggFact   = rcp(new CoupledAggregationFactory());
      // M.SetFactory("Aggregates", AggFact);
      // H.Keep("Aggregates", AggFact.get());

      // PTENT:
      RCP<Factory> PtentFact = rcp(new TentativePFactory());
      M.SetFactory("Ptent", PtentFact);
      H.Keep("P",           PtentFact.get());
    }

    RCP<Factory> AcFact = rcp(new RAPFactory());
    M.SetFactory("A", AcFact);

    if (optRecyclingRAPpattern) {
      H.Keep("RAP graph", AcFact.get());
    }
    if (optRecyclingAPpattern) {
      H.Keep("AP graph", AcFact.get());
    }
    //

    H.Setup(M);

    {
      RCP<Vector> X = VectorFactory::Build(map);
      RCP<Vector> B = VectorFactory::Build(map);

      X->putScalar((Scalar) 0.0);
      B->setSeed(846930886); B->randomize();

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

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

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

  //
  // MPI initialization
  //

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

  bool success = false;
  bool verbose = true;
  try {
    RCP< const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();

    //
    // Process command line arguments
    //
    Teuchos::CommandLineProcessor  clp(false);
    Galeri::Xpetra::Parameters<GO> matrixParameters(clp, 81); // manage parameters of the test case
    Xpetra::Parameters             xpetraParameters(clp);     // manage parameters of xpetra

    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;
      default:;
    }

    if (comm->getRank() == 0) std::cout << xpetraParameters << matrixParameters;

    //
    // Setup test case (Ax = b)
    //

    // Linear Algebra Library
    Xpetra::UnderlyingLib lib = xpetraParameters.GetLib();

    // Distribution
    RCP<const Map> map = MapFactory::Build(lib, matrixParameters.GetNumGlobalElements(), 0, comm);

    // Matrix
    RCP<Galeri::Xpetra::Problem<Map,CrsMatrixWrap,MultiVector> > Pr =
      Galeri::Xpetra::BuildProblem<SC, LO, GO, Map, CrsMatrixWrap, MultiVector>(matrixParameters.GetMatrixType(), map, matrixParameters.GetParameterList());
    RCP<Matrix> A = Pr->BuildMatrix();

    // User defined nullspace
    RCP<MultiVector> nullSpace = VectorFactory::Build(map,1); nullSpace->putScalar((SC) 1.0);

    // Define B
    RCP<Vector> X = VectorFactory::Build(map,1);
    RCP<Vector> B = VectorFactory::Build(map,1);
    X->setSeed(846930886);
    X->randomize();
    A->apply(*X, *B, Teuchos::NO_TRANS, (SC)1.0, (SC)0.0);

    // X = 0
    X->putScalar((SC) 0.0);

    //
    // Create a multigrid configuration
    //

    // Transfer operators
    RCP<TentativePFactory> TentativePFact = rcp( new TentativePFactory() );
    RCP<SaPFactory>        SaPFact        = rcp( new SaPFactory() );
    RCP<TransPFactory>     RFact          = rcp( new TransPFactory());

    FactoryManager M;
    M.SetFactory("Ptent", TentativePFact);
    M.SetFactory("P",     SaPFact);
    M.SetFactory("R",     RFact);

    M.SetFactory("Smoother", Teuchos::null);      //skips smoother setup
    M.SetFactory("CoarseSolver", Teuchos::null);  //skips coarsest solve setup

    //
    // Multigrid setup phase
    //

    int startLevel = 0;
    int maxLevels = 10;

    std::cout << "=============== Setup transfers only ====================" << std::endl;

    Hierarchy H;
    H.SetDefaultVerbLevel(MueLu::Medium);

    RCP<Level> finestLevel = H.GetLevel();
    finestLevel->Set("A", A);
    finestLevel->Set("Nullspace", nullSpace);

    // Indicate which Hierarchy operators we want to keep
    H.Keep("P", SaPFact.get());  //SaPFact is the generating factory for P.
    H.Keep("R", RFact.get());    //RFact is the generating factory for R.
    H.Keep("Ptent", TentativePFact.get());  //SaPFact is the generating factory for P.
//.........这里部分代码省略.........

  TEUCHOS_UNIT_TEST(GenericRFactory, SymmetricProblem)
  {
    out << "version: " << MueLu::Version() << std::endl;
    RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();

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

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

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

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

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

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

    RCP<SaPFactory>         Pfact = rcp( new SaPFactory());
    RCP<Factory>           Rfact = rcp( new GenericRFactory() );
    H->SetMaxCoarseSize(1);

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

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

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

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

    RCP<Level> coarseLevel = H->GetLevel(1);
    RCP<Matrix> P1 = coarseLevel->Get< RCP<Matrix> >("P");
    RCP<Matrix> R1 = coarseLevel->Get< RCP<Matrix> >("R");
    RCP<Level> coarseLevel2 = H->GetLevel(2);
    RCP<Matrix> P2 = coarseLevel2->Get< RCP<Matrix> >("P");
    RCP<Matrix> R2 = coarseLevel2->Get< RCP<Matrix> >("R");

    TEST_EQUALITY(Finest->IsAvailable("PreSmoother"), true);
    TEST_EQUALITY(Finest->IsAvailable("PostSmoother"), true);
    TEST_EQUALITY(coarseLevel->IsAvailable("PreSmoother"), true);
    TEST_EQUALITY(coarseLevel->IsAvailable("PostSmoother"), true);
    TEST_EQUALITY(coarseLevel2->IsAvailable("PreSmoother"), true);
    TEST_EQUALITY(coarseLevel2->IsAvailable("PostSmoother"), false);

    // test some basic multgrid data
    TEST_EQUALITY(P1->getGlobalNumEntries(), R1->getGlobalNumEntries());
    TEST_EQUALITY(P1->getGlobalNumRows(), R1->getGlobalNumCols());
    TEST_EQUALITY(P1->getGlobalNumCols(), R1->getGlobalNumRows());
    TEST_EQUALITY(P2->getGlobalNumEntries(), R2->getGlobalNumEntries());
    TEST_EQUALITY(P2->getGlobalNumRows(), R2->getGlobalNumCols());
    TEST_EQUALITY(P2->getGlobalNumCols(), R2->getGlobalNumRows());


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

    // since A is chosen symmetric, it is P^T = R
    // check P^T * P = R * P
    // note: the Epetra matrix-matrix multiplication using implicit transpose is buggy in parallel case
    //       (for multiplication of a square matrix with a rectangular matrix)
    //       however it seems to work for two rectangular matrices
    Teuchos::RCP<Xpetra::Matrix<Scalar,LO,GO> > RP = MueLu::Utils<Scalar,LO,GO>::Multiply(*R1,false,*P1,false,out);
    Teuchos::RCP<Xpetra::Matrix<Scalar,LO,GO> > PtP = MueLu::Utils<Scalar,LO,GO>::Multiply(*P1,true,*P1,false,out);

    RCP<Vector> x = VectorFactory::Build(RP->getDomainMap());
    RCP<Vector> bRP  = VectorFactory::Build(RP->getRangeMap());
    RCP<Vector> bPtP = VectorFactory::Build(PtP->getRangeMap());
//.........这里部分代码省略.........

  TEUCHOS_UNIT_TEST(MultiVectorTransferFactory, ThreeLevels)
  {
    out << "version: " << MueLu::Version() << std::endl;

    out << "Tests usage on a three-level hierarchy." << std::endl;

    GO nx = 199;
    RCP<Matrix> A = TestHelpers::TestFactory<SC, LO, GO, NO, LMO>::Build1DPoisson(nx);


    // Set up three level hierarchy.
    RCP<Hierarchy> H = rcp( new Hierarchy() );
    H->setDefaultVerbLevel(Teuchos::VERB_HIGH);

    RCP<Level> fineLevel = H->GetLevel();
    fineLevel->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    fineLevel->Set("A",A);                       // set fine level matrix
    RCP<MultiVector> nullSpace = MultiVectorFactory::Build(A->getRowMap(),1);
    nullSpace->putScalar( (SC) 1.0);
    fineLevel->Set("Nullspace",nullSpace);       // set null space information for finest level

    RCP<CoupledAggregationFactory> CoupledAggFact = rcp(new CoupledAggregationFactory());
    CoupledAggFact->SetMinNodesPerAggregate(3);
    CoupledAggFact->SetMaxNeighAlreadySelected(0);
    CoupledAggFact->SetOrdering(MueLu::AggOptions::NATURAL);
    CoupledAggFact->SetPhase3AggCreation(0.5);

    RCP<TentativePFactory> PFact  = rcp(new TentativePFactory()); //just using plain aggregation
    RCP<Factory>           RFact  = rcp(new TransPFactory());
    RCP<RAPFactory>        AcFact = rcp(new RAPFactory());
    H->SetMaxCoarseSize(1);

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

    FactoryManager M;
    M.SetFactory("Aggregates", CoupledAggFact);
    M.SetFactory("P", PFact);
    M.SetFactory("Ptent", PFact); // for nullspace
    M.SetFactory("R", RFact);
    M.SetFactory("A", AcFact);
    M.SetFactory("Smoother", SmooFact);
    M.SetFactory("CoarseSolver", SmooFact); // This line avoid dependency to Amesos/Amesos2 for this test.

    //set up the transfer factory
    RCP<MultiVector> fineOnes = MultiVectorFactory::Build(A->getRowMap(),1);
    fineOnes->putScalar(1.0);
    fineLevel->Set("onesVector",fineOnes);
    RCP<MueLu::MultiVectorTransferFactory<SC, LO, GO, NO, LMO> > mvtf = rcp(new MueLu::MultiVectorTransferFactory<SC, LO, GO, NO, LMO>("onesVector"));
    mvtf->SetFactory("R",RFact);
    M.SetFactory("onesVector",mvtf);
    AcFact->AddTransferFactory(mvtf);

    int maxLevels = 3;
    H->Setup(M, 0, maxLevels);

/*
    //FIXME we probably need to do some requests....
    coarseLevel.Request("onesVector",mvtf.get());
    coarseLevel.Request("R",RFact.get());
    coarseLevel.Request("P",TentativePFact.get());
*/

/*
    RCP<MultiVector> coarseOnes = coarseLevel.Get<RCP<MultiVector> >("onesVector",mvtf.get());
    Teuchos::Array<Teuchos::ScalarTraits<SC>::magnitudeType> vn(1);
    coarseOnes->norm2(vn);

    TEST_FLOATING_EQUALITY(vn[0]*vn[0],((SC)fineOnes->getGlobalLength()),1e-12);
*/
  } // ThreeLevels

//.........这里部分代码省略.........
    RCP<SubBlockAFactory> A11Fact = Teuchos::rcp(new SubBlockAFactory());
    A11Fact->SetFactory("A",MueLu::NoFactory::getRCP());
    A11Fact->SetParameter("block row",Teuchos::ParameterEntry(0));
    A11Fact->SetParameter("block col",Teuchos::ParameterEntry(0));
    RCP<SubBlockAFactory> A22Fact = Teuchos::rcp(new SubBlockAFactory());
    A22Fact->SetFactory("A",MueLu::NoFactory::getRCP());
    A22Fact->SetParameter("block row",Teuchos::ParameterEntry(1));
    A22Fact->SetParameter("block col",Teuchos::ParameterEntry(1));

    RCP<TentativePFactory> P11Fact = rcp(new TentativePFactory());
    RCP<TransPFactory> R11Fact = rcp(new TransPFactory());

    RCP<TentativePFactory> P22TentFact = rcp(new TentativePFactory());
    RCP<PgPFactory> P22Fact = rcp(new PgPFactory());
    RCP<GenericRFactory> R22Fact = rcp(new GenericRFactory());

    std::string ifpackType;
    Teuchos::ParameterList ifpackList;
    ifpackList.set("relaxation: sweeps", (LO) 5);
    ifpackList.set("relaxation: damping factor", (SC) 1.0);
    ifpackType = "RELAXATION";
    ifpackList.set("relaxation: type", "Symmetric Gauss-Seidel");
    RCP<SmootherPrototype> smoProto11     = rcp( new TrilinosSmoother(ifpackType, ifpackList, 0) );
    smoProto11->SetFactory("A", A11Fact);
    RCP<SmootherPrototype> smoProto22     = rcp( new TrilinosSmoother(ifpackType, ifpackList, 0) );
    smoProto22->SetFactory("A", A22Fact);

    //RCP<SmootherPrototype> smoProto11     = rcp( new DirectSolver("", Teuchos::ParameterList(), A11Fact) );
    //RCP<SmootherPrototype> smoProto22     = rcp( new DirectSolver("", Teuchos::ParameterList(), A22Fact) );

    RCP<SmootherFactory> Smoo11Fact = rcp( new SmootherFactory(smoProto11) );
    RCP<SmootherFactory> Smoo22Fact = rcp( new SmootherFactory(smoProto22) );

    RCP<FactoryManager> M11 = rcp(new FactoryManager());
    M11->SetFactory("A", A11Fact);
    M11->SetFactory("P", P11Fact);
    M11->SetFactory("Ptent", P11Fact); //for Nullspace
    M11->SetFactory("R", R11Fact);
    M11->SetFactory("Smoother", Smoo11Fact);
    M11->SetIgnoreUserData(true);

    RCP<FactoryManager> M22 = rcp(new FactoryManager());
    M22->SetFactory("A", A22Fact);
    M22->SetFactory("P", P22Fact);
    M22->SetFactory("R", R22Fact);
    M22->SetFactory("Ptent", P22TentFact); //for both P22 and Nullspace
    M22->SetFactory("Smoother", Smoo22Fact);
    M22->SetIgnoreUserData(true);

    RCP<BlockedPFactory> PFact = rcp(new BlockedPFactory());
    PFact->AddFactoryManager(M11);
    PFact->AddFactoryManager(M22);

    RCP<GenericRFactory> RFact = rcp(new GenericRFactory());

    RCP<Factory> AcFact = rcp(new BlockedRAPFactory());

    // Smoothers
    RCP<BlockedGaussSeidelSmoother> smootherPrototype     = rcp( new BlockedGaussSeidelSmoother() );
    smootherPrototype->SetParameter("Sweeps", Teuchos::ParameterEntry(2));
    smootherPrototype->SetParameter("Damping factor", Teuchos::ParameterEntry(1.0));
    smootherPrototype->AddFactoryManager(M11,0);
    smootherPrototype->AddFactoryManager(M22,1);
    RCP<SmootherFactory>   smootherFact          = rcp( new SmootherFactory(smootherPrototype) );

    // Coarse grid correction

//.........这里部分代码省略.........
    Op->getLocalRowView(0, indices, values);
    Op->replaceLocalValues(0, indices, newValues);
  }
  if (myRank == nCpus-1) { // JG TODO: can we use rowMap->isNodeLocalElement(lastRow) instead for more genericity?
    LO lastRow = Op->getNodeNumRows()-1;
    newValues[0] = -1.0; newValues[1] = 1.0;
    Op->getLocalRowView(lastRow, indices, values);
    Op->replaceLocalValues(lastRow, indices, newValues);
  }

  Op->fillComplete();
#endif // NEUMANN

  /**********************************************************************************/
  /*                                                                                */
  /**********************************************************************************/

  RCP<MultiVector> nullSpace = MultiVectorFactory::Build(map,1);
  nullSpace->putScalar( (SC) 1.0);
  Teuchos::Array<Teuchos::ScalarTraits<SC>::magnitudeType> norms(1);
  nullSpace->norm1(norms);
  if (comm->getRank() == 0)
    std::cout << "||NS|| = " << norms[0] << std::endl;

  RCP<MueLu::Hierarchy<SC,LO,GO,NO,LMO> > H = rcp( new Hierarchy() );
  H->SetDefaultVerbLevel(MueLu::Extreme);
  H->IsPreconditioner(false);

  RCP<MueLu::Level> Finest = H->GetLevel();
  Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
  Finest->Set("A", Op);
  Finest->Set("Nullspace", nullSpace);

  FactoryManager M;

  M.SetFactory("Aggregates", rcp(new CoupledAggregationFactory()));
  M.SetFactory("Ptent",      rcp(new TentativePFactory()));
  M.SetFactory("P",          rcp(new SaPFactory()));

#ifdef EMIN
  // Energy-minimization
  RCP<PatternFactory> PatternFact = rcp(new PatternFactory());
#if 0
  PatternFact->SetFactory("P", M.GetFactory("Ptent"));
#else
  PatternFact->SetFactory("P", M.GetFactory("P"));
#endif
  M.SetFactory("Ppattern",   PatternFact);

  RCP<EminPFactory> EminPFact = rcp(new EminPFactory());
  EminPFact->SetFactory("P", M.GetFactory("Ptent"));
  M.SetFactory("P",          EminPFact);

  RCP<NullspacePresmoothFactory> NullPreFact = rcp(new NullspacePresmoothFactory());
  NullPreFact->SetFactory("Nullspace", M.GetFactory("Nullspace"));
  M.SetFactory("Nullspace",  NullPreFact);
#endif

  RCP<SmootherPrototype> smooProto = gimmeMergedSmoother(nSmoothers, xpetraParameters.GetLib(), coarseSolver, comm->getRank());
  M.SetFactory("Smoother",   rcp(new SmootherFactory(smooProto)));

  Teuchos::ParameterList status;

  RCP<SmootherPrototype> coarseProto;
  if (maxLevels != 1)
    coarseProto = gimmeCoarseProto(xpetraParameters.GetLib(), coarseSolver, comm->getRank());

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

  {
    TimeMonitor tm(*TimeMonitor::getNewTimer("ScalingTest: 2 - MueLu Setup"));

    //
    // Hierarchy
    //

    // USER GUIDE     // instantiate new Hierarchy object
    H = rcp(new Hierarchy());
    H->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    H->SetMaxCoarseSize((GO) optMaxCoarseSize);
    // USER GUIDE     //

    //
    // Finest level
    //

    // USER GUIDE     // create a fine level object
    RCP<Level> Finest = H->GetLevel();
    Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    Finest->Set("A",           A);
    Finest->Set("Nullspace",   nullspace);
    Finest->Set("Coordinates", coordinates); //FIXME: XCoordinates, YCoordinates, ..
    // USER GUIDE //

    //
    // FactoryManager
    //

    // USER GUIDE     // define a factory manager
    FactoryManager M;
    // USER GUIDE     //

    //
    //
    // Aggregation
    //

    {
      RCP<UncoupledAggregationFactory> AggregationFact = rcp(new UncoupledAggregationFactory());
      *out << "========================= Aggregate option summary =========================" << std::endl;
      *out << "min DOFs per aggregate :                " << optMinPerAgg << std::endl;
      *out << "min # of root nbrs already aggregated : " << optMaxNbrSel << std::endl;
      AggregationFact->SetMinNodesPerAggregate(optMinPerAgg);  //TODO should increase if run anything othpermRFacter than 1D
      AggregationFact->SetMaxNeighAlreadySelected(optMaxNbrSel);
      std::transform(optAggOrdering.begin(), optAggOrdering.end(), optAggOrdering.begin(), ::tolower);
      if (optAggOrdering == "natural" || optAggOrdering == "random" || optAggOrdering == "graph") {
        *out << "aggregate ordering :                    " << optAggOrdering << std::endl;
        AggregationFact->SetOrdering(optAggOrdering);
      } else {
        std::string msg = "main: bad aggregation option """ + optAggOrdering + """.";
        throw(MueLu::Exceptions::RuntimeError(msg));
      }
      //AggregationFact->SetPhase3AggCreation(0.5);
      M.SetFactory("Aggregates", AggregationFact);

    *out << "=============================================================================" << std::endl;
    }

    //
    // Transfer
    //