On May 11, 2014, at 8:07 PM, LikunTan <[email protected]> wrote:
> Dear Petsc developers,
>
> I have a vector M which consists of a series of 3d vectors, and I want to
> reset M by normalizing each 3d vector. Here is my code:
This code is wrong in several ways
>
> /**********************************************************************
> VecGetArray(M, &aM);
This always returns an array that is indexed starting at 0, so the code
below when you access with aM[node*3+l] is like totally reading from the wrong
place.
Instead use
typedef struct {
PetscScalar x,y,z;
} Field;
Field *aM;
DMDAVecGetArray(da,M,&aM);
This routine returns an array whose index starts at xs
> DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
>
> for(node=xs; node<xs+xm; node++)
> {
> mag=PetscSqrtScalar(aM[xs].x*aM[xs].x + aM[xs].y*aM[xs].y
> aM[xs].z*aM[xs].z);
> if (mag != 0.0) {
> aM[xs].x /= mag;
> aM[xs].y /= mag;
> aM[xs].z /= mag;
}
> }
> DMDAVecRestoreArray(da,M, &aM);
>
When accessing the vector arrays directly you do not need
VecAssemblyBegin/End() they are only for use with VecSetValues()
Also use valgrind to find memory access problems:
http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
> VecView(M, PETSC_VIEWER_STDOUT_WORLD);
> **********************************************************************/
>
> but I got the error at the last step:
> --------------------------------------------------------------------------
> mpiexec noticed that process rank 3 with PID 17156 on node compute-21-8.local
> exited on signal 6 (Aborted).
> --------------------------------------------------------------------------
> and if I commented out VecView, and used vector M for other operations, e.g.
> KSPSolve(ksp, M, b), I got "memory corruption" message. Your comment on this
> issue is well appreciated.
