Declare location of include headers
In [1]:
.I /opt/install/do-conf-ep-serial/include/
Out[1]:
Load Trilinos libraries
In [2]:
.L libepetra
Out[2]:
In [3]:
.L libxpetra
Out[3]:
In [4]:
.L libmuelu
Out[4]:
In [5]:
.L libteuchoscore
Out[5]:
Include some standard header files to create a communicator and an Xpetra::Map
In [6]:
#include "Kokkos_DefaultNode.hpp"
#include "Epetra_SerialComm.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Teuchos_RCP.hpp"
#include "Teuchos_DefaultSerialComm.hpp"
#include "Xpetra_MapFactory.hpp"
#include "Xpetra_Map.hpp"
#include "Xpetra_EpetraCrsMatrix.hpp"
Out[6]:
Again, we use standard template parameters:
In [7]:
typedef double Scalar;
typedef int LocalOrdinal;
typedef int GlobalOrdinal;
typedef Kokkos::Compat::KokkosSerialWrapperNode EpetraNode;
typedef EpetraNode Node;
typedef Scalar SC;
typedef LocalOrdinal LO;
typedef GlobalOrdinal GO;
typedef Node NO;
Out[7]:
In [8]:
#include "Xpetra_CrsMatrixWrap.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_MultiVector.hpp"
#include "Xpetra_BlockedMultiVector.hpp"
#include "Xpetra_BlockedVector.hpp"
#include "Xpetra_MultiVectorFactory.hpp"
Out[8]:
Let's create a matrix and wrap it as an Xpetra::CrsMatrixWrap object
In [9]:
int NumMyElements = 30;
Out[9]:
In [10]:
Epetra_SerialComm Comm;
Epetra_Map Map(-1, NumMyElements, 0, Comm);
int NumGlobalElements = Map.NumGlobalElements();
Epetra_CrsMatrix A(Copy, Map, 3);
double negOne = -1.0;
double posTwo = 2.0;
double posOne = +1.0;
double zero = 0.0;
for (int i=0; i<NumMyElements; i++) {
int GlobalRow = A.GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1;
if (RowLess1!=-1) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);
if (RowPlus1!=NumGlobalElements) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);
A.InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);
};
int idxZero = 0;
int idxOne = 1;
A.ReplaceGlobalValues(0,1,&posOne, &idxZero);
A.ReplaceGlobalValues(0,1,&zero,&idxOne);
A.FillComplete();
Out[10]:
In [11]:
Teuchos::RCP<Epetra_CrsMatrix> rcpA = Teuchos::rcpFromRef(A);
Xpetra::EpetraCrsMatrixT<GO,NO> xA(rcpA);
Teuchos::RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > xCrsMat = Teuchos::rcpFromRef(xA);
Xpetra::CrsMatrixWrap<SC,LO,GO,NO> xCrsWrap = Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(xCrsMat);
Teuchos::RCP<Xpetra::Matrix<SC,LO,GO,NO> > xMat = Teuchos::rcpFromRef(xCrsWrap);
Out[11]:
The xCrsWrap object contains the matrix object that can be used as input for MueLu.
Smoother details are stored in a Teuchos::ParameterList. For demonstration purposes we create a Symmetric Gauss-Seidel smoother ($50$ sweeps, damping factor $0.45$):
In [12]:
#include "Teuchos_ParameterList.hpp"
Out[12]:
In [13]:
Teuchos::ParameterList smooParams;
Out[13]:
In [14]:
smooParams.set("type", Teuchos::ParameterEntry("RELAXATION"));
Teuchos::ParameterList& subSmooParams = smooParams.sublist("ParamterList");
subSmooParams.set("relaxation: type", "Symmetric Gauss-Seidel");
subSmooParams.set("relaxation: sweeps",50);
subSmooParams.set("relaxation: damping factor",0.45);
Out[14]:
You can easily print the content of a parameter list using
In [15]:
std::cout << smooParams << std::endl;
Out[15]:
The TrilinosSmoother class is a wrapper class for most relaxation-based and ILU type smoothers (provided by Ifpack or Ifpack2):
In [16]:
#include "MueLu_TrilinosSmoother.hpp"
Out[16]:
To create a new TrilinosSmoother you have to provide a type (e.g., RELAXATION) and the corresponding Ifpack/Ifpack2 parameter list:
In [17]:
MueLu::TrilinosSmoother<SC,LO,GO,NO> trSmoother("RELAXATION", subSmooParams);
Out[17]:
The smoother use variable A stored on the Level class as input. Therefore, let's create a Level container and store the matrix:
In [18]:
#include "MueLu_Level.hpp"
Out[18]:
In [19]:
MueLu::Level l;
Out[19]: