module sm
contains
 subroutine smooth( a, b, w0, w1, w2, n, m, niters )
  real, dimension(:,:) :: a,b
  real :: w0, w1, w2
  integer :: n, m, niters
  integer :: i, j, iter

  !$acc kernels present(a(:,2:n-1),b)
   do iter = 1,niters
    do i = 2,n-1
     do j = 2,m-1
      a(i,j) = w0 * b(i,j) + &
               w1 * (b(i-1,j) + b(i,j-1) + b(i+1,j) + b(i,j+1)) + &
               w2 * (b(i-1,j-1) + b(i-1,j+1) + b(i+1,j-1) + b(i+1,j+1))
     enddo
    enddo
    do i = 2,n-1
     do j = 2,m-1
      b(i,j) = a(i,j)
     enddo
    enddo
   enddo
  !$acc end kernels
 end subroutine
 subroutine smoothhost( a, b, w0, w1, w2, n, m, niters )
  real, dimension(:,:) :: a,b
  real :: w0, w1, w2
  integer :: n, m, niters
  integer :: i, j, iter
   do iter = 1,niters
    do i = 2,n-1
     do j = 2,m-1
      a(i,j) = w0 * b(i,j) + &
               w1 * (b(i-1,j) + b(i,j-1) + b(i+1,j) + b(i,j+1)) + &
               w2 * (b(i-1,j-1) + b(i-1,j+1) + b(i+1,j-1) + b(i+1,j+1))
     enddo
    enddo
    do i = 2,n-1
     do j = 2,m-1
      b(i,j) = a(i,j)
     enddo
    enddo
   enddo
 end subroutine
end module

program main
 use sm
 use accel_lib
 real,dimension(:,:),allocatable :: aa, bb
 real,dimension(:,:),allocatable :: aahost, bbhost
 real :: w0, w1, w2
 integer :: i,j,n,m
 integer :: errs, args
 real :: dif, tol

 n = 100
 m = 100

 allocate( aa(n,m) )
 allocate( bb(n,m) )
 allocate( aahost(n,m) )
 allocate( bbhost(n,m) )
 do i = 1,n
   do j = 1,m
     aa(i,j) = 0.
     bb(i,j) = i*1000 + j
     aahost(i,j) = 0.
     bbhost(i,j) = i*1000 + j
   enddo
 enddo

 w0 = 0.5
 w1 = 0.3
 w2 = 0.2

!$acc data copyout(aa(1:n,2:m-1)), copy(bb)
 call smooth( aa, bb, w0, w1, w2, n, m, 5 )
!$acc end data
 call smoothhost( aahost, bbhost, w0, w1, w2, n, m, 5 )

 ! check the results
 errs = 0
 tol = 0.000005
 do i = 2,n-1
  do j = 2,m-1
   dif = abs(aa(i,j) - aahost(i,j))
   if( aahost(i,j) .ne. 0 ) dif = abs(dif/aahost(i,j))
   if( dif .gt. tol )then
    errs = errs + 1
    if( errs .le. 10 )then
     print *, i, j, aa(i,j), aahost(i,j)
    endif
   endif
  enddo
 enddo
 print *, errs, ' errors found'
end program