On Fri, Oct 4, 2013 at 1:23 PM, Matthew Knepley <[email protected]> wrote:
> On Fri, Oct 4, 2013 at 6:11 AM, Bishesh Khanal <[email protected]>wrote: > >> Dear all, >> >> I have a part of a code that does NOT use Petsc, and contains an image K >> of dimensions mXnXr. >> It also provides me with a function to access any value at (i,j,k), e.g. >> using K.getVoxelAt(i,j,k). >> >> >> Now in my petsc code I create a DMDA of size mXnXr to solve a pde. The >> coefficients of the pde at each grid point should physically correspond to >> the one in the image K. >> The usual code I would use to create a matrix from DM is: >> >> DMDALocalInfo info; >> ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); >> >> for (PetscInt k = info.zs; k<info.zs+info.zm; ++k) { >> for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) { >> for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) { >> >> //Here (i,j,k) refers to the indexing/ordering that is internally >> organized/distributed to different processors. >> >> } >> } >> } >> >> Now, the problem is that I cannot use getVoxel(i,j,k) within the above >> loop because the ordering and distribution of the image K and the function >> getVoxel(i,j,k) knows nothing about the "da" I create in my part of the >> code. K might have been stored by each processor locally. So I believe if I >> want to use the coefficients from the image K in my code within the above >> loop, first I need to either: >> 1) create global vectors from the DMDA and copy all the values of image K >> to this global vector. Then access the values from this global vector. >> or >> 2) find a way to transform (i,j,k) within the above loop to a valid >> (i,j,k) for the function getVoxel(i,j,k) >> >> What should I do to resolve the problem ? >> > > This explanation is confused. I do not think you mean "global vector", but > rather you mean the K would be completely held by every process. There > is only one question here: > > Is your image K distributed across processes? > No, it is not. > > If so, just match your DA to that distribution. If not, its replicated on > every process and getVoxel() will work fine. > Yes, thanks it works since K is in every process. Sorry, got confused (debugging a code for too long makes sometimes people think in weird ways!) > > Matt > > >> Thanks, >> Bishesh >> >> >> > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener >
