feat(tests/cuda): Use cuda_tdsops_t in CUDA tests.

Author: semi-h

tests/cuda/test_cuda_tridiag.f90

@@ -6,38 +6,30 @@ program test_cuda_tridiag
    use m_common, only: dp, pi
    use m_cuda_common, only: SZ
    use m_cuda_kernels_dist, only: der_univ_dist, der_univ_subs
-   use m_derparams, only: der_1_vv, der_2_vv
+   use m_cuda_tdsops, only: cuda_tdsops_t, cuda_tdsops_init
 
    implicit none
 
    logical :: allpass = .true.
-   real(dp), allocatable, dimension(:,:,:) :: u, du, u_s, u_e
-   real(dp), device, allocatable, dimension(:,:,:) :: u_dev, du_dev
-   real(dp), device, allocatable, dimension(:,:,:) :: u_recv_s_dev, u_recv_e_dev, &
-                                                      u_send_s_dev, u_send_e_dev
+   real(dp), allocatable, dimension(:, :, :) :: u, du
+   real(dp), device, allocatable, dimension(:, :, :) :: u_dev, du_dev
+   real(dp), device, allocatable, dimension(:, :, :) :: &
+      u_recv_s_dev, u_recv_e_dev, u_send_s_dev, u_send_e_dev
 
-   real(dp), device, allocatable, dimension(:,:,:) :: send_s_dev, send_e_dev, &
-                                                      recv_s_dev, recv_e_dev
+   real(dp), device, allocatable, dimension(:, :, :) :: &
+      send_s_dev, send_e_dev, recv_s_dev, recv_e_dev
 
-   real(dp), allocatable, dimension(:,:) :: coeffs_s, coeffs_e
-   real(dp), allocatable, dimension(:) :: coeffs, dist_fr, dist_bc, dist_af, &
-                                          dist_sa, dist_sc
+   type(cuda_tdsops_t) :: tdsops
 
-   real(dp), device, allocatable, dimension(:,:) :: coeffs_s_dev, coeffs_e_dev
-   real(dp), device, allocatable, dimension(:) :: coeffs_dev, &
-                                                  dist_fr_dev, dist_bc_dev, &
-                                                  dist_af_dev, &
-                                                  dist_sa_dev, dist_sc_dev
-
-   integer :: n, n_block, i, j, k, n_halo, n_stencil, n_iters
+   integer :: n, n_block, i, j, k, n_halo, n_iters
    integer :: n_glob
    integer :: nrank, nproc, pprev, pnext, tag1=1234, tag2=1234
-   integer, allocatable :: srerr(:), mpireq(:)
+   integer :: srerr(4), mpireq(4)
    integer :: ierr, ndevs, devnum, memClockRt, memBusWidth
 
    type(dim3) :: blocks, threads
-   real(dp) :: dx, dx2, alfa, norm_du, tol = 1d-8, tstart, tend
-   real(dp) :: achievedBW, deviceBW
+   real(dp) :: dx, dx_per, norm_du, tol = 1d-8, tstart, tend
+   real(dp) :: achievedBW, deviceBW, achievedBWmax, achievedBWmin
 
    call MPI_Init(ierr)
    call MPI_Comm_rank(MPI_COMM_WORLD, nrank, ierr)
@@ -48,50 +40,34 @@ program test_cuda_tridiag
    ierr = cudaGetDeviceCount(ndevs)
    ierr = cudaSetDevice(mod(nrank, ndevs)) ! round-robin
    ierr = cudaGetDevice(devnum)
-   print*, 'I am rank', nrank, 'I am running on device', devnum
-   pnext = modulo(nrank-nproc+1, nproc)
-   pprev = modulo(nrank-1, nproc)
-   print*, 'rank', nrank, 'pnext', pnext, 'pprev', pprev
-   allocate(srerr(nproc), mpireq(nproc))
+
+   !print*, 'I am rank', nrank, 'I am running on device', devnum
+   pnext = modulo(nrank - nproc + 1, nproc)
+   pprev = modulo(nrank - 1, nproc)
 
    n_glob = 512*4
    n = n_glob/nproc
    n_block = 512*512/SZ
-   n_iters = 1000
+   n_iters = 100
 
    allocate(u(SZ, n, n_block), du(SZ, n, n_block))
    allocate(u_dev(SZ, n, n_block), du_dev(SZ, n, n_block))
 
-   dx = 2*pi/n_glob
-   dx2 = dx*dx
+   dx_per = 2*pi/n_glob
+   dx = 2*pi/(n_glob - 1)
 
    do k = 1, n_block
       do j = 1, n
          do i = 1, SZ
-            u(i, j, k) = sin((j-1+nrank*n)*dx)
+            u(i, j, k) = sin((j - 1 + nrank*n)*dx_per)
          end do
       end do
    end do
 
    ! move data to device
    u_dev = u
 
-   ! set up the tridiagonal solver coeffs
-   call der_2_vv(coeffs, coeffs_s, coeffs_e, dist_fr, dist_bc, dist_af, &
-                 dist_sa, dist_sc, n_halo, dx2, n, 'periodic')
-
-   n_stencil = n_halo*2 + 1
-
-   allocate(coeffs_s_dev(n_stencil, n_halo), coeffs_e_dev(n_stencil, n_halo))
-   allocate(coeffs_dev(n_stencil))
-   coeffs_s_dev(:, :) = coeffs_s(:, :); coeffs_e_dev(:, :) = coeffs_e(:, :)
-   coeffs_dev(:) = coeffs(:)
-
-   allocate(dist_fr_dev(n), dist_bc_dev(n), dist_af_dev(n), &
-            dist_sa_dev(n), dist_sc_dev(n))
-   dist_fr_dev(:) = dist_fr(:); dist_bc_dev(:) = dist_bc(:)
-   dist_af_dev(:) = dist_af(:)
-   dist_sa_dev(:) = dist_sa(:); dist_sc_dev(:) = dist_sc(:)
+   n_halo = 4
 
    ! arrays for exchanging data between ranks
    allocate(u_send_s_dev(SZ, n_halo, n_block))
@@ -102,13 +78,17 @@ program test_cuda_tridiag
    allocate(send_s_dev(SZ, 1, n_block), send_e_dev(SZ, 1, n_block))
    allocate(recv_s_dev(SZ, 1, n_block), recv_e_dev(SZ, 1, n_block))
 
+   ! preprocess the operator and coefficient arrays
+   tdsops = cuda_tdsops_init(n, dx_per, operation='second-deriv', &
+                             scheme='compact6')
+
    blocks = dim3(n_block, 1, 1)
    threads = dim3(SZ, 1, 1)
 
    call cpu_time(tstart)
    do i = 1, n_iters
-      u_send_s_dev(:,:,:) = u_dev(:,1:4,:)
-      u_send_e_dev(:,:,:) = u_dev(:,n-n_halo+1:n,:)
+      u_send_s_dev(:, :, :) = u_dev(:, 1:4, :)
+      u_send_e_dev(:, :, :) = u_dev(:, n - n_halo + 1:n, :)
 
       ! halo exchange
       if (nproc == 1) then
@@ -135,8 +115,8 @@ program test_cuda_tridiag
 
       call der_univ_dist<<<blocks, threads>>>( &
          du_dev, send_s_dev, send_e_dev, u_dev, u_recv_s_dev, u_recv_e_dev, &
-         coeffs_s_dev, coeffs_e_dev, coeffs_dev, &
-         n, dist_fr_dev, dist_bc_dev, dist_af_dev &
+         tdsops%coeffs_s_dev, tdsops%coeffs_e_dev, tdsops%coeffs_dev, &
+         n, tdsops%dist_fw_dev, tdsops%dist_bw_dev, tdsops%dist_af_dev &
       )
 
       ! halo exchange for 2x2 systems
@@ -163,36 +143,58 @@ program test_cuda_tridiag
 
       call der_univ_subs<<<blocks, threads>>>( &
          du_dev, recv_s_dev, recv_e_dev, &
-         n, dist_sa_dev, dist_sc_dev &
+         n, tdsops%dist_sa_dev, tdsops%dist_sc_dev &
       )
    end do
-   call cpu_time(tend)
-   print*, 'Total time', tend-tstart
 
-   achievedBW = 6._dp*n_iters*n*n_block*SZ*dp/(tend-tstart)
-   print'(a, f8.3, a)', 'Achieved BW: ', achievedBW/2**30, ' GiB/s'
+   call cpu_time(tend)
+   if (nrank == 0) print*, 'Total time', tend - tstart
+
+   ! BW utilisation and performance checks
+   ! 4 in the first phase, 2 in the second phase, 6 in total
+   achievedBW = 6._dp*n_iters*n*n_block*SZ*dp/(tend - tstart)
+   call MPI_Allreduce(achievedBW, achievedBWmax, 1, MPI_DOUBLE_PRECISION, &
+                      MPI_MAX, MPI_COMM_WORLD, ierr)
+   call MPI_Allreduce(achievedBW, achievedBWmin, 1, MPI_DOUBLE_PRECISION, &
+                      MPI_MIN, MPI_COMM_WORLD, ierr)
+
+   if (nrank == 0) then
+      print'(a, f8.3, a)', 'Achieved BW min: ', achievedBWmin/2**30, ' GiB/s'
+      print'(a, f8.3, a)', 'Achieved BW max: ', achievedBWmax/2**30, ' GiB/s'
+   end if
 
    ierr = cudaDeviceGetAttribute(memClockRt, cudaDevAttrMemoryClockRate, 0)
-   ierr = cudaDeviceGetAttribute(memBusWidth, cudaDevAttrGlobalMemoryBusWidth, 0)
+   ierr = cudaDeviceGetAttribute(memBusWidth, &
+                                 cudaDevAttrGlobalMemoryBusWidth, 0)
    deviceBW = 2*memBusWidth/8._dp*memClockRt*1000
 
-   print'(a, f8.3, a)', 'Device BW:   ', deviceBW/2**30, ' GiB/s'
-   print'(a, f5.2)', 'Effective BW utilization: %', achievedBW/deviceBW*100
+   if (nrank == 0) then
+      print'(a, f8.3, a)', 'Device BW:   ', deviceBW/2**30, ' GiB/s'
+      print'(a, f5.2)', 'Effective BW util min: %', achievedBWmin/deviceBW*100
+      print'(a, f5.2)', 'Effective BW util max: %', achievedBWmax/deviceBW*100
+   end if
 
    ! check error
    du = du_dev
-   norm_du = norm2(u+du)/n_block
-   print*, 'error norm', norm_du
-
-   if ( norm_du > tol ) then
-      allpass = .false.
-      write(stderr, '(a)') 'Check second derivatives... failed'
-   else
-      write(stderr, '(a)') 'Check second derivatives... passed'
+   norm_du = norm2(u + du)
+   norm_du = norm_du*norm_du/n_glob/n_block/SZ
+   call MPI_Allreduce(MPI_IN_PLACE, norm_du, 1, MPI_DOUBLE_PRECISION, &
+                      MPI_SUM, MPI_COMM_WORLD, ierr)
+   norm_du = sqrt(norm_du)
+
+   if (nrank == 0) print*, 'error norm', norm_du
+
+   if (nrank == 0) then
+      if ( norm_du > tol ) then
+         allpass = .false.
+         write(stderr, '(a)') 'Check second derivatives... failed'
+      else
+         write(stderr, '(a)') 'Check second derivatives... passed'
+      end if
    end if
 
    if (allpass) then
-      write(stderr, '(a)') 'ALL TESTS PASSED SUCCESSFULLY.'
+      if (nrank == 0) write(stderr, '(a)') 'ALL TESTS PASSED SUCCESSFULLY.'
    else
       error stop 'SOME TESTS FAILED.'
    end if