[petsc-users] Problems with VecGetArray under sub-communicators

[Smith, Barry F.]

   DMDAVecGetArray() uses global indices while VecGetArray() uses local indices. 

  Barry


> On Apr 22, 2018, at 2:34 PM, Zin Lin <zinlin.zinlin at gmail.com> wrote:
> 
> Thanks a lot!
> That clears things up.
> I do have one more question though.
> The behavior of DMDAVecGetArray is different then? From ex3.c (below), it doesn't seem like the raw array needs the local indices?
> 
> From ex3.c, 
>   c(cda,global,&coors);
> ...
>   DMDAGetCorners(cda,&mstart,&nstart,&pstart,&m,&n,&p);
> ...
>   for (i=mstart; i<mstart+m; i++) {
>     for (j=nstart; j<nstart+n; j++) {
>       for (k=pstart; k<pstart+p; k++) {
>         ...
>           coors[k][j][i].x = ...
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> On Sun, Apr 22, 2018 at 3:14 PM, Dave May <dave.mayhem23 at gmail.com> wrote:
> 
> On Sun, 22 Apr 2018 at 20:13, Zin Lin <zinlin.zinlin at gmail.com> wrote:
> Hi 
> I am experiencing possible memory issues with VecGetArray when it is used under sub-communicators (when I split the PETSC_COMM_WORLD to multiple subcomms). The following is the minimal code. Basically, you can test that if you parallelize the vector to more than one processor under a subcomm, the array obtained from the VecGetArray call doesn't seem to be correct.
> Please test it with 
> 
> 1)  mpirun -np 1 ./example -ncomms 1
> 2)  mpirun -np 2 ./example -ncomms 2 
> 3)  mpirun -np 2 ./example -ncomms 1   
> 4)  mpirun -np 4 ./example -ncomms 2
> 
> you will 1) and 2) work as expected while in 3) and 4) some entries of the array are assigned erroneous values.
> 
> First - your call to PetscPrintf contains a mistake. The second instance of i is missing.
> 
> Second (more important), you cannot access the raw array obtained via VecGetArray() using global indices. You must use local indices.
> Change the access to _u[i-ns] and the code should work.
> 
> Also, debugging based on what printf() tells you can be misleading. It's much better to use valgrind - see here 
> 
> https://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
> 
> 
> Thanks,
>   Dave
> 
> 
> Any input will be appreciated.
> Thanks
> Zin
> 
> Minimal Code
> 
> PetscErrorCode main(int argc, char **argv)
> {
> 
>   MPI_Init(NULL, NULL);
>   PetscInitialize(&argc,&argv,NULL,NULL);
>   int size;
>   MPI_Comm_size(MPI_COMM_WORLD, &size);
>   PetscPrintf(PETSC_COMM_WORLD,"\tThe total number of processors is %d\n",size);
> 
>   //check if the number of processors is divisible by the number of subcomms
>   int ncomms, np_per_comm;
>   PetscOptionsGetInt(NULL,"-ncomms",&ncomms,NULL);
>   if(!(size%ncomms==0)) SETERRQ(PETSC_COMM_WORLD,1,"The number of processes must be a multiple of ncomms so that it is divisible by the number of subcomms.");
>   np_per_comm=size/ncomms;
>   PetscPrintf(PETSC_COMM_WORLD,"\tThe number of subcomms is %d.\n\tEach subcomm has %d processors.\n",ncomms,np_per_comm);
> 
>   //calculate the colour of each subcomm ( = rank of each processor / number of processors in each subcomm )
>   //note once calculated, the colour is fixed throughout the entire run
>   int rank;
>   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
>   MPI_Comm subcomm;
>   int colour = rank/np_per_comm;
>   MPI_Comm_split(MPI_COMM_WORLD, colour, rank, &subcomm);
> 
>   Vec u;
>   PetscScalar *_u;
>   int i,ns,ne;
>   PetscScalar tmp;
> 
>   VecCreateMPI(subcomm,PETSC_DECIDE,10,&u);
>   VecSet(u,1.0+PETSC_i*0);
> 
>   VecGetArray(u,&_u);
>   VecGetOwnershipRange(u,&ns,&ne);
>   for(i=ns;i<ne;i++){
>     VecGetValues(u,1,&i,&tmp);
>     PetscPrintf(PETSC_COMM_SELF,"colour %d, u[%d]_array = %g + i * (%g), u[%d]_vec = %g + i * %g \n",colour,i,creal(_u[i]),cimag(_u[i]),
> 
> Insert i here
> 
> creal(tmp),cimag(tmp));
>   }
>   VecRestoreArray(u,&_u);
> 
>   PetscFinalize();
>   return 0;
> 
> }
> 
> 
> 
> 
> -- 
> Zin Lin
>