device_math_parallel.pf

module device_math_parallel
  use mpi
  use math
  use pfunit
  use device
  use device_math  
  use neko_config
  use num_types
  use comm, only : NEKO_COMM, pe_rank, pe_size
  implicit none

  @TestCase
  type, extends(MPITestCase) :: test_device_math
   contains
     procedure :: setUp
     procedure :: tearDown
  end type test_device_math
  
contains

  subroutine setUp(this)
    class(test_device_math), intent(inout) :: this
    if (NEKO_BCKND_DEVICE .eq. 1) then
       call device_init
    end if
  end subroutine setUp

  subroutine tearDown(this)
    class(test_device_math), intent(inout) :: this
    if (NEKO_BCKND_DEVICE .eq. 1) then
         call device_finalize
   end if
  end subroutine tearDown

  @test
  subroutine test_device_math_rzero(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp) :: a(n)
    type(c_ptr) :: a_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       a = 1.0_rp
       
       call device_map(a, a_d, size(a))
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_rzero(a_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(0.0_rp, a(i), tolerance=NEKO_EPS)
       end do
       
       do i = 1, n
          a(i) = real(i, rp)
       end do
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_rzero(a_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(0.0_rp, a(i), tolerance=NEKO_EPS)
       end do
    end if
    
  end subroutine test_device_math_rzero

  @test
  subroutine test_device_math_rone(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp) :: a(n)
    type(c_ptr) :: a_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       a = 0.0_rp
       
       call device_map(a, a_d, size(a))
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_rone(a_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(1.0_rp, a(i), tolerance=NEKO_EPS)
       end do
       
       do i = 1, n
          a(i) = real(i, rp)
       end do
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_rone(a_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(1.0_rp, a(i), tolerance=NEKO_EPS)
       end do
    end if
    
  end subroutine test_device_math_rone

  @test
  subroutine test_device_math_copy(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp) :: a(n), b(n)
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       
       call rone(a, n)
       call rzero(b, n)
       do i = 1, n
          @assertNotEqual(a(i), b(i))
       end do

       call device_map(a, a_d, size(a))
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
    
       call device_map(b, b_d, size(b))
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_copy(b_d, a_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.false.)
       call device_memcpy(b, b_d, size(a), DEVICE_TO_HOST, sync=.true.)
    
       do i = 1, n
          @assertEqual(a(i), b(i))
       end do
    end if
    
  end subroutine test_device_math_copy

  @test
  subroutine test_device_math_cfill(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp) :: a(n)
    real(kind=rp), parameter :: c = 42.0_rp
    type(c_ptr) :: a_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       a = 1.0_rp
       
       call device_map(a, a_d, size(a))
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_cfill(a_d, c, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i), tolerance=NEKO_EPS)
       end do

    end if
    
  end subroutine test_device_math_cfill

    @test
  subroutine test_device_math_cadd(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp) :: a(n)
    real(kind=rp), parameter :: c = 2.0_rp
    type(c_ptr) :: a_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       a = 40.0_rp
       
       call device_map(a, a_d, size(a))
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_cadd(a_d, c, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i), tolerance=NEKO_EPS)
       end do
       
    end if
    
  end subroutine test_device_math_cadd

  @test
  subroutine test_device_math_add2(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 40.0_rp
          b(i) = 2.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_add2(a_d, b_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i))
       end do
    end if
    
  end subroutine test_device_math_add2
  
  @test
  subroutine test_device_math_add2s1(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    real(kind=rp), parameter :: c1 = 2.0_rp
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 20.0_rp
          b(i) = 2.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_add2s1(a_d, b_d, c1, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i))
       end do
    end if
    
  end subroutine test_device_math_add2s1

  @test
  subroutine test_device_math_add2s2(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    real(kind=rp), parameter :: c1 = 2_rp
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 20.0_rp
          b(i) = 2.0_rp
       end do

       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_add2s2(a_d, b_d, c1, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)

       do i = 1, n
          @assertEqual(24.0_rp, a(i))
       end do
    end if
    
  end subroutine test_device_math_add2s2

  @test
  subroutine test_device_math_add3s2(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, c, b
    real(kind=rp), parameter :: c1 = 2_rp
    real(kind=rp), parameter :: c2 = 4_rp
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    type(c_ptr) :: c_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 0.0_rp
          b(i) = 2.0_rp
          c(i) = 20.0_rp
       end do

       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       call device_map(c, c_d, size(c))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(c, c_d, size(c), HOST_TO_DEVICE, sync=.false.)
       
       call device_add3s2(a_d, b_d, c_d, c1, c2, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)

       do i = 1, n
          @assertEqual(84.0_rp, a(i))
       end do
    end if
    
  end subroutine test_device_math_add3s2

  @test
  subroutine test_device_math_invcol1(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a
    type(c_ptr) :: a_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 2.0_rp
       end do
       
       call device_map(a, a_d, size(a))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       call device_invcol1(a_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(0.5_rp, a(i))
       end do
    end if

  end subroutine test_device_math_invcol1
  
  @test
  subroutine test_device_math_invcol2(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 21.0_rp
          b(i) = 2.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_invcol2(a_d, b_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(10.5_rp, a(i))
       end do
    end if

  end subroutine test_device_math_invcol2
  
  @test
  subroutine test_device_math_col2(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 21.0_rp
          b(i) = 2.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_col2(a_d, b_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i))
       end do
    end if

  end subroutine test_device_math_col2

  @test
  subroutine test_device_math_col3(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b, c
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    type(c_ptr) :: c_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 17.0_rp
          b(i) = 2.0_rp
          c(i) = 21.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       call device_map(c, c_d, size(c))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(c, c_d, size(c), HOST_TO_DEVICE, sync=.false.)
       
       call device_col3(a_d, b_d, c_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i))
       end do
    end if

  end subroutine test_device_math_col3

  @test
  subroutine test_device_math_sub2(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 44.0_rp
          b(i) = 2.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       call device_sub2(a_d, b_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(42.0_rp, a(i))
       end do
    end if

  end subroutine test_device_math_sub2
  
  @test
  subroutine test_device_math_sub3(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b, c
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    type(c_ptr) :: c_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 17.0_rp
          b(i) = 2.0_rp
          c(i) = 21.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       call device_map(c, c_d, size(c))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(c, c_d, size(c), HOST_TO_DEVICE, sync=.false.)
       
       call device_sub3(a_d, b_d, c_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(-19.0_rp, a(i))
       end do
    end if

  end subroutine test_device_math_sub3

  @test
  subroutine test_device_math_addcol3(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b, c
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    type(c_ptr) :: c_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 4669.0_rp
          b(i) = 2.0_rp
          c(i) = 21.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       call device_map(c, c_d, size(c))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(c, c_d, size(c), HOST_TO_DEVICE, sync=.false.)
       
       call device_addcol3(a_d, b_d, c_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(4711.0_rp, a(i))
       end do
    end if

  end subroutine test_device_math_addcol3

  @test
  subroutine test_device_math_addcol4(this)
    class(test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b, c, d
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    type(c_ptr) :: c_d = C_NULL_PTR
    type(c_ptr) :: d_d = C_NULL_PTR
    integer :: i

    if (NEKO_BCKND_DEVICE .eq. 1) then
       do i = 1, n
          a(i) = 4669.0_rp
          b(i) = 2.0_rp
          c(i) = 21.0_rp
          d(i) = 1.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       call device_map(c, c_d, size(c))
       call device_map(d, d_d, size(d))

       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(c, c_d, size(c), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(d, d_d, size(d), HOST_TO_DEVICE, sync=.false.)
       
       call device_addcol4(a_d, b_d, c_d, d_d, n)
       
       call device_memcpy(a, a_d, size(a), DEVICE_TO_HOST, sync=.true.)
       
       do i = 1, n
          @assertEqual(4711.0_rp, a(i))
       end do
    end if

  end subroutine test_device_math_addcol4

  @test(npes=[1])
  subroutine test_device_math_glsc3(this)
    class (test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b, c
    real(kind=rp) :: res, expected
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    type(c_ptr) :: c_d = C_NULL_PTR
    integer :: i, ierr

    if (NEKO_BCKND_DEVICE .eq. 1) then
       
       call MPI_Comm_dup(this%getMpiCommunicator(), NEKO_COMM%mpi_val, ierr)
       
       do i = 1, n
          a(i) = 1.0_rp
          b(i) = 1.0_rp
          c(i) = 1.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       call device_map(c, c_d, size(c))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(c, c_d, size(c), HOST_TO_DEVICE, sync=.false.)
       
       expected = real(this%getNumProcesses(), rp)  * real(n, rp)
       res = device_glsc3(a_d, b_d, c_d, n)
       @assertEqual(expected, res)
       @MPIassertEqual(expected, res)
    end if
    
  end subroutine test_device_math_glsc3

  @test(npes=[1])
  subroutine test_device_math_glsc2(this)
    class (test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a, b
    real(kind=rp) :: res, expected
    type(c_ptr) :: a_d = C_NULL_PTR
    type(c_ptr) :: b_d = C_NULL_PTR
    integer :: i, ierr

    if (NEKO_BCKND_DEVICE .eq. 1) then
       
       call MPI_Comm_dup(this%getMpiCommunicator(), NEKO_COMM%mpi_val, ierr)
       
       do i = 1, n
          a(i) = 1.0_rp
          b(i) = 1.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       call device_map(b, b_d, size(b))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       call device_memcpy(b, b_d, size(b), HOST_TO_DEVICE, sync=.false.)
       
       expected = real(this%getNumProcesses(), rp)  * real(n, rp)
       res = device_glsc2(a_d, b_d, n)
       @assertEqual(expected, res)
       @MPIassertEqual(expected, res)
    end if
    
  end subroutine test_device_math_glsc2

  @test(npes=[1])
  subroutine test_device_math_glsum(this)
    class (test_device_math), intent(inout) :: this
    integer, parameter :: n = 4711
    real(kind=rp), dimension(n) :: a
    real(kind=rp) :: res, expected
    type(c_ptr) :: a_d = C_NULL_PTR
    integer :: i, ierr

    if (NEKO_BCKND_DEVICE .eq. 1) then
       
       call MPI_Comm_dup(this%getMpiCommunicator(), NEKO_COMM%mpi_val, ierr)
       
       do i = 1, n
          a(i) = 1.0_rp
       end do
       
       call device_map(a, a_d, size(a))
       
       call device_memcpy(a, a_d, size(a), HOST_TO_DEVICE, sync=.false.)
       
       expected = real(this%getNumProcesses(), rp)  * real(n, rp)
       res = device_glsum(a_d, n)
       @assertEqual(expected, res)
       @MPIassertEqual(expected, res)
    end if
    
  end subroutine test_device_math_glsum

end module device_math_parallel