Kernels

Every source file with kernels starts with #include "fundal.H" (see Macros).

• A nested loop
• A reduction
• An iterative solver with swapped buffers
• A kernel in a procedure
• A procedure working on a mapped array
• Which clause for which build

A nested loop

Collapse two loops into one kernel; the inner loop runs over the first index, contiguous in memory.

!$acc parallel loop collapse(2) DEVICEVAR(a_dev)
!$omp OMPLOOP collapse(2) DEVICEPTR(a_dev)
do j=0, nj-1     ! outer loop: the last index
   do i=0, ni-1  ! inner loop: the first index, contiguous in memory
      a_dev(i,j) = real(10 * i + j, R8P)
   enddo
enddo

$ kernels collapse
  0.0  1.0  2.0
 10.0 11.0 12.0
 20.0 21.0 22.0
 30.0 31.0 32.0
a(i,j) == 10 i + j: T

A reduction

Sum and maximum of a device array: the clause goes before the macros, in both directives.

total = 0._R8P
peak = -huge(1._R8P)
!$acc parallel loop reduction(+:total) reduction(max:peak) DEVICEVAR(a_dev)
!$omp OMPLOOP reduction(+:total) reduction(max:peak) DEVICEPTR(a_dev)
do i=1, n
   total = total + a_dev(i)
   peak = max(peak, a_dev(i))
enddo
print '(A,F9.1,A,F7.1)', 'sum: ', total, ', max: ', peak
print '(A,L1)', 'sum == n (n+1) / 2: ', total == real(n * (n + 1) / 2, R8P)

$ kernels reduction
sum:  500500.0, max:  1000.0
sum == n (n+1) / 2: T

An iterative solver with swapped buffers

Jacobi iterations for the Laplace equation until the largest change is below a tolerance: the new iterate is written to a second buffer, and the two pointers are swapped after each iteration.

do iter=1, max_iter
   change = 0._R8P
   !$acc parallel loop collapse(2) reduction(max:change) DEVICEVAR(t_dev, tn_dev)
   !$omp OMPLOOP collapse(2) reduction(max:change) DEVICEPTR(t_dev, tn_dev)
   do j=1, n
      do i=1, n
         tn_dev(i,j) = 0.25_R8P * (t_dev(i-1,j) + t_dev(i+1,j) + t_dev(i,j-1) + t_dev(i,j+1))
         change = max(change, abs(tn_dev(i,j) - t_dev(i,j)))
      enddo
   enddo
   swap => t_dev ; t_dev => tn_dev ; tn_dev => swap ! the new iterate becomes the current one
   if (change < tol) exit
enddo

$ kernels jacobi
converged in 864 iterations
max |T - 1| < 1e-6: T

WARNING

Without the swap, every iteration recomputes the same tn_dev from the same t_dev and the loop never converges. The swap exchanges the host pointers only; do not copy one array into the other.

A kernel in a procedure

Pass the device arrays as dummy arguments and name the dummies in the clauses, as tutorial chapter 4 does:

subroutine diffuse(n, r, t, tn)
!< One explicit step: the device arrays are dummy arguments, as OpenACC deviceptr requires.
integer(I4P), intent(in)    :: n      ! interior cells
real(R8P),    intent(in)    :: r      ! alpha*dt/dx**2
real(R8P),    intent(in)    :: t(0:)  ! temperature at the current step, device memory
real(R8P),    intent(inout) :: tn(0:) ! temperature at the next step, device memory
integer(I4P)                :: i      ! counter

!$acc parallel loop DEVICEVAR(t, tn)
!$omp OMPLOOP DEVICEPTR(t, tn)
do i=1, n
   tn(i) = t(i) + r * (t(i-1) - 2._R8P * t(i) + t(i+1))
enddo
endsubroutine diffuse

A procedure working on a mapped array

An array mapped by dev_alloc_unstr is found by OpenACC with present; OpenMP needs no clause.

subroutine scale_mapped(a, factor)
!< Scale an array mapped to the device by dev_alloc_unstr: present, not deviceptr.
real(R8P), intent(inout) :: a(:)   ! host array, mapped to the device
real(R8P), intent(in)    :: factor ! scale factor
integer(I4P)             :: i      ! counter

!$acc parallel loop present(a)
!$omp OMPLOOP
do i=1, size(a)
   a(i) = factor * a(i)
enddo
endsubroutine scale_mapped

allocate(a(5)) ; a = [(real(i, R8P), i=1, 5)]
call dev_alloc_unstr(fptr_dev=a)
call dev_memcpy_to_device_unstr(dst=a)
call scale_mapped(a=a, factor=10._R8P)
call dev_memcpy_from_device_unstr(dst=a)
call dev_free_unstr(fptr=a)
print '(A,*(F6.1))', 'a:', a

$ kernels unstructured
a:  10.0  20.0  30.0  40.0  50.0

WARNING

Do not use DEVICEVAR/DEVICEPTR on mapped arrays, nor present on dev_alloc pointers with nvfortran: the first are host arrays with a device copy, the second are device addresses unknown to the present table.

Which clause for which build

Memory	OpenACC, nvfortran	OpenACC, gfortran	OpenMP offload	CPU mode
dev_alloc pointer	DEVICEVAR = deviceptr	DEVICEVAR = present	DEVICEPTR = has_device_addr	shared
dev_alloc_unstr array	present	present	none	none