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? Sylvain 2012/2/15 Matthew Knepley <[email protected]>:
On Wed, Feb 15, 2012 at 5:25 PM, Sylvain Barbot <[email protected]> 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/VecGhostUpdateBe..., 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