[petsc-users] Does an array obtained with VecGetArray remain valid after vec's values changed?

[Jed Brown]

On Mon, Oct 31, 2011 at 09:19, Satish Balay <balay at mcs.anl.gov> wrote:

> Also - after you are done using the array - you should call
> VecRestoreArray().  And you should be changing the values only between
> these two calls. [i.e do not stash the pointer for later use - and
> keep changing values with this pointer - after the call to
> VecRestroeArray()].  If you need to change the vec again - call
> VecGetArray() again.
>

This seems to come up frequently. People check that the implementation is
mostly correct when they stash the pointer and then want to do this,
perhaps because of some false sense that it is an "optimization". Don't do
that. Changing the values in this way can invalidate norms and artificially
restricts what can be done with other Vec types or possible future
multicore memory optimizations. Just call VecGetArray() and
VecRestoreArray() around any place that you modify values.

Some people want to wrap PETSc Vecs in their own C++ type that implements
operator[]. I don't recommend wrapping because it makes implementing
callbacks from PETSc more complicated. (You get a plain PETSc Vec and have
to wrap it before you can use it, which means that any extra semantic
information that you wanted to have in the wrapper needs to be
reconstructed (where?). You should be able to compose that extra semantic
information with a PETSc Vec (or DM) so that all the callbacks get that
information.) But if you insist on wrapping and implementing operator[],
you have a few choices:

1. Put VecGetArray() in the operator[] implementation, return a proxy
object with VecRestoreArray() in the proxy's destructor. This is the finest
grain option and I don't like its semantics. To have multiple proxies live
at once, you would need to hold a reference count so that VecRestoreArray()
is only called when the last one goes out of scope. Although unfortunately
a common idiom in C++, returning proxy objects from operator[] is
complicated and unavoidably causes lots of painful semantics, see
std::vector<bool>.

2. Call VecGetArray() when needed in the operator[] implementation and have
an explicit function in your wrapper that ends the access phase by calling
VecRestoreArray(). I don't like this because there isn't a good way to
declare whether you will be modifying the vector or whether access is
read-only (can use VecGetArrayRead()) and I think some explicit statement
that you are about to modify a Vec is good documentation.

3. Have explicit functions in your wrappers to gain and restore access,
with operator[] only allowed between these calls. This has the same
semantics as using the native PETSc API, but it lets you access without
needing to declare local pointers (you still needed your wrappers).

4. Have a function in your wrapper that returns an accessor object that
implements operator[] (in this model, your wrapper does not implement
operator[] itself). The accessor's destructor can call VecRestoreArray(). A
variant of this is to have an accessor whose constructor takes your wrapper
or a PETSc Vec directly. This is my favorite option if you are trying to
avoid raw pointers in your internal interfaces.

PetscErrorCode YourCallback(...,Vec vec) {
VecAccessor a(vec, access_mode); // Can associate grid information and
possibly update ghost values.
for (int i=a.xs; i<a.xs+a.xm; i++) { // Some grid traversal
 a[i] = ...;
}
// Destructor calls VecRestoreArray(). It could also accumulate ghost
values depending on access_mode.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20111031/d8ba46e9/attachment.htm>