Map, MultiVector, CrsMatrix,...)Factory class to create new instances of the linear algebra objects without writing {E,T}petra specific codeBlockedMap, BlockedMultiVector,...)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]:
Create a serial communicator...
In [3]:
#include "Epetra_SerialComm.h"
Out[3]:
In [4]:
Epetra_SerialComm Comm;
Out[4]:
... and define the matrix size.
In [5]:
int NumMyElements = 20;
Out[5]:
Create a new map and fill the matrix.
In [6]:
#include "Epetra_Map.h"
#include "Epetra_CrsMatrix.h"
Out[6]:
In [7]:
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;
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);
};
A.FillComplete();
Out[7]:
The following code demonstrates how to wrap the Epetra_CrsMatrix object into an Xpetra::Matrix object that serves as abstract input class for most of the algebraic multigrid algorithms in MueLu.
Even though the example uses Epetra, the procedure for Tpetra is similar.
In [8]:
.L libxpetra
Out[8]:
First, we have to include some header files:
In [9]:
#include <Kokkos_DefaultNode.hpp>
Out[9]:
Next, we declare some typedefs:
In [10]:
typedef double Scalar;
typedef int LocalOrdinal;
typedef int GlobalOrdinal;
typedef Kokkos::Compat::KokkosSerialWrapperNode EpetraNode;
typedef EpetraNode Node;
Out[10]:
If you know Tpetra you are certainly familiar with most of the template parameters. Xpetra is templated on the exact same template parameters. Therefore, we have to provide compatible template parameters for the Epetra branch in Xpetra.
Please note the EpetraNode typedef which can be either the SerialNode or OpenMPNode depending on the configuration of Epetra on your machine. This notebook only contains a standard Epetra installation using GO=int on the SerialNode.
To improve the readability of the source code, you should include the Xpetra_UseShortNames.hpp file after the typedefs. It contains typedefs of all Xpetra classes with the corresponding template parameters defined above. This allows us to write, e.g., Map instead of Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> later in the code. However, this only works if your code lives in the Xpetra namespace. This is not the case in this notebook. So, we have to write the long version. At least, we can shorten the template parameters using SC, LO, GO, NO. The corresponding typedefs come with the inclusion of Xpetra_UseShortNames.hpp
In [11]:
#include "Xpetra_UseShortNames.hpp"
Out[11]:
Xpetra_EpetraCrsMatrix.hpp defines the class which encapsulates an Epetra_CrsMatrix object in Xpetra.
In [12]:
#include "Xpetra_EpetraCrsMatrix.hpp"
Out[12]:
The constructor of Xpetra::EpetraCrsMatrix only accepts an RCP pointer to an Epetra_CrsMatrix object. So, let's wrap A into an rcp pointer:
In [13]:
#include "Teuchos_RCP.hpp"
Out[13]:
In [14]:
Teuchos::RCP<Epetra_CrsMatrix> rcpA = Teuchos::rcpFromRef(A);
std::cout << rcpA << std::endl;
Out[14]:
Then call the constructor
In [15]:
Xpetra::EpetraCrsMatrixT<GO,NO> xA(rcpA);
Out[15]:
In [16]:
std::cout << *rcpA << std::endl;
Out[16]:
The {E,T}petraCrsMatrix class are package-specific wrappers for Epetra and Tpetra which encapsulate an Epetra Epetra_CrsMatrix or an Tpetra Tpetra::CrsMatrix object.
The EpetraCrsMatrix object is templated on the GlobalOrdinal and NodeType. GlobalOrdinal can be either int (default) or long long depending on the Epetra configuration. In this container we only support GO=int. The NodeType template parameter should be Xpetra::EpetraNode which translates to the serial or OpenMP node, depending on the configuration of Epetra. Note, that Epetra can only be compiled for one choice of GlobalOrdinal and NodeType.
The TpetraCrsMatrix object is templated on the Tpetra template parameters, which include Scalar, LocalOrdinal, GlobalOrdinal and Node.
Xpetra consequently follows the Tpetra style and naming conventions for functions.
Next, one has to wrap the Xpetra::CrsMatrix object into an Xpetra::Matrix object. This can be done with the Xpetra::CrsMatrixWrap class.
In [17]:
#include "Xpetra_CrsMatrixWrap.hpp"
#include "Xpetra_CrsMatrix.hpp"
Out[17]:
Again, we have to use RCP pointers from the Teuchos package.
In [18]:
Teuchos::RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > xCrsMat = Teuchos::rcpFromRef(xA);
Out[18]:
In [19]:
Teuchos::RCP<Xpetra::Matrix<SC,LO,GO,NO> > xMat = Teuchos::rcp(new Xpetra::CrsMatrixWrap<SC,LO,GO,NO>(xCrsMat));
GO numRows = xMat->getGlobalNumRows();
auto numEntries = xMat->getNodeNumEntries();
std::cout << "The matrix has " << numRows << " rows and " << numEntries << " nonzeros." << std::endl;
Out[19]:
The object xMat is an RCP pointer to a Xpetra::Matrix object that can be used as (abstract) input to MueLu.
It supports the Tpetra style interface, i.e., you can basically use it like an Tpetra object.
First of all: you never should need that!
Xpetra is meant to provide a common wrapper for {E,T}petra which allows you to write algorithms independent from the concrete implementation of {E,T}petra, but providing a Petra like API! This is a big advantage compared to, e.g., Thyra, which has a similar intention, but is too abstract for writing algorithms which operate on the low-level data (like matrix entries).
Nevertheless, here is how you can access the underlying Epetra_CrsMatrix object given an Xpetra::Matrix object:
In [20]:
Teuchos::RCP<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> > xCrsMatWrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO> >(xMat);
Teuchos::RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > xCrsMatTmp = xCrsMatWrap->getCrsMatrix();
Teuchos::RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > xEpCrsMatTmp = Teuchos::rcp_dynamic_cast<Xpetra::EpetraCrsMatrixT<GO,NO>>(xCrsMatTmp);
Teuchos::RCP<const Epetra_CrsMatrix> epMat = xEpCrsMatTmp->getEpetra_CrsMatrix();
std::cout << "The Epetra matrix has " << epMat->NumGlobalRows()<< " rows." << std::endl;
Out[20]:
Again: you should not need this. As you can see from the last code snippet, it introduces {E,T}petra specific code which is against the Xpetra philosophy. You use Xpetra to write code independent from {E,T}petra!
Again, let's assume we have an Epetra_Vector:
In [21]:
Epetra_Vector x(Map);
x.Random();
Out[21]:
We can wrap it into a Xpetra::EpetraVectorT class using
In [22]:
#include "Xpetra_EpetraVector.hpp"
Out[22]:
In [23]:
Xpetra::EpetraVectorT<GO,NO> xx(Teuchos::rcpFromRef(x));
Out[23]:
Let's print the content of the vector. It contains the random values of the underlying Epetra vector.
In [24]:
std::cout << xx << std::endl;
Out[24]:
The Xpetra::EpetraVectorT class is derived from the Xpetra::Vector class:
In [25]:
Teuchos::RCP<Xpetra::Vector<SC,LO,GO,NO> > xv = Teuchos::rcpFromRef(xx);
std::cout << *xv << std::endl;
Out[25]:
To wrap multivectors for {E,T}petra requires a similar procedure.
If you write an algorithm in Xpetra (independent from the concrete {E,T}petra implementation), you may need to create objects such as (multi)vectors or matrices. You cannot use the {E,T}petra wrappers given by, e.g., Xpetra::EpetraMultiVector or Xpetra::TpetraCrsMatrix, since this would make your code {E,T}petra specific.
Instead, you use the Xpetra factory classes to create new instances of maps, vectors, multivectors, matrices, etc...
In [26]:
#include "Xpetra_MapFactory.hpp"
#include "Xpetra_Map.hpp"
#include "Teuchos_DefaultSerialComm.hpp"
Out[26]:
Xpetra uses the communicator object from Teuchos:
In [27]:
Teuchos::RCP<const Teuchos::Comm<int> > comm = Xpetra::toXpetra(Comm);
Out[27]:
The factory classes have Build routines which accept a flag Xpetra::Use{E,T}petra and user parameters for the constructor of the Map object. The code, e.g., could look like:
In [28]:
Xpetra::UnderlyingLib lib = Xpetra::UseEpetra;
Out[28]:
In [29]:
Teuchos::RCP<Xpetra::Map<LO,GO,NO> > xm = Xpetra::MapFactory<LO,GO,NO>::Build(lib,10,0,comm);
Out[29]:
By consequently using lib when creating map objects through Xpetra factories one can quickly switch from an Epetra to a Tpetra implementation by setting lib = Xpetra::UseTpetra. No further code changes necessary.
In [30]:
std::cout << *xm <<std::endl;
Out[30]:
There is a Xpetra::MultiVectorFactory which works the same way:
In [31]:
#include "Xpetra_MultiVector.hpp"
#include "Xpetra_Vector.hpp"
#include "Xpetra_EpetraMultiVector.hpp"
#include "Xpetra_BlockedMultiVector.hpp"
#include "Xpetra_BlockedVector.hpp"
Out[31]:
In [32]:
#include "Xpetra_MultiVectorFactory.hpp"
Out[32]:
The underlying implementation (Epetra or Tpetra) is already specified by the type of the map.
In [33]:
Teuchos::RCP<Xpetra::MultiVector<SC,LO,GO,NO> > xmv = Xpetra::MultiVectorFactory<SC,LO,GO,NO>::Build(xm, 3);
xmv->putScalar(42.0);
std::cout << *xmv << std::endl;
Out[33]:
In [ ]: