Ghost point communication for Red-Black Gauss-Siedel
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
On Wed, Feb 15, 2012 at 5:25 PM, Sylvain Barbot <[email protected]>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/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
-- 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
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 <[email protected]>:
On Wed, Feb 15, 2012 at 5:25 PM, Sylvain Barbot <[email protected]> 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/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
-- 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
On Wed, Feb 15, 2012 at 8:53 PM, Sylvain Barbot <[email protected]>wrote:
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 <[email protected]>:
On Wed, Feb 15, 2012 at 5:25 PM, Sylvain Barbot < [email protected]> 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/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
-- 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
-- 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
R-B has crappy cache performance since you alternately pull data but ignore half of it.
Understood.
Why not just use the GMG in PETSc? Its very flexible, scalable, and efficient.
I tried to use Petsc's implementation of multigrid, cf thread "V-cycle multigrid with matrix shells". My concern is performance. I'd like the whole procedure to use matrix shells, or "matrix-free" matrices, because my stencil can have up to 21 off-diagonal terms. Petsc 3.1 did not allow this functionality. I do not know what Petsc 3.2 has to offer in that regard. I have a working multi-grid method that uses the low functionality of Petsc. I would obviously prefer if I could use something more advanced. Cheers, Sylvain
On Wed, Feb 15, 2012 at 9:14 PM, Sylvain Barbot <[email protected]>wrote:
R-B has crappy cache performance since you alternately pull data but ignore half of it.
Understood.
Why not just use the GMG in PETSc? Its very flexible, scalable, and efficient.
I tried to use Petsc's implementation of multigrid, cf thread "V-cycle multigrid with matrix shells". My concern is performance. I'd like the whole procedure to use matrix shells, or "matrix-free" matrices, because my stencil can have up to 21 off-diagonal terms. Petsc 3.1 did not allow this functionality. I do not know what Petsc 3.2 has to offer in that regard.
So you want to calculate the action to avoid the memory bandwidth limit? If so, you can still use MG by specifying the coarse operators as MatShells, rather than using Galerkin. You might want to tweak the interpolator later if it is inadequate. Matt
I have a working multi-grid method that uses the low functionality of Petsc. I would obviously prefer if I could use something more advanced.
Cheers, Sylvain
-- 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
On Wed, Feb 15, 2012 at 18:25, Sylvain Barbot <[email protected]>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 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.
VecGhost does not work with a DA (pretty obvious since the structured grid layout doesn't put those points near each other). Writing just the updated points into the local array will cost as much as updating the whole thing, so don't worry about it. Your red-black Gauss-Seidel will have horrible cache performance anyway, so no need to sweat about the part that isn't the bottleneck.
participants (3)
-
Jed Brown -
Matthew Knepley -
Sylvain Barbot