使用 DMDAVecGetArray

问题描述

我正在尝试使用具有 2 个自由度的 PETSc 的 DMDA 向量，并使用结构 like in the manual 访问它们。 但是，当我尝试使用 DMDAVecGetArray 时，即使有一个自由度（如下例所示），我也会收到内存 double free or corruption 错误。当我用 DMDAVecGetArray 替换 VecGetArray 时，一切正常。

导致此错误的原因是什么？

---

编译MWE

> gcc -Wall -Wextra -O0 -g $(pkg-config --cflags petsc mpi) -o mwe mwe.c $(pkg-config --libs petsc mpi)

并运行

> 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,&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,&farr)); // with this
    E(DMLocalToGlobal(domain,fg));

    E(VecDestroy(&f));
    E(VecDestroy(&fg));
    E(DMDestroy(&domain));

    E(PetscFinalize());
    return ierr;
}

Debian 靶心； petsc-dev 3.14.4+dfsg1-1； libopenmpi-dev 4.1.0-7

解决方法

好的，找到问题了。这里的主要区别在于，使用 VecGetArray 时使用本地索引对底层数组进行索引，而使用 DMDAVecGetArray 时使用全局索引。

因此，为了修复代码，我不仅要获取数组的大小，还要获取局部向量部分的第一个索引。

使用函数

PetscInt start,size;
E(DMDAGetCorners(domain,&start,NULL,&size,NULL));

然后从 start 循环到 start + size 修复了内存损坏错误。