11. On the GPU

probe computes the signed distance of 100000 points from the slanted face inside an OpenACC parallel loop, then compares the result with the method of chapter 5 on the host.

program probe
!< Tutorial, chapter 11: the signed distances of many points to the slanted face, on the GPU.
use penf, only : R8P
use vecfor, only : assign_vector_oac, crossproduct_oac, dotproduct_oac, ex, ey, ez, R8P_mul_vector_oac, vector, &
                   vector_mul_R8P_oac, vector_sub_vector_oac, vector_sum_vector_oac
implicit none
integer, parameter :: n = 100000 ! points
type(vector)       :: p(4)       ! the vertices
type(vector)       :: g(n)       ! the points
type(vector)       :: s12, s13   ! two edges of the slanted face
type(vector)       :: nrm        ! its normal
type(vector)       :: r          ! from the face to a point
type(vector)       :: t          ! a temporary
real(R8P)          :: d(n)       ! signed distances, on the device
real(R8P)          :: h(n)       ! signed distances, on the host
integer            :: i          ! counter

p(2) = ex ; p(3) = ey ; p(4) = ez
do i = 1, n
   g(i) = (real(mod(i, 97), R8P) * ex + real(mod(i, 89), R8P) * ey + real(mod(i, 83), R8P) * ez) / 50
enddo
!$acc parallel loop copyin(g, p) copyout(d) private(s12, s13, nrm, r, t)
do i = 1, n
   call vector_sub_vector_oac(p(3), p(2), s12)
   call vector_sub_vector_oac(p(4), p(2), s13)
   call crossproduct_oac(s12, s13, nrm)
   call vector_mul_R8P_oac(nrm, 1._R8P / sqrt(dotproduct_oac(nrm, nrm)), t)   ! the unit normal
   call assign_vector_oac(nrm, t)
   call R8P_mul_vector_oac(-1._R8P, p(2), t)
   call vector_sum_vector_oac(g(i), t, r)                                      ! r = g(i) - p(2)
   d(i) = dotproduct_oac(nrm, r)
enddo
h = g%distance_to_plane(pt1=p(2), pt2=p(3), pt3=p(4))
print '(A,I0,A)', 'distances of ', n, ' points'
print '(A,I0)',    'points outside            ', count(d > 0._R8P)
print '(A,ES9.1)', 'largest difference to host', maxval(abs(d - h))
endprogram probe

$ probe
distances of 100000 points
points outside            96820
largest difference to host  0.0E+00

The output above comes from gfortran without -fopenacc: the !$acc lines are comments, and the loop runs on the CPU. Built with NVIDIA nvfortran and -acc (the FoBiS mode tests-nvf-acc), the same source runs on the GPU.

Why a separate API

The operators and the methods of the vector cannot be called inside an OpenACC device routine with nvfortran (26.1, the compiler VecFor is validated on): the type-bound operators dispatch through a table of procedures, an indirect call the device does not support; the passed object of a method is polymorphic, with the same problem; the assignment between two vectors is a user-defined assignment, rejected on the device; and a function returning a vector made every thread of a loop share one temporary. So VecFor exposes a few operations as plain procedures, each an !$acc routine seq, that take and return type(vector) arguments:

Operation	Operator	Device routine
sum	c = a + b	call vector_sum_vector_oac(a, b, c)
difference	c = a - b	call vector_sub_vector_oac(a, b, c)
number times vector	c = s * a	call R8P_mul_vector_oac(s, a, c)
vector times number	c = a * s	call vector_mul_R8P_oac(a, s, c)
cross product	c = a .cross. b	call crossproduct_oac(a, b, c)
dot product	d = a .dot. b	d = dotproduct_oac(a, b)
assignment	a = b	call assign_vector_oac(a, b)

The result of a vector operation is the last argument of a subroutine; only the dot product, a number, is a function. The scalar of the products is real(R8P) whatever the kind of the vector. Each kind has its own set: crossproduct_R4P_oac, dotproduct_R16P_oac, ...; the names without a suffix are those of the default vector.

Writing a kernel

The loop body of probe is the method distance_to_plane written out with these seven routines: two edges, their cross product, the unit normal by a product with the inverse of its length, the vector from the face to the point, and its dot product with the normal. The temporaries are private to each iteration. Here the differences with the method are zero; in general, expect differences of the order of the last digit, since the two compute in a different order.

The OpenACC API covers what its first user, a distance kernel, needed: the other operators and methods, the comparisons, the rotations and mirrors are not device-callable (see Device-callable API).

What you learned

The seven _oac routines, their calling convention, an OpenACC loop that uses them; why the operators cannot run on the device. Reference: Device-callable API.

This is the end of the tutorial. The cookbook has a short recipe for each task, the reference every operator, method and argument.