saxpy.cuf

!**************************************************************************
!
!  Copyright (C) Codeplay Software Ltd.
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!  Codeplay's SYCL-For-CUDA-Examples
!
!  saxpy.cuf
!
!  Description:
!    CUDA Fortran code calling SAXPY from SYCL
!*************************************************************************/

module mathOps
contains
  attributes(global) subroutine saxpy(x, y, a)
    implicit none
    real :: x(:), y(:)
    real, value :: a
    integer :: i, n
    n = size(x)
    i = blockDim%x * (blockIdx%x - 1) + threadIdx%x
    if (i <= n) y(i) = y(i) + a*x(i)
  end subroutine saxpy 
end module mathOps

program testSaxpy
  use mathOps
  use cudafor

  implicit none

interface saxpy_sycl
  subroutine saxpy_call(x, y, a, N) &
    bind(C,name='saxpy_sycl_cuda_wrapper')
    implicit none
    real :: x(:), y(:)
    real, value :: a
    integer, value :: N
  end subroutine
end interface


  integer, parameter :: N = 1024
  real :: x(N), y(N), a
  real, device :: x_d(N), y_d(N)
  type(dim3) :: grid, tBlock

  tBlock = dim3(256,1,1)
  grid = dim3(ceiling(real(N)/tBlock%x),1,1)

  write (*,*) 'CUDA version: '
  x = 1.0; y = 2.0; a = 2.0
  x_d = x
  y_d = y
  call saxpy<<<grid, tBlock>>>(x_d, y_d, a)
  y = y_d
  write(*,*) 'Max error: ', maxval(abs(y-4.0))
  write(*,*) 'N ', N

  write (*,*) 'SYCL version: '
  y = 2.0;
  call saxpy_call(x, y, a, N);
  write(*,*) 'Max error: ', maxval(abs(y-4.0))

end program testSaxpy