Using MueLu in User Applications

This tutorial demonstrates how to use MueLu from within user applications in C++. We will use the Tpetra linear algebra stack. We will read the configuration for MueLu from an XML file and then create a MueLu::Hierarchy. This will then be used as a preconditioner within Belos as well as a standalone solver.

Note

There is also support for Stratimikos. Please refer to the examples in the MueLu folder for more details.

XML Interface using CreateTpetraPreconditioner

The most comfortable way to declare the multigrid parameters for MueLu is using the XML interface. In fact, MueLu provides two different XML interfaces. There is a simplified XML interface for multigrid users and a more advanced XML interface for experts, which allows to make use of all features of MueLu as a multigrid framework. Both XML file formats are introduced in the previous sections of this hands on tutorial. However, for the C++ code it makes no difference which type of XML interface is used.

Note

In the next sections, we give some code snippets. THey are borrowed form the TutorialDriver.cpp file in the tutorial tests.

Preparations

First of all, we need to grab a communicator object. Therefore, it is easy to use some utilities from the Teuchos package:

      RCP<const Teuchos::Comm<int>> comm = Teuchos::DefaultComm<int>::getComm();
      [[maybe_unused]] int MyPID         = comm->getRank();
      [[maybe_unused]] int NumProc       = comm->getSize();

For the multigrid method, we need a linear operator \(A\). We can generate a problem-dependent matrix operator using the Galeri package (see Example problem), with Tpetra as our underlying linear algebra framework (through Xpetra wrappers). Below is an example for the Laplace2D Problem:

      // Create the matrix
      RCP<Galeri::Xpetra::Problem<Map, CrsMatrixWrap, MultiVector>> galeriProblem =
          Galeri::Xpetra::BuildProblem<SC, LO, GO, Map, CrsMatrixWrap, MultiVector>("Laplace2D", xpetra_map, galeriList);
      RCP<Matrix> A = galeriProblem->BuildMatrix();

For aggregation-based algebraic multigrid methods, one has to provide a valid set of near null space vectors to produce transfer operators. In case of a Laplace problem, we just use a constant vector.

      // Build null space vector for a scalar problem, i.e. vector with all ones
      RCP<MultiVector> nullspace = rcp(new MultiVector(dofMap, 1));
      nullspace->putScalar(one);

Setup phase

With a fine level operator \(A\) available as Tpetra::CrsMatrix object and a set of near null space vectors (available as Tpetra::MultiVector), all minimum requirements are fulfilled for generating an algebraic multigrid hierarchy. There are two different ways to setup a multigrid hierarchy in MueLu. One can either use a parameter list driven setup process which accepts either Teuchos::ParameterList objects or XML files in two different XML file formats. Alternatively, one can use the MueLu C++ API to define the multigrid setup at compile time. In the next sections we show both variants.

The most comfortable way to declare the multigrid parameters for MueLu is using the XML interface. In fact, MueLu provides two different XML interfaces. There is a simplified XML interface for multigrid users and a more advanced XML interface for expert which allows to make use of all features of MueLu as a multigrid framework. Both XML file formats are introduced in the previous sections of this hands on tutorial. However, for the C++ code it makes no difference which type of XML interface is used.

We first read the MueLu configuration from the XML file:

      // Read MueLu parameter list from xml file
      RCP<ParameterList> mueluParams = Teuchos::getParametersFromXmlFile(xmlFileName);

Then, we store the near null space vectors in the "user data" sublist of this parameter list:

      // Register nullspace as user data in the MueLu parameter list
      ParameterList& userDataList = mueluParams->sublist("user data");
      userDataList.set<RCP<MultiVector>>("Nullspace", nullspace);

We then can create a MueLu object ready to be used as a preconditioner:

      // Create the MueLu preconditioner based on the Tpetra stack
      RCP<MueLu::TpetraOperator<SC, LO, GO, NO>> mueLuPreconditioner =
          MueLu::CreateTpetraPreconditioner(Teuchos::rcp_dynamic_cast<Operator>(matrix), *mueluParams);

Iteration Phase

Once the setup phase is completed, the MueLu multigrid hierarchy is ready for being used.

There are several ways how to use the multigrid method. One can apply the multigrid method as standalone solver for linear systems. Multigrid methods are also known to be efficient preconditioners within iterative (Krylov) solvers such as CG or GMRES methods.

In the next subsections it is demonstrated how to use MueLu as standalone solver and as preconditioner for iterative solvers from the Belos and AztecOO package in Trilinos.

For solving a linear system \(Ax=b\), we need a right hand side vector \(b\). When using iterative solvers we also need an initial guess for the solution vector.

In this example we just create Tpetra::MultiVectors (with just a single vector each). The right-hand side vector \(b\) is initialized with ones and the solution vector \(x\) is filled with random values.

MueLu as multigrid solver

MueLu can be used as standalone solver. First, we create and initialize a solution vector:

      // Create solution vector and set it to an initial guess (X0)
      RCP<MultiVector> multigridSolVec = rcp(new MultiVector(dofMap, 1, false));
      Tpetra::deep_copy(*multigridSolVec, *X0);

If necessary, extract the MueLu::Hiearchy from the Tpetra preconditioner object:

        // Extract the underlying MueLu hierarchy
        RCP<MueLu::Hierarchy<SC, LO, GO, NO>> hierarchy = mueLuPreconditioner->GetHierarchy();

Then, the MueLu::Hierarchy object is set to the non-preconditioner mode:

        // Configure MueLu to be used as solver
        hierarchy->IsPreconditioner(false);

Finally, we solve the system by calling the Iterate() routine to perform mgridSweeps sweeps with the chosen multigrid cycle.

        // Solve
        hierarchy->Iterate(*Xpetra::toXpetra(B), *Xpetra::toXpetra(multigridSolVec), mgridSweeps);

If successful, the multigridSolVec vector contains the solution.

MueLu as preconditioner for Belos

Belos is the Krylov package in Trilinos and works both for Epetra and Tpetra. Here, we demonstrate how to use MueLu as preconditioner for Belos solvers using Tpetra.

First, we create and initialize a solution vector:

      // Create solution vector and set it to an initial guess (X0)
      RCP<MultiVector> precSolVec = rcp(new MultiVector(dofMap, 1, false));
      Tpetra::deep_copy(*precSolVec, *X0);

Then, we create a Belos::LinearProblem and hand in the MueLu preconditioner object:

        // Construct a Belos LinearProblem object and hand-in the MueLu preconditioner
        RCP<Belos::LinearProblem<SC, MultiVector, Operator>> belosProblem =
            rcp(new Belos::LinearProblem<SC, MultiVector, Operator>(matrix, precSolVec, B));
        belosProblem->setLeftPrec(mueLuPreconditioner);
        bool set = belosProblem->setProblem();

Now, we define the linear solver configuration in a Teuchos::ParameterList and create the Belos solver with this configuration:

        // Belos parameter list
        RCP<ParameterList> belosList = Teuchos::parameterList();
        belosList->set("Maximum Iterations", 500);     // Maximum number of iterations allowed
        belosList->set("Convergence Tolerance", tol);  // Relative convergence tolerance requested
        belosList->set("Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails);
        belosList->set("Output Frequency", 1);
        belosList->set("Output Style", Belos::Brief);

        // Create an iterative solver manager
        Belos::SolverFactory<SC, MultiVector, Operator> solverFactory;
        RCP<Belos::SolverManager<SC, MultiVector, Operator>> solver;
        if (dsolveType == "gmres")
          solver = solverFactory.create("Block GMRES", belosList);
        else if (dsolveType == "cg")
          solver = solverFactory.create("CG", belosList);
        else {
          fancyout << "\nERROR:  Please select a valid solver!" << std::endl;
          return EXIT_FAILURE;
        }

        // Pass the linear problem to the solver
        solver->setProblem(belosProblem);

Finally, we solve the linear system:

        // Perform solve
        Belos::ReturnType retStatus = Belos::Unconverged;
        retStatus                   = solver->solve();

Full example using the XML interface

The reader may refer to TutorialDriver.cpp for a working example to study the source code. This demonstration program has some more features that are not discussed in this tutorial.

Exercise 1

Compile the example in TutorialDriver.cpp and then run the program in parallel using two processors

mpirun -np 2 ./MueLu_TutorialDriver.exe --help

Study the screen output and try to run the example with an XML file as input for the multigrid setup.

Exercise 2

Create large scale examples using the --nx and --ny parameters for a finer mesh.

Choose reasonable numbers for --nx and --ny for your machine and make use of your knowledge about MueLu for generating efficient preconditioners.

C++ Interface

As an alternative to the XML interfaces, the user can also define the multigrid hierarchy using the C++ API directly. In contrary to the XML interface, which allows to build the layout of the multigrid preconditioner at runtime, the preconditioner is fully defined at compile time when using the C++ interface.

First, a MueLu::Hierarchy object has to be defined, which manages the multigrid hierarchy including all multigrid levels. It provides routines for the multigrid setup and the multigrid cycle algorithms (such as V-cycle and W-cycle).

      // create new hierarchy
      RCP<MueLu::Hierarchy<SC, LO, GO, NO>> H;

There are some member functions which can be used to describe the basic multigrid hierarchy. The SetMaxCoarseSize member function is used to set the maximum size of the coarse level problem before the coarsening process can be stopped.

      // Instantiate new Hierarchy object
      H = rcp(new Hierarchy());
      H->setDefaultVerbLevel(Teuchos::VERB_HIGH);
      H->SetMaxCoarseSize((GO)50);

Next, one defines an empty MueLu::Level object for the finest level. The MueLu::Level objects represent a data container storing the internal variables on each multigrid level. The user has to provide and fill the level container for the finest level only. The MueLu::Hierarchy object then automatically generates the coarse levels using the multigrid parameters. The absolute minimum requirements for the finest level that the user has to provide is the fine level operator \(A\) which represents the fine level matrix. MueLu is based on Xpetra. So, the matrix \(A\) has to be of type Xpetra::Matrix. In addition, the user should also provide a valid set of near null space vectors. For a Laplace problem we can just use the constant nullspace vector that has previously been defined. Some routines need additional information. For example, the user has to provide the node coordinates for repartitioning.

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

Note

When including the MueLu_UseShortNames.hpp header file,

the template parameters usually can be dropped for compiling. The most important template parameters are SC for the scalar type, LO for the local ordinal type and GO for the global ordinal type. For a detailed description of the template parameters, the reader may refer to the Tpetra documentation.

A MueLu::FactoryManager object is used for the internal management of data dependencies and generating algorithms of the multigrid setup. Even though not absolutely necessary, we show the usage of the MueLu::FactoryManager object as it allows for user-specific enhancements of the multigrid code.

      MueLu::FactoryManager M;
      M.SetKokkosRefactor(false);

The user can define its own factories for performing different tasks in the setup process. The following code shows how to define a smoothed aggregation transfer operator and a restriction operator. The MueLu::RAPFactory is used for the (standard) Galerkin product to generate the coarse level matrix \(A\).

      // declare some factories (potentially overwrite default factories)
      RCP<SaPFactory> PFact = rcp(new SaPFactory());
      PFact->SetParameter("sa: damping factor", Teuchos::ParameterEntry(4. / 3));

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

      RCP<RAPFactory> AcFact = rcp(new RAPFactory());
      AcFact->setVerbLevel(Teuchos::VERB_HIGH);

      H->SetImplicitTranspose(true);
      Teuchos::ParameterList Aclist = *(AcFact->GetValidParameterList());
      Aclist.set("transpose: use implicit", true);
      AcFact->SetParameterList(Aclist);

The user-defined factories have to be registered in the FactoryManager using the lines

      // configure factory manager (no repartitioning)
      M.SetFactory("P", PFact);
      M.SetFactory("R", RFact);
      M.SetFactory("A", AcFact);

Warning

If you forget to register the new factories, the FactoryManager will use some internal default factories for being responsible to create the corresponding variables.

Then your user-specified factories are just ignored during the multigrid setup!

Note

The FactoryManager is also responsible for resolving all dependencies between different factories.

That is, after the user-defined factories have been registered, all factories that request variable \(P\) are provided with the prolongation operator \(P\) that has been generated by the registered factory PFact. If there is some data requested for which no factory has been registered by the user, the FactoryManager manages an internal list for reasonable default choices and default factories.

Next, the user has to declare a level smoother. The following code can be used to define a symmetric Gauss-Seidel smoother. Other methods can be set up in a similar way.

      // define smoother object
      std::string ifpackType;
      Teuchos::ParameterList ifpackList;
      ifpackList.set("relaxation: sweeps", (LO)2);
      ifpackList.set("relaxation: damping factor", (SC)1.0);
      ifpackType = "RELAXATION";
      ifpackList.set("relaxation: type", "Symmetric Gauss-Seidel");

Before the level smoother can be used, a MueLu::SmootherFactory has to be defined for the smoother factory. The SmootherFactory is used in the multigrid setup to generate level smoothers for the corresponding levels using the prototyping design pattern. Note, that the SmootherFactory has also to be registered in the FactoryManager object. If the user forgets this, the multigrid setup will use some kind of default smoother, i.e., the user-chosen smoother options are just ignored.

      // create smoother factory
      RCP<SmootherPrototype> smootherPrototype = rcp(new TrilinosSmoother(ifpackType, ifpackList));
      M.SetFactory("Smoother", rcp(new SmootherFactory(smootherPrototype)));

Once the FactoryManager is set up, it can be used with the Hierarchy::Setup routine to initiate the coarsening process and set up the multigrid hierarchy.

      // setup multigrid hierarchy
      int startLevel = 0;
      H->Setup(M, startLevel, 10);

Footnotes

[1]

12. Berger-Vergiat, C. A. Glusa, G. Harper, J. J. Hu, M. Mayr, P. Ohm, A. Prokopenko, C. M. Siefert, R. S. Tuminaro, and T. A. Wiesner. MueLu User’s Guide. Technical Report SAND2023-12265, Sandia National Laboratories, Albuquerque, NM (USA) 87185, 2023.