[Trilinos-Users] Handling of ghosts in distributed vectors

[Hoemmen, Mark]

On Apr 19, 2013, at 12:00 PM, <trilinos-users-request at software.sandia.gov>
 wrote:
> Message: 1
> Date: Fri, 19 Apr 2013 14:04:02 +0200
> From: "Bart Janssens" <bart.janssens at lid.kviv.be>
> Subject: [Trilinos-Users] Handling of ghosts in distributed vectors
> To: "trilinos-users at software.sandia.gov"
> 	<trilinos-users at software.sandia.gov>
> Message-ID:
> 	<CAJoBg_VgLifMpFpxGH+R5r9r13vR12hMciZRCnYUjHnJz8Hd2g at mail.gmail.com>
> Content-Type: text/plain; charset="iso-8859-1"
> 
> Hello,
> 
> We are using Epetra_Vectors to store the solution and RHS for our linear
> systems. Up to now, we stored the ghost nodes in the vector (all ghosts at
> the end). However, I am now trying to directly use the Apply method from
> the Epetra_Operator interface, and it seems our vectors are not compatible
> with the domain of the matrix.

That's correct.  The domain Map of an Epetra_CrsMatrix must be one-to-one, which means that if a process owns a global index (row of the vector), no other process may own that index.  Most Epetra users do not need to worry about storing or communicating ghost nodes, because Epetra handles it automatically.  However, it is possible for expert users to handle storing and/or communicating ghost nodes themselves.  If you want to do this, please see the directions below.

> We call FillComplete() without arguments when building the matrix, which
> seems to set both range and domain for the matrix to the rowmap that is
> used at construction, which excludes the ghosts.

If you want to handle storing ghost nodes in the Vectors yourself, but want to let Epetra handle communication of ghost nodes, here is what you do.

1. Begin by creating Vectors with the CrsMatrix's column Map.  

That corresponds to "storing the ghost nodes at the end," since by default, the CrsMatrix creates a column Map exactly like this (with ghost nodes at the end, in increasing order of process rank).  

2. Give "domain Map views" of these column Map Vectors to the Apply() method.

Before you give a Vector to the Apply() method, you would want either to change the Map in place (using the ReplaceMap() method inherited from Epetra_MultiVector), or create a separate View of the original Vector's data using the domain Map, and give the view to the Apply() method.  There are a couple different ways to do the latter.  I recommend the latter approach, as that way you don't have to remember to put the Vector's original Map back after calling Apply().

> If I leave the ghost nodes out of the solution and RHS vectors, both
> Epetra_Operator::Apply and the thyra solve appear to work. So it seems it's
> better not to include ghost nodes in the vectors, is that correct?

See above.  Epetra does not require users to include ghost nodes in the vectors, but it is possible.  I only recommend it for expert users, since it complicates the code, and mistakes can cause undefined behavior (e.g., segfaults).

> This also makes me wonder: how is the data at the ghost nodes -needed for matrix
> application- obtained? Does this involve extra MPI communication?

It involves no extra MPI communication, other than what a distributed-memory sparse matrix-vector multiply necessarily requires (that is, communicating the ghost nodes from the processes that own them).  If you have already done that communication before calling the Apply() method, then none of my comments above apply.  Instead, you could use a "local apply" technique to skip the communication step.  I am not an Epetra expert (I am a Tpetra developer) so I would have to ask a local expert about how to do that with Epetra.

> Also one final question: If we do store ghost nodes in a vector, is there
> any way to force an update of their values through MPI?

You don't normally have to do that for a sparse matrix-vector multiply with Epetra_CrsMatrix, because its Apply() method does that automatically.  However, if you want to do it for some other reason, you can use the matrix's Epetra_Import object.  Call its Importer() method to get the object, then call the Import() method (inherited from Epetra_DistObject) on the column Map Vector, using the domain Map Vector (which may be a view of the column Map Vector, but need not be) as the source of the Import().  For example:

Epetra_CrsMatrix A (yourRowMap, …); // your matrix
// … fill your matrix ...
A.fillComplete (); // otherwise A doesn't have a column Map
Epetra_Vector sourceVector (A.DomainMap ());
Epetra_Vector destinationVector (A.ColMap ());
// … fill sourceVector with your data ...
const int err = destinationVector.Import (sourceVector, * (A.Importer ()), INSERT);
if (err != 0) {
  std::runtime_error ("Epetra Import failed");
}

Please let us know if you have any further questions.

Best,
mfh