C++ FactoryManager::SetFactory方法代码示例

  TEUCHOS_UNIT_TEST(Hierarchy, KeepAggregates)
  {

    out << "version: " << MueLu::Version() << std::endl;

    RCP<const Teuchos::Comm<int> > comm = TestHelpers::Parameters::getDefaultComm();
    RCP<Matrix> A = TestHelpers::TestFactory<SC, LO, GO, NO>::Build1DPoisson(399*comm->getSize());

    Hierarchy H(A);
    H.SetMaxCoarseSize(1);

    RCP<CoupledAggregationFactory> CoupledAggFact = rcp(new CoupledAggregationFactory());
    FactoryManager M;
    M.SetFactory("Aggregates", CoupledAggFact);

    H.GetLevel(0)->Keep("Aggregates", CoupledAggFact.get());
    H.Setup(M, 0, 2);

    for (LocalOrdinal l=0; l<H.GetNumLevels()-1;l++) {
      TEST_EQUALITY(H.GetLevel(l)->IsAvailable("Aggregates", CoupledAggFact.get()), true);
    }

  } //FullPopulate_KeepAggregates

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

  ////////////////////////////////////////// prepare null space for A11
  RCP<MultiVector> nullspace11 = MultiVectorFactory::Build(xstridedvelmap, 2);  // this is a 2D standard null space

  for (int i=0; i<nDofsPerNode-1; ++i) {
    Teuchos::ArrayRCP<Scalar> nsValues = nullspace11->getDataNonConst(i);
    int numBlocks = nsValues.size() / (nDofsPerNode - 1);
    for (int j=0; j< numBlocks; ++j) {
      nsValues[j*(nDofsPerNode - 1) + i] = 1.0;
    }
  }

  Finest->Set("Nullspace1",nullspace11);

  ////////////////////////////////////////// prepare null space for A22
  RCP<MultiVector> nullspace22 = MultiVectorFactory::Build(xstridedpremap, 1);  // this is a 2D standard null space
  Teuchos::ArrayRCP<Scalar> nsValues22 = nullspace22->getDataNonConst(0);
  for (int j=0; j< nsValues22.size(); ++j) {
    nsValues22[j] = 1.0;
  }

  Finest->Set("Nullspace2",nullspace22);


  /////////////////////////////////////////// define rebalanced block AC factory
  // This is the main factory for "A" and defines the input for
  //   - the SubBlockAFactory objects
  //   - the rebalanced block Ac factory
  RCP<RebalanceBlockAcFactory> RebalancedAcFact = rcp(new RebalanceBlockAcFactory());

  /////////////////////////////////////////// define non-rebalanced blocked transfer ops
  RCP<BlockedPFactory> PFact = rcp(new BlockedPFactory()); // use row map index base from bOp
  RCP<GenericRFactory> RFact = rcp(new GenericRFactory());
  RFact->SetFactory("P", PFact);

  // non-rebalanced block coarse matrix factory
  // output is non-rebalanced coarse block matrix Ac
  // used as input for rebalanced block coarse factory RebalancedAcFact
  RCP<Factory> AcFact = rcp(new BlockedRAPFactory());
  AcFact->SetFactory("A", MueLu::NoFactory::getRCP());
  AcFact->SetFactory("P", PFact);  // use non-rebalanced block prolongator as input
  AcFact->SetFactory("R", RFact);  // use non-rebalanced block restrictor as input

  // define matrix sub-blocks of possibly rebalanced block matrix A
  // These are used as input for
  //   - the sub blocks of the transfer operators
  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));

  /////////////////////////////////////////// define rebalancing factories
  // define sub blocks of the coarse non-rebalanced block matrix Ac
  // input is the block operator generated by AcFact
  RCP<SubBlockAFactory> rebA11Fact = Teuchos::rcp(new SubBlockAFactory());
  rebA11Fact->SetFactory("A",AcFact);
  rebA11Fact->SetParameter("block row",Teuchos::ParameterEntry(0));
  rebA11Fact->SetParameter("block col",Teuchos::ParameterEntry(0));
  RCP<SubBlockAFactory> rebA22Fact = Teuchos::rcp(new SubBlockAFactory());
  rebA22Fact->SetFactory("A",AcFact);
  rebA22Fact->SetParameter("block row",Teuchos::ParameterEntry(1));
  rebA22Fact->SetParameter("block col",Teuchos::ParameterEntry(1));

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

    FactoryManager M;

    //
    //
    // 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") {
        *out << "aggregate ordering :                    NATURAL" << std::endl;
        AggregationFact->SetOrdering(MueLu::AggOptions::NATURAL);
      } else if (optAggOrdering == "random") {
        *out << "aggregate ordering :                    RANDOM" << std::endl;
        AggregationFact->SetOrdering(MueLu::AggOptions::RANDOM);
      } else if (optAggOrdering == "graph") {
        *out << "aggregate ordering :                    GRAPH" << std::endl;
        AggregationFact->SetOrdering(MueLu::AggOptions::GRAPH);
      } 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
    //

    {
      //
      // Non rebalanced factories
      //

      RCP<SaPFactory> PFact = rcp(new SaPFactory());
      PFact->SetDampingFactor(optSaDamping);

      RCP<Factory>    RFact = rcp(new TransPFactory());

      RCP<RAPFactory> AFact = rcp(new RAPFactory());
      AFact->setVerbLevel(Teuchos::VERB_HIGH);
      if (!optExplicitR) {
        H->SetImplicitTranspose(true);
        ParameterList Aclist = *(AFact->GetValidParameterList());
        Aclist.set("implicit transpose", true);
        AFact->SetParameterList(Aclist);
        if (comm->getRank() == 0) std::cout << "\n\n* ***** USING IMPLICIT RESTRICTION OPERATOR ***** *\n" << std::endl;
      }

      //
      // Repartitioning (if needed)
      //

  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,Node> > RP = Xpetra::MatrixMatrix<Scalar,LO,GO,Node>::Multiply(*R1,false,*P1,false,out);
    Teuchos::RCP<Xpetra::Matrix<Scalar,LO,GO,Node> > PtP = Xpetra::MatrixMatrix<Scalar,LO,GO,Node>::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());
//.........这里部分代码省略.........

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

    for (int i=0; i<numPDEs; ++i) {
      Teuchos::ArrayRCP<Scalar> nsValues = nullspace->getDataNonConst(i);
      int numBlocks = nsValues.size() / numPDEs;
      for (int j=0; j< numBlocks; ++j) {
        nsValues[j*numPDEs + i] = 1.0;
      }
    }

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

    fancyout << "========================================================\nGaleri complete.\n========================================================" << std::endl;

    // =========================================================================
    // Preconditioner construction
    // =========================================================================
    comm->barrier();
    tm = rcp(new TimeMonitor(*TimeMonitor::getNewTimer("ScalingTest: 1.5 - MueLu read XML")));

    RCP<Hierarchy> H = rcp ( new Hierarchy() );
    H->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    H->SetMaxCoarseSize(maxCoarseSize);

    // build finest Level
    RCP<MueLu::Level> Finest = H->GetLevel();
    Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    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)
    // =========================================================================

int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int argc, char *argv[]) {
#include <MueLu_UseShortNames.hpp>

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

  //
  // MPI initialization
  //

  Teuchos::oblackholestream blackhole;

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

    //
    // Process command line arguments
    //
    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)
    //

    // 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.

    H.Setup(M,startLevel,maxLevels);

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

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

//.........这里部分代码省略.........
    /////////////////////////////////////// generate MapExtractor object

    std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > xmaps;
    xmaps.push_back(xstridedvelmap);
    xmaps.push_back(xstridedpremap);

    Teuchos::RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > map_extractor = Xpetra::MapExtractorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(xstridedfullmap,xmaps);

    /////////////////////////////////////// build blocked transfer operator
    // using the map extractor
    Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bOp = Teuchos::rcp(new Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map_extractor,map_extractor,10));
    bOp->setMatrix(0,0,xA11);
    bOp->setMatrix(0,1,xA12);
    bOp->setMatrix(1,0,xA21);
    bOp->setMatrix(1,1,xA22);

    bOp->fillComplete();

    //////////////////////////////////////////////////// create Hierarchy
    RCP<Hierarchy> H = rcp ( new Hierarchy() );
    H->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    //H->setDefaultVerbLevel(Teuchos::VERB_NONE);
    H->SetMaxCoarseSize(maxCoarseSize);

    //////////////////////////////////////////////////////// finest Level
    RCP<MueLu::Level> Finest = H->GetLevel();
    Finest->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    Finest->Set("A",Teuchos::rcp_dynamic_cast<Matrix>(bOp));

    /////////////////////////////////////////////// define subblocks of A
    // make A11 block and A22 block available as variable "A" generated
    // by A11Fact and A22Fact
    RCP<SubBlockAFactory> A11Fact = 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 = rcp(new SubBlockAFactory());
    A22Fact->SetFactory("A",MueLu::NoFactory::getRCP());
    A22Fact->SetParameter("block row",Teuchos::ParameterEntry(1));
    A22Fact->SetParameter("block col",Teuchos::ParameterEntry(1));

    ////////////////////////////////////////// prepare null space for A11
    RCP<MultiVector> nullspace11 = MultiVectorFactory::Build(xstridedvelmap, 2);  // this is a 2D standard null space

    for (int i=0; i<nDofsPerNode-1; ++i) {
      Teuchos::ArrayRCP<Scalar> nsValues = nullspace11->getDataNonConst(i);
      int numBlocks = nsValues.size() / (nDofsPerNode - 1);
      for (int j=0; j< numBlocks; ++j) {
        nsValues[j*(nDofsPerNode - 1) + i] = 1.0;
      }
    }

    Finest->Set("Nullspace1",nullspace11);

    ///////////////////////////////////////// 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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

int main(int argc, char *argv[]) {
#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();

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

Teuchos::RCP<Vector> runExample(std::vector<size_t> stridingInfo, LocalOrdinal stridedBlockId, GlobalOrdinal offset) {
    using Teuchos::RCP;
    using Teuchos::rcp;

    RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
    RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    out->setOutputToRootOnly(0);
    *out << MueLu::MemUtils::PrintMemoryUsage() << std::endl;

    // Timing
    Teuchos::Time myTime("global");
    Teuchos::TimeMonitor Mt(myTime);

#ifndef HAVE_TEUCHOS_LONG_LONG_INT
    *out << "Warning: scaling test was not compiled with long long int support" << std::endl;
#endif

    // custom parameters
    LO maxLevels = 4;

    GO maxCoarseSize=1; //FIXME clp doesn't like long long int
    std::string aggOrdering = "natural";
    int minPerAgg=3;
    int maxNbrAlreadySelected=0;

    ////////////////////////////////////////////////////////////////////////////////////////
    // prepare redistribution of matrix (parallelization)
    int globalNumDofs = 7020;
    int nProcs = comm->getSize();
    int nDofsPerNode = 2;

    int nLocalDofs = (int) globalNumDofs / nProcs;
    nLocalDofs = nLocalDofs - (nLocalDofs % nDofsPerNode);
    int nCumulatedDofs = 0;
    sumAll(comm,nLocalDofs, nCumulatedDofs);

    if(comm->getRank() == nProcs-1) {
        nLocalDofs += globalNumDofs - nCumulatedDofs;
    }

    std::cout << "PROC: " << comm->getRank() << " numLocalDofs=" << nLocalDofs << std::endl;

    ////////////////////////////////////////////////////////////////////////////////////////
    // read in problem
    Epetra_Map emap (globalNumDofs, nLocalDofs, 0, *Xpetra::toEpetra(comm));
    Epetra_CrsMatrix * ptrA = 0;
    Epetra_Vector * ptrf = 0;
    Epetra_MultiVector* ptrNS = 0;

    std::cout << "Reading matrix market file" << std::endl;
    EpetraExt::MatrixMarketFileToCrsMatrix("stru2d_A.txt",emap,emap,emap,ptrA);
    EpetraExt::MatrixMarketFileToVector("stru2d_b.txt",emap,ptrf);
    EpetraExt::MatrixMarketFileToMultiVector( "stru2d_ns.txt", emap, ptrNS);
    RCP<Epetra_CrsMatrix> epA = Teuchos::rcp(ptrA);
    RCP<Epetra_Vector> epv = Teuchos::rcp(ptrf);
    RCP<Epetra_MultiVector> epNS = Teuchos::rcp(ptrNS);

    ////////////////////////////////////////////
    // Epetra_CrsMatrix -> Xpetra::Matrix
    RCP<CrsMatrix> exA = Teuchos::rcp(new Xpetra::EpetraCrsMatrix(epA));
    RCP<CrsMatrixWrap> crsOp = Teuchos::rcp(new CrsMatrixWrap(exA));
    RCP<Matrix> Op = Teuchos::rcp_dynamic_cast<Matrix>(crsOp);
    Op->SetFixedBlockSize(nDofsPerNode);

    // Epetra_Vector -> Xpetra::Vector
    RCP<Vector> xRhs = Teuchos::rcp(new Xpetra::EpetraVector(epv));

    RCP<MultiVector> xNS = Teuchos::rcp(new Xpetra::EpetraMultiVector(epNS));

    // Epetra_Map -> Xpetra::Map
    const RCP< const Map> map = Xpetra::toXpetra(emap);

    ////////////////////////////////////////////
    // create new MueLu hierarchy
    RCP<Hierarchy> H = rcp ( new Hierarchy() );
    H->setDefaultVerbLevel(Teuchos::VERB_HIGH);
    H->SetMaxCoarseSize(maxCoarseSize);

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

    // prepare CoalesceDropFactory
    RCP<CoalesceDropFactory> dropFact = rcp(new CoalesceDropFactory());
    //dropFact->SetVariableBlockSize();

    // prepare aggregation strategy
    RCP<CoupledAggregationFactory> CoupledAggFact = rcp(new CoupledAggregationFactory());
    CoupledAggFact->SetFactory("Graph", dropFact);
    *out << "========================= Aggregate option summary  =========================" << std::endl;
    *out << "min DOFs per aggregate :                " << minPerAgg << std::endl;
    *out << "min # of root nbrs already aggregated : " << maxNbrAlreadySelected << std::endl;
    CoupledAggFact->SetMinNodesPerAggregate(minPerAgg); //TODO should increase if run anything other than 1D
    CoupledAggFact->SetMaxNeighAlreadySelected(maxNbrAlreadySelected);
    std::transform(aggOrdering.begin(), aggOrdering.end(), aggOrdering.begin(), ::tolower);
    if (aggOrdering == "natural") {
        *out << "aggregate ordering :                    NATURAL" << std::endl;
        CoupledAggFact->SetOrdering(MueLu::AggOptions::NATURAL);
//.........这里部分代码省略.........

  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

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

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

  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();
    RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
    out->setOutputToRootOnly(0);
    *out << MueLu::MemUtils::PrintMemoryUsage() << std::endl;

    //#ifndef HAVE_XPETRA_INT_LONG_LONG
    *out << "Warning: scaling test was not compiled with long long int support" << std::endl;
    //#endif

    /**********************************************************************************/
    /* SET TEST PARAMETERS                                                            */
    /**********************************************************************************/
    // Note: use --help to list available options.
    Teuchos::CommandLineProcessor clp(false);

    Xpetra::Parameters xpetraParameters(clp);             // manage parameters of xpetra

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

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

    xpetraParameters.check();

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

    if (comm->getRank() == 0) {
      std::cout << xpetraParameters;
      // TODO: print custom parameters // Or use paramList::print()!
    }

    /**********************************************************************************/
    /* CREATE INITIAL MATRIX                                                          */
    /**********************************************************************************/
    RCP<const Map> bigMap;
    RCP<const Map> map1;
    RCP<const Map> map2;
    GO numElements = 500;
    GO numElements1 = 400;
    GO numElements2 = 100;

    //bigMap = MapFactory::Build(Xpetra::UseEpetra, numElements,  0, comm); // ok this is the problem :-)
    std::vector<size_t> stridingInfo;
    stridingInfo.push_back(1);
    map1   = StridedMapFactory::Build(lib, numElements1, 0, stridingInfo, comm, -1);
    map2   = StridedMapFactory::Build(lib, numElements2, numElements1, stridingInfo, comm, -1);

    std::vector<GlobalOrdinal> localGids; // vector with all local GIDs on cur proc
    Teuchos::ArrayView< const GlobalOrdinal > map1eleList = map1->getNodeElementList(); // append all local gids from map1 and map2
    localGids.insert(localGids.end(), map1eleList.begin(), map1eleList.end());
    Teuchos::ArrayView< const GlobalOrdinal > map2eleList = map2->getNodeElementList();
    localGids.insert(localGids.end(), map2eleList.begin(), map2eleList.end());
    Teuchos::ArrayView<GlobalOrdinal> eleList(&localGids[0],localGids.size());
    bigMap = StridedMapFactory::Build(lib, numElements, eleList, 0, stridingInfo, comm); // create full big map (concatenation of map1 and map2)

    std::vector<Teuchos::RCP<const Map> > maps;
    maps.push_back(map1); maps.push_back(map2);

    Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LO, GO, Node> > mapExtractor = Xpetra::MapExtractorFactory<Scalar,LO,GO,Node>::Build(bigMap, maps);

    RCP<CrsMatrixWrap> Op11 = MueLuTests::GenerateProblemMatrix(map1,2,-1,-1);
    RCP<CrsMatrixWrap> Op22 = MueLuTests::GenerateProblemMatrix(map2,3,-2,-1);

    /*Op11->describe(*out,Teuchos::VERB_EXTREME);
      Op22->describe(*out,Teuchos::VERB_EXTREME);*/

    // build blocked operator
    Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar,LO,GO,Node> > bOp = Teuchos::rcp(new Xpetra::BlockedCrsMatrix<Scalar,LO,GO>(mapExtractor,mapExtractor,10));

    Teuchos::RCP<Xpetra::CrsMatrix<Scalar,LO,GO,Node> > crsMat11 = Op11->getCrsMatrix();
    Teuchos::RCP<Xpetra::CrsMatrix<Scalar,LO,GO,Node> > crsMat22 = Op22->getCrsMatrix();
    bOp->setMatrix(0,0,crsMat11);
    bOp->setMatrix(1,1,crsMat22);

    bOp->fillComplete();

    // build hierarchy
    Hierarchy H;
    H.SetMaxCoarseSize(50);
    RCP<Level> levelOne = H.GetLevel();
    levelOne->Set("A", Teuchos::rcp_dynamic_cast<Matrix>(bOp)); // set blocked operator

    RCP<SubBlockAFactory> A11Fact = Teuchos::rcp(new SubBlockAFactory());
    A11Fact->SetFactory("A",MueLu::NoFactory::getRCP());
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  }
  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());
  else
    coarseProto = gimmeMergedSmoother(nSmoothers, xpetraParameters.GetLib(), coarseSolver, comm->getRank());