[petsc-users] Ghost point communication for Red-Black Gauss-Siedel
Sylvain Barbot
sylbar.vainbot at gmail.com
Wed Feb 15 20:53:14 CST 2012
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
More information about the petsc-users
mailing list