EpetraLesson05

Lesson topics

The Lesson 02 Vector example introduces and describes the Epetra_Map class, which is Epetra's representation of a data distribution. This example builds on that by showing how to use Maps and the Epetra_Export class to redistribute data. In this case, we build a sparse matrix on one MPI process, and redistribute it to a sparse matrix stored in block row fashion, with an equal number of rows per process.

Epetra data redistribution: Import and Export

The Epetra_Map class describes a data distribution over one or more distributed-memory parallel processes. It "maps" global indices (unique labels for the elements of a data structure) to parallel processes. This ability to describe a data distribution calls for a redistribution capability, that is, to reorganize or remap data from one distribution to another. Epetra provides this capability through the Epetra_Import and Epetra_Export classes.

Import redistributes from a uniquely owned (one-to-one) Map to a possibly not uniquely owned Map. Export redistributes from a possibly not uniquely owned to a uniquely owned Map. We distinguish between these cases both for historical reasons and for performance reasons.

Import and Export objects encapsulate and remember a communication pattern for reuse. Computing the computation pattern requires nontrivial work, but keeping around the Import or Export object lets you reuse that work. This is very important for operations that are performed frequently, such as the Import and Export operations in Epetra's sparse matrix-vector multiply.

In both cases, Import and Export let the user specify how to combine incoming new data with existing data that has the same global index. For example, one may replace old data with new data or sum them together.

Code example

This example shows how to migrate the data in Epetra objects (sparse matrices and vectors) between two different parallel distributions.

// This defines useful macros like HAVE_MPI, which is defined if and
// only if Epetra was built with MPI enabled.
#include <Epetra_config.h>

#ifdef HAVE_MPI
#  include <mpi.h>
// Epetra's wrapper for MPI_Comm.  This header file only exists if
// Epetra was built with MPI enabled.
#  include <Epetra_MpiComm.h>
#else
#  include <Epetra_SerialComm.h>
#endif // HAVE_MPI

#include <Epetra_CrsMatrix.h>
#include <Epetra_Export.h>
#include <Epetra_Map.h>
#include <Epetra_Vector.h>
#include <Epetra_Version.h>

#include <sstream>
#include <stdexcept>


// The type of global indices.  You could just set this to int,
// but we want the example to work for Epetra64 as well.
#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
// Epetra was compiled only with 64-bit global index support,
// so use 64-bit global indices.
#    define EXAMPLE_USES_64BIT_GLOBAL_INDICES 1
typedef long long global_ordinal_type;
#else
#  ifdef EPETRA_NO_64BIT_GLOBAL_INDICES
// Epetra was compiled only with 32-bit global index support,
// so use 32-bit global indices.
typedef int global_ordinal_type;
#  else
#    define EXAMPLE_USES_64BIT_GLOBAL_INDICES 1
// Epetra was compiled with both 64-bit and 32-bit global index
// support.  Use 64-bit global indices for maximum generality.
typedef long long global_ordinal_type;
#  endif // EPETRA_NO_64BIT_GLOBAL_INDICES
#endif // EPETRA_NO_32BIT_GLOBAL_INDICES


// Create and return a pointer to an example CrsMatrix, with row
// distribution over the given Map.  The caller is responsible for
// freeing the result.
Epetra_CrsMatrix*
createCrsMatrix (const Epetra_Map& map)
{
  const Epetra_Comm& comm = map.Comm ();

  // Create an Epetra_CrsMatrix using the Map, with dynamic allocation.
  Epetra_CrsMatrix* A = new Epetra_CrsMatrix (Copy, map, 3);

  // The list of global indices owned by this MPI process.

  const global_ordinal_type* myGblElts = NULL;
  global_ordinal_type numGblElts = 0;
#ifdef EXAMPLE_USES_64BIT_GLOBAL_INDICES
  myGblElts = map.MyGlobalElements64 ();
  numGblElts = map.NumGlobalElements64 ();
#else
  myGblElts = map.MyGlobalElements ();
  numGblElts = map.NumGlobalElements ();
#endif // EXAMPLE_USES_64BIT_GLOBAL_INDICES

  // The number of global indices owned by this MPI process.
  const int numMyElts = map.NumMyElements ();

  // In general, tests like this really should synchronize across all
  // processes.  However, the likely cause for this case is a
  // misconfiguration of Epetra, so we expect it to happen on all
  // processes, if it happens at all.
  if (numMyElts > 0 && myGblElts == NULL) {
    throw std::logic_error ("Failed to get the list of global indices");
  }

  // Local error code for use below.
  int lclerr = 0;

  // Fill the sparse matrix, one row at a time.
  double tempVals[3];
  global_ordinal_type tempGblInds[3];
  for (int i = 0; i < numMyElts; ++i) {
    // A(0, 0:1) = [2, -1]
    if (myGblElts[i] == 0) {
      tempVals[0] = 2.0;
      tempVals[1] = -1.0;
      tempGblInds[0] = myGblElts[i];
      tempGblInds[1] = myGblElts[i] + 1;
      if (lclerr == 0) {
        lclerr = A->InsertGlobalValues (myGblElts[i], 2, tempVals, tempGblInds);
      }
      if (lclerr != 0) {
        break;
      }
    }
    // A(N-1, N-2:N-1) = [-1, 2]
    else if (myGblElts[i] == numGblElts - 1) {
      tempVals[0] = -1.0;
      tempVals[1] = 2.0;
      tempGblInds[0] = myGblElts[i] - 1;
      tempGblInds[1] = myGblElts[i];
      if (lclerr == 0) {
        lclerr = A->InsertGlobalValues (myGblElts[i], 2, tempVals, tempGblInds);
      }
      if (lclerr != 0) {
        break;
      }
    }
    // A(i, i-1:i+1) = [-1, 2, -1]
    else {
      tempVals[0] = -1.0;
      tempVals[1] = 2.0;
      tempVals[2] = -1.0;
      tempGblInds[0] = myGblElts[i] - 1;
      tempGblInds[1] = myGblElts[i];
      tempGblInds[2] = myGblElts[i] + 1;
      if (lclerr == 0) {
        lclerr = A->InsertGlobalValues (myGblElts[i], 3, tempVals, tempGblInds);
      }
      if (lclerr != 0) {
        break;
      }
    }
  }

  // If any process failed to insert at least one entry, throw.
  int gblerr = 0;
  (void) comm.MaxAll (&lclerr, &gblerr, 1);
  if (gblerr != 0) {
    if (A != NULL) {
      delete A;
    }
    throw std::runtime_error ("Some process failed to insert an entry.");
  }

  // Tell the sparse matrix that we are done adding entries to it.
  gblerr = A->FillComplete ();
  if (gblerr != 0) {
    if (A != NULL) {
      delete A;
    }
    std::ostringstream os;
    os << "A->FillComplete() failed with error code " << gblerr << ".";
    throw std::runtime_error (os.str ());
  }

  return A;
}


void
example (const Epetra_Comm& comm)
{
  // The global number of rows in the matrix A to create.  We scale
  // this relative to the number of (MPI) processes, so that no matter
  // how many MPI processes you run, every process will have 10 rows.
  const global_ordinal_type numGblElts = 10 * comm.NumProc ();
  // The global min global index in all the Maps here.
  const global_ordinal_type indexBase = 0;

  // Local error code for use below.
  //
  // In the ideal case, we would use this to emulate behavior like
  // that of Haskell's Maybe in the context of MPI.  That is, if one
  // process experiences an error, we don't want to abort early and
  // cause the other processes to deadlock on MPI communication
  // operators.  Rather, we want to chain along the local error state,
  // until we reach a point where it's natural to pass along that
  // state with other processes.  For example, if one is doing an
  // MPI_Allreduce anyway, it makes sense to pass along one more bit
  // of information: whether the calling process is in a local error
  // state.  Epetra's interface doesn't let one chain the local error
  // state in this way, so we use extra collectives below to propagate
  // that state.  The code below uses very conservative error checks;
  // typical user code would not need to be so conservative and could
  // therefore avoid all the all-reduces.
  int lclerr = 0;

  // Construct a Map that is global (not locally replicated), but puts
  // all the equations on MPI Proc 0.
  const int procZeroMapNumLclElts = (comm.MyPID () == 0) ?
    numGblElts :
    static_cast<global_ordinal_type> (0);
  Epetra_Map procZeroMap (numGblElts, procZeroMapNumLclElts, indexBase, comm);

  // Construct a Map that puts approximately the same number of
  // equations on each processor.
  Epetra_Map globalMap (numGblElts, indexBase, comm);

  // Create a sparse matrix using procZeroMap.
  Epetra_CrsMatrix* A = createCrsMatrix (procZeroMap);
  if (A == NULL) {
    lclerr = 1;
  }

  // Make sure that sparse matrix creation succeeded.  Normally you
  // don't have to check this; we are being extra conservative because
  // this example is also a test.  Even though the matrix's rows live
  // entirely on Process 0, the matrix is nonnull on all processes in
  // its Map's communicator.
  int gblerr = 0;
  (void) comm.MaxAll (&lclerr, &gblerr, 1);
  if (gblerr != 0) {
    throw std::runtime_error ("createCrsMatrix returned NULL on at least one "
                              "process.");
  }

  //
  // We've created a sparse matrix whose rows live entirely on MPI
  // Process 0.  Now we want to distribute it over all the processes.
  //

  // Redistribute the matrix.  Since both the source and target Maps
  // are one-to-one, we could use either an Import or an Export.  If
  // only the source Map were one-to-one, we would have to use an
  // Import; if only the target Map were one-to-one, we would have to
  // use an Export.  We do not allow redistribution using Import or
  // Export if neither source nor target Map is one-to-one.

  // Make an export object with procZeroMap as the source Map, and
  // globalMap as the target Map.  The Export type has the same
  // template parameters as a Map.  Note that Export does not depend
  // on the Scalar template parameter of the objects it
  // redistributes.  You can reuse the same Export for different
  // Tpetra object types, or for Tpetra objects of the same type but
  // different Scalar template parameters (e.g., Scalar=float or
  // Scalar=double).
  Epetra_Export exporter (procZeroMap, globalMap);

  // Make a new sparse matrix whose row map is the global Map.
  Epetra_CrsMatrix B (Copy, globalMap, 0);

  // Redistribute the data, NOT in place, from matrix A (which lives
  // entirely on Proc 0) to matrix B (which is distributed evenly over
  // the processes).
  //
  // Export() has collective semantics, so we must always call it on
  // all processes collectively.  This is why we don't select on
  // lclerr, as we do for the local operations above.
  lclerr = B.Export (*A, exporter, Insert);

  // Make sure that the Export succeeded.  Normally you don't have to
  // check this; we are being extra conservative because this example
  // example is also a test.  We test both min and max, since lclerr
  // may be negative, zero, or positive.
  gblerr = 0;
  (void) comm.MinAll (&lclerr, &gblerr, 1);
  if (gblerr != 0) {
    throw std::runtime_error ("Export() failed on at least one process.");
  }
  (void) comm.MaxAll (&lclerr, &gblerr, 1);
  if (gblerr != 0) {
    throw std::runtime_error ("Export() failed on at least one process.");
  }

  // FillComplete has collective semantics, so we must always call it
  // on all processes collectively.  This is why we don't select on
  // lclerr, as we do for the local operations above.
  lclerr = B.FillComplete ();

  // Make sure that FillComplete succeeded.  Normally you don't have
  // to check this; we are being extra conservative because this
  // example is also a test.  We test both min and max, since lclerr
  // may be negative, zero, or positive.
  gblerr = 0;
  (void) comm.MinAll (&lclerr, &gblerr, 1);
  if (gblerr != 0) {
    throw std::runtime_error ("B.FillComplete() failed on at least one process.");
  }
  (void) comm.MaxAll (&lclerr, &gblerr, 1);
  if (gblerr != 0) {
    throw std::runtime_error ("B.FillComplete() failed on at least one process.");
  }

  if (A != NULL) {
    delete A;
  }
}


int
main (int argc, char *argv[])
{
  using std::cout;
  using std::endl;

#ifdef HAVE_MPI
  MPI_Init (&argc, &argv);
  Epetra_MpiComm comm (MPI_COMM_WORLD);
#else
  Epetra_SerialComm comm;
#endif // HAVE_MPI

  const int myRank = comm.MyPID ();
  const int numProcs = comm.NumProc ();

  if (myRank == 0) {
    // Print out the Epetra software version.
    cout << Epetra_Version () << endl << endl
         << "Total number of processes: " << numProcs << endl;
  }

  example (comm); // Run the whole example.

  // This tells the Trilinos test framework that the test passed.
  if (myRank == 0) {
    cout << "End Result: TEST PASSED" << endl;
  }

#ifdef HAVE_MPI
  (void) MPI_Finalize ();
#endif // HAVE_MPI

  return 0;
}