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

[Matthew Knepley]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20120215/875f3218/attachment.htm>