On Wed, Feb 15, 2012 at 8:53 PM, Sylvain Barbot <sylbar.vainbot at gmail.com>wrote:
> Hi Matt, > > You're right, I'm not worried about the interior points, that why I'm > looking for updating the ghost points only - the interior points > unaffected by the process. Jed seems to indicate that it's not a > bottle neck. > > Jed, > > Would you recommend a particular step to improve cache performance? > Can you be a bit more specific? > R-B has crappy cache performance since you alternately pull data but ignore half of it. Why not just use the GMG in PETSc? Its very flexible, scalable, and efficient. Matt > Sylvain > > 2012/2/15 Matthew Knepley <knepley at gmail.com>: > > On Wed, Feb 15, 2012 at 5:25 PM, Sylvain Barbot < > sylbar.vainbot at gmail.com> > > wrote: > >> > >> Hi All, > >> > >> I'm implementing an elasticity multi-grid solver with Petsc with > >> matrix shells. I am using the Red-Black Gauss-Siedel smoother. In the > >> middle of the Red-Black passes, I need to update the ghost points > >> between processors. To do that effectively, I'd like to update only > >> the ghost points between the different processors, instead of the > >> whole array. The goal is to do the most possible operations in place > > > > > > What does that mean? Updating ghost dofs means taking the value of that > > dof on the process that owns it, and copying it to all the processes > which > > hold > > that as a ghost dof. This operation has no meaning for interior points. > > > > Matt > > > >> > >> in local array lx, and having to update global vector x only once. The > >> working matrix multiply smoothing operation is shown below. I'm using > >> Petsc 3.1. I read about > >> > >> > http://www.mcs.anl.gov/petsc/petsc-3.1/docs/manualpages/Vec/VecGhostUpdateBegin.html > , > >> but I'm not entirely clear whether this does what I want or not. > >> > >> > >> SUBROUTINE matrixsmooth(A,b,omega,flag,shift,its,level,x,ierr) > >> USE types > >> USE heap > >> USE petscsys > >> USE petscda > >> USE petscis > >> USE petscvec > >> USE petscmat > >> USE petscpc > >> USE petscksp > >> > >> IMPLICIT NONE > >> > >> Mat, INTENT(IN) :: A > >> Vec, INTENT(IN) :: b > >> Vec, INTENT(INOUT) :: x > >> PetscReal, INTENT(IN) :: omega,shift > >> MatSORType, INTENT(IN) :: flag > >> PetscInt, INTENT(IN) :: its, level ! iterations, multi-grid level > >> PetscErrorCode, INTENT(OUT) :: ierr > >> > >> Vec :: lx,lb > >> PetscInt :: istart,iend > >> PetscScalar, POINTER :: px(:) ! pointer to solution array > >> PetscScalar, POINTER :: pb(:) ! pointer to body-force array > >> INTEGER :: rank,isize,i,k > >> INTEGER :: sw,off,p > >> INTEGER :: i2,i3,i2i,i3i, & > >> i000,i0p0,i00p, & > >> i0m0,i00m > >> TYPE(DALOCALINFOF90) :: info > >> > >> CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) > >> CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) > >> > >> CALL VecGetOwnershipRange(x,istart,iend,ierr) > >> > >> ! allocate memory for local vector with ghost points > >> CALL DACreateLocalVector(c%daul(1+level),lx,ierr) > >> CALL VecDuplicate(lx,lb,ierr) > >> > >> ! retrieve forcing term b with ghost points > >> CALL DAGlobalToLocalBegin(c%daul(1+level),b,INSERT_VALUES,lb,ierr) > >> CALL DAGlobalToLocalEnd(c%daul(1+level),b,INSERT_VALUES,lb,ierr) > >> > >> ! obtain a pointer to local data > >> CALL VecGetArrayF90(lx,px,ierr); > >> CALL VecGetArrayF90(lb,pb,ierr); > >> > >> ! geometry info about local vector and padding > >> CALL DAGetLocalInfoF90(c%daul(1+level),info,ierr); > >> > >> ! retrieve stencil width (ghost-point padding) > >> CALL > DAGetInfo(c%daul(1+level),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & > >> > >> > >> > PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, > >> & > >> PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & > >> sw,PETSC_NULL_INTEGER, & > >> PETSC_NULL_INTEGER,ierr); > >> > >> ! offset due to node topology > >> off=MOD(info%xm+info%ym,2) > >> > >> ! fast smoothing (relaxation) its=number of iteration > >> DO k=1,its > >> ! Red-Black Gauss-Siedel scheme > >> DO p=0,1 > >> ! retrieve initial guess x with ghost points > >> CALL > >> DAGlobalToLocalBegin(c%daul(1+level),x,INSERT_VALUES,lx,ierr) > >> CALL > DAGlobalToLocalEnd(c%daul(1+level),x,INSERT_VALUES,lx,ierr) > >> ! smoothing (relaxation) > >> DO i=istart+p+off,iend-1,info%dof*2 > >> i3=(i-istart)/(info%xm*info%dof) > >> i2=(i-istart-i3*info%xm*info%dof)/info%dof > >> i3=i3+i3i ! i3 in ( 0 .. sx3-1 ) > >> i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) > >> > >> i000=1+((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+0)*info%dof > >> > >> i0p0=1+((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1)*info%dof > >> i00p=1+((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0)*info%dof > >> i0m0=1+((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1)*info%dof > >> i00m=1+((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0)*info%dof > >> > >> px(i000+IU1)=px(i000+IU1)+ & > >> (pb(i000+IU1)-((+(px(i0p0+IU1)-px(i000+IU1) > >> )-(px(i000+IU1)-px(i0m0+IU1) )) & > >> > >> +(+(px(i00p+IU1)-px(i000+IU1) )-(px(i000+IU1)-px(i00m+IU1) )))) / > >> (-4._8) > >> END DO > >> > >> ! publish new values of x for its global vector > >> CALL DALocalToGlobal(c%daul(1+level),lx,INSERT_VALUES,x,ierr) > >> > >> END DO > >> END DO > >> > >> ! dispose of the local vectors > >> CALL VecRestoreArrayF90(lx,px,ierr) > >> CALL VecRestoreArrayF90(lb,pb,ierr) > >> CALL DARestoreLocalVector(c%daul(1+level),lx,ierr) > >> CALL DARestoreLocalVector(c%daul(1+level),lb,ierr) > >> > >> END SUBROUTINE matrixsmooth > >> > >> > >> > >> > >> > >> Could this be changed to something like: > >> > >> > >> > >> ... INITIALIZE... > >> > >> CALL DAGlobalToLocalBegin(c%daul(1+c%level),x,INSERT_VALUES,lx,ierr) > >> CALL DAGlobalToLocalEnd(c%daul(1+c%level),x,INSERT_VALUES,lx,ierr) > >> > >> ! fast smoothing (relaxation) its=number of iteration > >> DO k=1,its > >> ! Red-Black Gauss-Siedel scheme > >> DO p=0,1 > >> ! smoothing (relaxation) > >> DO i=istart+p+off,iend-1,info%dof*2 > >> ...STENCIL OPERATION... > >> END DO > >> > >> ...UPDATE GHOST POINTS... > >> > >> END DO > >> END DO > >> > >> ! publish new values of x for its global vector > >> CALL DALocalToGlobal(c%daul(1+level),lx,INSERT_VALUES,x,ierr) > >> > >> ...CLEAN MEMORY... > >> > >> > >> > >> > >> > >> Any recommendation? > >> > >> Best wishes, > >> Sylvain Barbot > > > > > > > > > > -- > > 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 > -- 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 -------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120215/d2d48463/attachment.htm>
