I'm trying to use PETSc's DMDA Vectors with 2 degrees of freedom and access them using struct, like in the manual.
However, when I try to use DMDAVecGetArray
even with one degree of freedom (like in example below) I get memory double free or corruption
error. When I replace DMDAVecGetArray
with VecGetArray
everything works just fine.
What causes this error?
Compile MWE with
> gcc -Wall -Wextra -O0 -g $(pkg-config --cflags petsc mpi) -o mwe mwe.c $(pkg-config --libs petsc mpi)
and run with
> mpiexec --host localhost:$(nproc) -n $(nproc) mwe
MWE:
#include <petscdm.h>
#include <petscdmda.h>
#include <petscvec.h>
int
main(int argc, char **argv)
{
PetscErrorCode ierr;
#define E(x) \
do { \
ierr = x; \
CHKERRQ(ierr); \
} while (0)
PetscInt N = 10;
Vec fg, f;
ierr = PetscInitialize(&argc, &argv, NULL, NULL);
if (ierr) {
return ierr;
}
DM domain;
E(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, N, 1, 1, NULL, &domain));
E(DMSetUp(domain));
E(DMCreateGlobalVector(domain, &fg));
E(DMCreateLocalVector(domain, &f));
E(VecZeroEntries(fg));
E(VecAssemblyBegin(fg));
E(VecAssemblyBegin(f));
E(VecAssemblyEnd(fg));
E(VecAssemblyEnd(f));
PetscInt size;
E(VecGetSize(f, &size));
E(DMGlobalToLocal(domain, fg, INSERT_VALUES, f));
PetscScalar *farr;
E(VecGetArray(f, &farr)); // replace this
// E(DMDAVecGetArrayWrite(domain, f, &farr)); // with this
for (PetscInt i = 0; i < size; i++) {
farr[i] += 1;
}
E(VecRestoreArray(f, &farr)); // replace this
// E(DMDAVecRestoreArrayWrite(domain, f, &farr)); // with this
E(DMLocalToGlobal(domain, f, INSERT_VALUES, fg));
E(VecDestroy(&f));
E(VecDestroy(&fg));
E(DMDestroy(&domain));
E(PetscFinalize());
return ierr;
}
Debian bullseye; petsc-dev 3.14.4+dfsg1-1; libopenmpi-dev 4.1.0-7
Ok, found the problem. The key difference here is that when using VecGetArray
underlying array is indexed using local indexes, and when using DMDAVecGetArray
global indexes are used.
So, in order to fix the code I have to get not only the size of the array, but also the first index of the local vector portion.
Using function
PetscInt start, size;
E(DMDAGetCorners(domain, &start, NULL, NULL, &size, NULL, NULL));
and then looping from start
to start + size
fixes the memory corruption error.