[petsc-users] Ghost point communication for Red-Black Gauss-Siedel

Matthew Knepley knepley at gmail.com
Wed Feb 15 20:55:44 CST 2012 

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>

More information about the petsc-users mailing list