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@gmail.com>:
> On Wed, Feb 15, 2012 at 5:25 PM, Sylvain Barbot <sylbar.vainbot@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