<div class="gmail_quote">On Mon, Oct 31, 2011 at 09:19, Satish Balay <span dir="ltr"><<a href="mailto:balay@mcs.anl.gov">balay@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div id=":1y1">Also - after you are done using the array - you should call<br>
VecRestoreArray(). And you should be changing the values only between<br>
these two calls. [i.e do not stash the pointer for later use - and<br>
keep changing values with this pointer - after the call to<br>
VecRestroeArray()]. If you need to change the vec again - call<br>
VecGetArray() again.</div></blockquote></div><br><div>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.</div>
<div><br></div><div>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:</div>
<div><br></div><div>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>.</div>
<div><br></div><div>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.</div>
<div><br></div><div>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).</div>
<div><br></div><div>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.</div>
<div><br></div><div>PetscErrorCode YourCallback(...,Vec vec) {</div><div>VecAccessor a(vec, access_mode); // Can associate grid information and possibly update ghost values.</div><div>for (int i=a.xs; i<a.xs+a.xm; i++) { // Some grid traversal</div>
<div> a[i] = ...;</div><div>}</div><div>// Destructor calls VecRestoreArray(). It could also accumulate ghost values depending on access_mode.</div>