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

15 Feb 2012

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