Skip to content

Use existing memory in fftw vec when doing VecLoad_HDF5 to avoid memory leaks and double-free

Junchao Zhang requested to merge jczhang/fix-vec-hdf5-fftw-memory-error into maint

The MR tries to fix memory leaks and double-free errors when VecLoad_HDF5 on fftw vectors.

fftw vectors have their data array allocated through fftw_malloc (instead of PetscMalloc).

fftw_malloc ... These are functions that behave identically to malloc and free, except that they guarantee that the returned pointer obeys any special alignment restrictions imposed by any algorithm in FFTW (e.g. for SIMD acceleration).

fftw vectors are created with VecCreateMPIWithArray(...,XXX,&v) and have v->data->array = XXX and v->data->array_allocated = NULL. They are destroyed by

PetscErrorCode VecDestroy_MPIFFTW(Vec v)
{
  ...
  PetscFunctionBegin;
  ierr = VecGetArray(v,&array);CHKERRQ(ierr);
  fftw_free((fftw_complex*)array);CHKERRQ(ierr);
  ierr = VecRestoreArray(v,&array);CHKERRQ(ierr);
  ierr = VecDestroy_MPI(v);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

VecLoad_HDF5(xin,viewer) calls PetscViewerHDF5Load, which PetscMallocs a new array x, then calls VecReplaceArray()

  ierr = PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void**)&x);CHKERRQ(ierr);
  ierr = VecSetUp(xin);CHKERRQ(ierr);
  ierr = VecReplaceArray(xin, x);CHKERRQ(ierr);

VecReplaceArray() turns out to be

PetscErrorCode VecReplaceArray_Seq(Vec vin,const PetscScalar *a)
{
  ...
  PetscFunctionBegin;
  ierr = PetscFree(v->array_allocated);CHKERRQ(ierr);
  v->array_allocated = v->array = (PetscScalar*)a;
  PetscFunctionReturn(0);
}

So, if the vector is an fftw vector, we have

PetscErrorCode VecReplaceArray_Seq(Vec vin,const PetscScalar *a)
{
  ...
  PetscFunctionBegin;
  ierr = PetscFree(v->array_allocated);CHKERRQ(ierr); // PetscFree(NULL)
  v->array_allocated = v->array = (PetscScalar*)a; // fftw_malloc'ed memory (v->array) will be lost here.
  PetscFunctionReturn(0);
}

And, when destroying the vector

PetscErrorCode VecDestroy_MPIFFTW(Vec v)
{
  ...
  PetscFunctionBegin;
  ierr = VecGetArray(v,&array);CHKERRQ(ierr);
  fftw_free((fftw_complex*)array);CHKERRQ(ierr); // fftw_free(v->data->array), which was allocated by PetscMalloc
  ierr = VecRestoreArray(v,&array);CHKERRQ(ierr);
  ierr = VecDestroy_MPI(v);CHKERRQ(ierr); // PetscFree(v->data->array_allocated), which was already (incorrectly) freed.
  PetscFunctionReturn(0);
}

The solution is in VecLoad_HDF5 checking whether the vector is established. If yes, reusing existing memory (might be allocated by third party libraries like fftw). I have to add a requirement that the vector loaded from VecLoad_HDF5 has the same length as the established one. Otherwise, I don't know how to properly allocate memory for the vector.

Edited by Junchao Zhang

Merge request reports