[Trilinos-Users] Epetra_MultiVector - construction with fixed stride

[Heroux, Mike]

Martin,

This constructor takes in:
- CV: Either defined as View or Copy.  If it is Copy, new storage will be
allocated, and your data will be copied into the new storage.  You can
safely delete your data object afterward.  If it is View, Epetra retains a
pointer to your data.  It is an unmanaged reference, so you are in charge
of guaranteeing the integrity of data.  There are a few ways to do this,
but it requires some care on your part.

- Map: A constructed Epetra_Map or Epetra_BlockMap (most people use
Epetra_Map, which is a specialization of Epetra_Map) The number of rows,
both global and local are defined by the Map.NumGlobalElements() and
Map.NumMyElements(), respectively.

- A pointer to a pre-allocated contiguous section of memory, which is
logically treated as a 2-dimensional Fortran-style array on each MPI
process.    

- MyLDA: The local ³stride² of the array A follows the concept defined for
BLAS and LAPACK.

- NumVectors: The number of columns in the 2D Fortan array, again
following the BLAS convention.

If you want to minimize memory use, you should use the View option, but
then you need to manage the pointer to A yourself.  There are some reasons
why you might do this for performance, especially if you are using an
array in Epetra and in another part of your code and you don¹t want to
copy data back and forth.  Beyond this, it is typically (IMO) better to
use copy mode, unless there is a severe memory use penalty.

The actual physical storage for A is a one-dimensional array, so in your
example, you would have an array of 500 doubles.  How they are interpreted
depends on your values in the Map.NumMyElements(), NumVectors and MyLDA.
Typically we would set MyLDA = 100, NumVectors = 5.

With Fortran storage, column values are contiguously ordered and row
values are MyLDA*8 bytes apart, so accessing row i is done as follows:

 double * curRow = A+i;
 for (j=0; j<NumVectors; ++j) curRow[j] = curRow[j*MyLDA];

This method of accessing a row is simple but not necessarily cache
efficient, since you are accessing data through potentially large memory
strides, so you may consider a different approach if this is a performance
bottleneck.

I hope this is helpful.

Mike

P.S. The Tpetra/Kokkos stack offers a more advanced approach to handling
these issues, but requires fairly advanced understanding C++ templates and
depends on code that is still evolving.  It may be something to keep in
mind for the future, or now if you have a long view of your activities.


On 10/23/14, 10:18 AM, "Martin Vymazal" <martin.vymazal at vki.ac.be> wrote:

>Hello,
>
> could someone please provide a code snippet demonstrating the use of the
>following Epetra_MultiVector constructor
>
>Epetra_MultiVector::Epetra_MultiVector	(Epetra_DataAccess CV,
>							         const Epetra_BlockMap & Map,
>								double * 	A,
>								int 	MyLDA,
>								int 	NumVectors
>								)
>
>Say I would like to store 5 vectors of length 100 aligned in memory
>directly 
>after each other with no gaps. Does this mean that I am responsible for
>allocating an array of doubles with length 5 x 100 to A ? I assume that
>MyLDA 
>should be 100 and NumVectors should be 5.
>
>What should 'CV' be set to? The documentation says:
> Copy - user data will be copied at construction (but then the memory is
>allocated twice - once by me and then again in the constructor - waste of
>space?)
> View - user data will be encapsulated and used throughout the life of
>the 
>object - does this mean that I am responsible for pre-allocating the raw
>array 
>prior to construction of the MultiVector and for deallocating A after the
>corresponding Epetra_MultiVector is destructed? And how would I know if
>that 
>is safe if the Epetra_MultiVector is held by an RCP from Teuchos?
>
>I'm looking for a way of storing the vectors so that they mimic a dense
>matrix 
>and I can efficiently access one row of this matrix. If the stride is
>constant, 
>I only need to know A, the row index and MyLDA. If this is not the
>correct 
>approach, I'm open to suggestions of course. Thank you for any advice.
>
>Best regards,
>
>  Martin Vymazal
>_______________________________________________
>Trilinos-Users mailing list
>Trilinos-Users at software.sandia.gov
>https://software.sandia.gov/mailman/listinfo/trilinos-users