Is PETSc's DMDA relying on undefined behavior?

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?

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

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.
• 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. – Jed Brown Oct 5 '13 at 0:32
• 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. – jjramsey Oct 5 '13 at 15:41
• 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. – Jed Brown Oct 5 '13 at 16:34