4
$\begingroup$

I was looking at the documentation for DMDAVecGetArray and was surprised that it could create a plain ordinary C array whose indices somehow ranged from, say, istart to istart + size - 1, rather than from 0 to size - 1, or a multidimensional array (also an ordinary C array) whose indices ranged from istart to istart + isize - 1, jstart to jstart + jsize - 1, etc. I wasn't sure how that worked, so I looked into the code and found that DMDAVecGetArray uses the functions VecGetArray1d, VecGetArray2d, etc. behind the scenes to create these oddly indexed arrays. The code for VecGetArray1d looks like this:

PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
 {
   PetscInt       N;

   VecGetLocalSize(x,&N);
   if (m != N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
   VecGetArray(x,a);
   *a  -= mstart;
   return(0);
 }

At least superficially, this looks similar to the trick for doing 1-indexed arrays in C, that is,

int realarray[10];
int *array = &realarray[-1];

except that the pointer a in the PETSc code above is directly decremented instead of having a second pointer variable like array above.

Am I looking at what I think I'm looking at? I find it hard to believe that a well-used library code like PETSc would relying on undefined behavior, so I'm not sure if that PETSc code is really illegal.

What is actually going on?

|cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ Pointer arithmetic isn't illegal in C. In fact, bracket notation and pointer arithmetic are identical. $\endgroup$ – Bill Barth Oct 4 '13 at 13:47
  • $\begingroup$ As long as the pointer isn't dereferenced outside the valid values, there shouldn't be any undefined behavior. $\endgroup$ – Godric Seer Oct 4 '13 at 14:27
  • $\begingroup$ From the link in the question: "the behavior is undefined, even if the pointer is not dereferenced". Is this in the C standard? $\endgroup$ – Dr_Sam Oct 4 '13 at 14:35
  • $\begingroup$ Related: stackoverflow.com/questions/2289548/… $\endgroup$ – Dr_Sam Oct 4 '13 at 14:43
4
$\begingroup$

Heh, fun question. If you read the C standard carefully, you'll find wording like (C99, 6.5.6.8) about pointer arithmetic.

If both the pointer operand and the result point to elements of the same array object, or one past the last element of the array object, the evaluation shall not produce an overflow; otherwise, the behavior is undefined.

As I understand it, the rationale for this restriction is segmented memory models, for which pointers reside in different registers and merely creating a pointer to an invalid memory location (without dereferencing it) could crash or corrupt the program.

We discussed the issue on petsc-dev a few years ago. The conclusion was that segmented memory systems are more-or-less dead and that the "trick" was too useful to give up. If you need to run on a system with segmented memory in which this pointer arithmetic will fail, you should use plain VecGetArray() and do the indexing yourself. You might use C99 VLA-pointers to preserve structured array indexing, but those arrays will necessarily start at zero, so you wouldn't get uniform indexing between the local and global vectors.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ I noticed that the implementation for MatSetValuesStencil and similar functions already had code for translating arrays of MatStencil structs to arrays of local indices, which were then fed to functions such as MatSetValuesLocal. You could get rid of the dependence on that "trick" by putting that code in a new function, say, DMDAStencilToLocalIndices, which could then be used with both matrices and vectors. $\endgroup$ – jjramsey Oct 4 '13 at 21:08
  • $\begingroup$ Sort of, but that's a much less convenient interface and there would have to be a different mapping for the owned part of a global vector and for a local vector (which includes ghost points). Since about 2000 when it was added, we have yet to get a user report of DMDAVecGetArray failing on their system, so we feel comfortable recommending it for general use. $\endgroup$ – Jed Brown Oct 5 '13 at 0:32
  • 1
    $\begingroup$ True, the interface is less convenient. There are probably ways that it could be tweaked, and I'm thinking of trying out some test code following the suggestion you had for using VecGetArray. Anyway, it may be wise for PETSc to have an official alternative to DMDAVecGetArray that doesn't rely on undefined behavior, just in case, even if that alternative isn't as convenient. $\endgroup$ – jjramsey Oct 5 '13 at 15:41
  • $\begingroup$ Yes, the alternative is to use VecGetArray. Fortran callers effectively have to do this anyway. The indexing within a local or owned part is simple and well-defined, just less convenient than the pointer-pointer trick. If you think you have a better way, I suggest writing a simple example that uses that mechanism and sending it as a patch to petsc-dev@mcs.anl.gov or as a pull-request on petsc.git. Thanks. $\endgroup$ – Jed Brown Oct 5 '13 at 16:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.