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 <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
