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