Thanks, as always. Interesting that -Wall (with clang) didn’t catch that.
-gideon > On Feb 27, 2017, at 10:28 PM, Barry Smith <[email protected]> wrote: > > > You are using the info struct without ever having set values into it You > are, at a minimum, missing a call to DMDAGetLocalInfo(). Valgrind would > immediately have found the use of uninitialized values. > > Note that DMDAGetCorners() provides some of the information that > DMDAGetLocalInfo() provides but in a slightly different format. > > > Barry > > > >> On Feb 27, 2017, at 9:18 PM, Gideon Simpson <[email protected]> wrote: >> >> Not sure if this is a bug or not, but, for the following code: >> >> #include <petscvec.h> >> #include <petscdm.h> >> #include <petscdmda.h> >> >> typedef struct{ >> PetscScalar u; >> PetscScalar v; >> }b_node; >> >> >> int main(int argc, char **args) >> { >> char dataname[PETSC_MAX_PATH_LEN] = "testvec.bin"; >> PetscInt N = 8; >> DM dm; >> Vec b; >> PetscScalar mass, eng; >> PetscViewer viewer; >> DMDALocalInfo info; >> b_node * b_array; >> Vec b_local; >> PetscInt local_index, local_width, i; >> PetscScalar u, v; >> >> >> PetscInitialize(&argc,&args,NULL,NULL); >> >> PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL); >> >> DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_GHOSTED, N, 2 , 1 , NULL, &dm); >> DMCreateGlobalVector(dm ,&b); >> >> /* Load vector from disk */ >> PetscViewerBinaryOpen(PETSC_COMM_WORLD,dataname, FILE_MODE_READ, &viewer); >> VecLoad(b, viewer); >> PetscViewerDestroy(&viewer); >> >> /* Inspect with viewer */ >> VecView(b, PETSC_VIEWER_STDOUT_WORLD); >> >> DMGetLocalVector(dm,&b_local); >> DMGlobalToLocalBegin(dm,b,INSERT_VALUES,b_local); >> DMGlobalToLocalEnd(dm,b,INSERT_VALUES,b_local); >> >> DMDAVecGetArray(dm, b_local, &b_array); >> DMDAGetCorners (dm, &local_index, NULL, NULL,&local_width, NULL, NULL); >> >> /* Zero out the ghost points */ >> if(info.xs == 0){ >> b_array[-1].u =0.0; >> b_array[-1].v =0.0; >> } >> if(info.xs + info.xm == info.mx){ >> b_array[info.mx].u =0.0; >> b_array[info.mx].v =0.0; >> } >> >> /* Manually inspect */ >> for(i=local_index;i<local_index+local_width;i++){ >> u = b_array[i].u; >> v = b_array[i].v; >> PetscPrintf(PETSC_COMM_WORLD," %i: u = %g, v = %g\n", i, u, v); >> } >> >> DMDAVecRestoreArray (dm,b_local,&b_array); >> DMRestoreLocalVector(dm,&b_local); >> >> VecDestroy(&b); >> DMDestroy(&dm) ; >> >> PetscFinalize(); >> return 0; >> >> } >> >> I am getting the output >> >> Vec Object: 1 MPI processes >> type: seq >> Vec Object:Vec_0x84000000_0 1 MPI processes >> type: mpi >> Process [0] >> 1. >> 0. >> 0.707107 >> 0.707107 >> 6.12323e-17 >> 1. >> -0.707107 >> 0.707107 >> -1. >> 1.22465e-16 >> -0.707107 >> -0.707107 >> -1.83697e-16 >> -1. >> 0.707107 >> -0.707107 >> 0: u = 0., v = 0. >> 1: u = 0.707107, v = 0.707107 >> 2: u = 6.12323e-17, v = 1. >> 3: u = -0.707107, v = 0.707107 >> 4: u = -1., v = 1.22465e-16 >> 5: u = -0.707107, v = -0.707107 >> 6: u = -1.83697e-16, v = -1. >> 7: u = 0.707107, v = -0.707107 >> >> The input vector is an alternating sequence of cos( pi / 4 * j), sin(pi/4 * >> j), for j = 0,1,2…7. What I want to point out is that the first element of >> this vector, as it is read in, is 1. But when we do the manually looping, >> it spits out u = 0 at index i = 0. If i turn off the piece of code that >> zeros out the ghost points, the problem disappears. The reason I am >> concerned is that this kind of looping through the local vector with ghost >> points shows up in a piece of my code, and I am now concerned about the >> output. Do I have a bug here? >> >> >> >> -gideon >> >
