Tpetra::CrsGraph creation, from a graph perspective

Tpetra::CrsGraph represents a bipartite graph, with edges leaving the set of rows and going to the set of columns. The graph is distributed over one or more MPI processes in a given MPI communicator. I use the words "rows" and "columns" because the graph will eventually represent the sparsity structure of a sparse matrix. However, you can just think of "rows" and "columns" as two possibly different sets of vertices.

Let’s say you want a graph with 5 rows, where the global indices of those rows are [0, 1, 3, 7, 10]. I deliberately made the set noncontiguous in this example; this is a perfectly valid graph or matrix in Tpetra. (I won’t say anything about the graph’s columns yet, but see below.) You want this graph to be distributed over exactly two processes, with rows [0, 1, 3] on Process 0, and rows [7, 10] on Process 1. Let’s first make a communicator that is guaranteed to have two processes:

// Teuchos::Comm: see Trilinos/packages/teuchos/comm/src/Teuchos_Comm.hpp
auto initialComm = Tpetra::getDefaultComm ();
TEUCHOS_TEST_FOR_EXCEPTION(initialComm.getSize () < 2, std::runtime_error, “This test requires at least two processes”);

// Select processes whose rank is less than 2.
const int color = (initialComm.getRank() < 2) ? 0 : 1;
const int key = 0; // see MPI_Comm_split documentation
auto comm = initialComm->split (color, key);

if (color == 0) {
  run_the_actual_test (comm); // exactly 2 processes
}

The rest of the code in this example is happening inside run_the_actual_test(comm), with comm as the communicator.

Now let’s make the row Map. If you look at the Doxygen documentation or at Trilinos/packages/tpetra/core/src/Tpetra_Map_decl.hpp, you’ll see that Map has three kinds of constructors. The third kind of constructor takes an arbitrary list of global indices on each process. That’s the constructor you want.

using map_type = Tpetra::Map<>;
using GO = map_type::global_ordinal_type;
using LO = map_type::local_ordinal_type;

const GO gblNumInds = 5;
const LO lclNumInds = (myRank == 0) ? 3 : 2;
std::vector<GO> myGIDs;
const int myRank = comm->getRank ();
if (myRank == 0) {
  myGIDs = {0, 1, 3};
}
else if (myRank == 1) {
  myGIDs = {7, 10};
}
// Make the row Map.
RCP<const map_type> rowMap (new map_type (gblNumInds, myGIDs.data (), lclNumInds, 0, comm));

Now we need to make the graph. You have many options, but they boil down to two:

1. Create graph with a row Map only
2. Create graph with a row Map and a column Map

CrsGraph is a bipartite graph, but it’s a distributed bipartite graph. If you create the graph with a column Map, you assert that for every edge on the calling process, the target vertex of that edge is in the column Map on the calling process. The column Map can have “extra” entries in it, but you pay for those in terms of communication. Thus, if you don’t know what should go in the column Map on each process, you should not specify the column Map when creating the CrsGraph. CrsGraph’s fillComplete method will create the column Map for you, based on the column indices you insert into the graph on each process.

Let’s suppose you want to insert the following (row, column) edges into the graph. Note how I arrange them by row index:

(0, 0), (0, 11)
(1, 3), (1, 4)
(3, 2)
(7, 5) (7, 7)
(10, 6)

Here’s how you could create the CrsGraph (this code works with StaticProfile):

Teuchos::ArrayRCP<size> numEntPerRow (lclNumInds); // indexed by local index
If (myRank == 0) {
  numEntPerRow[0] = 2; // (0, 0), (0, 11)
  numEntPerRow[1] = 2; // (1, 3), (1, 4)
  numEntPerRow[2] = 1; // (3, 2) – note that 3 is a global index here
}
else if (myRank == 1) {
  numEntPerRow[0] = 2; // (7, 5), (7, 7)
  numEntPerRow[1] = 1; // (10, 6) 
}

using crs_graph_type = Tpetra::CrsGraph<>;
RCP<crs_graph_type> G (new crs_graph_type (rowMap, numEntPerRow()));

And here’s how you would insert the edges:

std::vector<GO> myGblInds (2);
if (myRank == 0) {
  myGblInds[0] = 0; myGblInds[1] = 11;
  G->insertGlobalIndices (0, myGblInds.data (), 2); // (0, 0), (0, 11)
  myGblInds[0] = 3; myGblInds[1] = 4;
  G->insertGlobalIndices (1, myGblInds.data (), 2); // (1, 3), (1, 4)
  myGblInds[0] = 2;
  G->insertGlobalIndices (3, myGblInds.data (), 1); // (3, 2)
}
else if (myRank == 1) {
  myGblInds[0] = 5; myGblInds[1] = 7;
  G->insertGlobalIndices (7, myGblInds.data (), 2); // (7, 5), (7, 7)
  myGblInds[0] = 6;
  G->insertGlobalIndices (10, myGblInds.data (), 1); // (10, 6)
}

At this point, you are ready to call fillComplete on the graph. This does a few things:

1. Create a column Map (if the graph doesn’t already have one). The resulting column Map will have global indices [0, 3, 2, 4, 11] on Process 0, and [7, 5, 6] on Process 1. (Don’t worry about the very specific order of those column Map entries for now.)
2. Take in the domain and range Maps. Domain and range Maps interpret a graph as a matrix in the linear algebra sense. For example, the domain Map (NOT the column Map) defines the matrix’s number of columns.
3. Compute communication patterns for sparse matrix-vector multiply, using all four Maps.

For example, suppose I would like to interpret this graph as a non-square matrix with 5 rows (globally), and 12 columns (globally). I’ll create domain and range Maps as follows:

// rowMap works fine here as a range Map, because it has the right number of rows, and because it’s one to one.
RCP<const map_type> rangeMap = rowMap;

// I can’t use the column Map here, because it doesn’t have all the indices that actually would live in the range of the matrix.
// This is the Map constructor that takes nothing but a global number of indices, and an index base (the starting / globally smallest index).
const GO indexBase = 0; 
RCP<const map_type> domainMap (new map_type (12, indexBase, comm));

And now I’ll call fillComplete on the graph. This is a collective, in the MPI sense.

G->fillComplete (domainMap, rangeMap);