[Trilinos-Users] Epetra_MultiVector - construction with fixed stride

[Martin Vymazal]

Hello Mike,

 this is definitely helpful. Thank you for your reply.

 Martin


On Friday 24 of October 2014 20:53:33 Heroux, Mike wrote:
> 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