pipecg_device.F90

Go to the documentation of this file.

! Copyright (c) 2021-2024, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module pipecg_device
  use krylov, only : ksp_t, ksp_monitor_t, ksp_max_iter
  use precon,  only : pc_t
  use ax_product, only : ax_t
  use num_types, only: rp, c_rp
  use field, only : field_t
  use coefs, only : coef_t
  use gather_scatter, only : gs_t, gs_op_add
  use bc_list, only : bc_list_t
  use math, only : glsc3, rzero, copy, abscmp
  use device_math, only : device_rzero, device_copy, &
       device_glsc3, device_vlsc3
  use device
  use utils, only : neko_error
  use comm, only : neko_comm, pe_size, mpi_iallreduce, mpi_status, &
       mpi_real_precision, mpi_sum, mpi_in_place, mpi_request, &
       mpi_wait
  use, intrinsic :: iso_c_binding, only :  c_ptr, c_null_ptr, &
       c_associated, c_size_t, c_sizeof, c_int, c_loc
  implicit none
  private
 
  integer, parameter :: device_pipecg_p_space = 10
 
  type, public, extends(ksp_t) :: pipecg_device_t
     real(kind=rp), allocatable :: p(:)
     real(kind=rp), allocatable :: q(:)
     real(kind=rp), allocatable :: r(:)
     real(kind=rp), allocatable :: s(:)
     real(kind=rp), allocatable :: u(:,:)
     real(kind=rp), allocatable :: w(:)
     real(kind=rp), allocatable :: z(:)
     real(kind=rp), allocatable :: mi(:)
     real(kind=rp), allocatable :: ni(:)
     real(kind=rp), allocatable :: alpha(:)
     real(kind=rp), allocatable :: beta(:)
     type(c_ptr) :: p_d = c_null_ptr
     type(c_ptr) :: q_d = c_null_ptr
     type(c_ptr) :: r_d = c_null_ptr
     type(c_ptr) :: s_d = c_null_ptr
     type(c_ptr) :: u_d_d = c_null_ptr
     type(c_ptr) :: w_d = c_null_ptr
     type(c_ptr) :: z_d = c_null_ptr
     type(c_ptr) :: mi_d = c_null_ptr
     type(c_ptr) :: ni_d = c_null_ptr
     type(c_ptr) :: alpha_d = c_null_ptr
     type(c_ptr) :: beta_d = c_null_ptr
     type(c_ptr), allocatable :: u_d(:)
     type(c_ptr) :: gs_event = c_null_ptr
   contains
     procedure, pass(this) :: init => pipecg_device_init
     procedure, pass(this) :: free => pipecg_device_free
     procedure, pass(this) :: solve => pipecg_device_solve
     procedure, pass(this) :: solve_coupled => pipecg_device_solve_coupled
  end type pipecg_device_t
 
#ifdef HAVE_CUDA
  interface
     subroutine cuda_pipecg_vecops(p_d, q_d, r_d, s_d, u_d1, u_d2, &
          w_d, z_d, ni_d, mi_d, alpha, beta, mult_d, reduction,n) &
          bind(c, name='cuda_pipecg_vecops')
       use, intrinsic :: iso_c_binding
       import c_rp
       implicit none
       type(c_ptr), value :: p_d, q_d, r_d, s_d, u_d1, u_d2
       type(c_ptr), value :: w_d, ni_d, mi_d, z_d, mult_d
       integer(c_int) :: n
       real(c_rp) :: alpha, beta, reduction(3)
     end subroutine cuda_pipecg_vecops
  end interface
 
  interface
     subroutine cuda_cg_update_xp(x_d, p_d, u_d_d, alpha, beta, &
          p_cur, p_space, n) &
          bind(c, name='cuda_cg_update_xp')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x_d, p_d, u_d_d, alpha, beta
       integer(c_int) :: p_cur, n, p_space
     end subroutine cuda_cg_update_xp
  end interface
#elif HAVE_HIP
  interface
     subroutine hip_pipecg_vecops(p_d, q_d, r_d, s_d, u_d1, u_d2, &
          w_d, z_d, ni_d, mi_d, alpha, beta, mult_d, reduction,n) &
          bind(c, name='hip_pipecg_vecops')
       use, intrinsic :: iso_c_binding
       import c_rp
       implicit none
       type(c_ptr), value :: p_d, q_d, r_d, s_d, u_d1, u_d2
       type(c_ptr), value :: w_d, ni_d, mi_d, z_d, mult_d
       integer(c_int) :: n
       real(c_rp) :: alpha, beta, reduction(3)
     end subroutine hip_pipecg_vecops
  end interface
 
  interface
     subroutine hip_cg_update_xp(x_d, p_d, u_d_d, alpha, beta, &
          p_cur, p_space, n) &
          bind(c, name='hip_cg_update_xp')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: x_d, p_d, u_d_d, alpha, beta
       integer(c_int) :: p_cur, n, p_space
     end subroutine hip_cg_update_xp
  end interface
#endif
 
contains
 
  subroutine device_pipecg_vecops(p_d, q_d, r_d, s_d, u_d1, u_d2, &
       w_d, z_d, ni_d, mi_d, alpha, beta, mult_d, reduction,n)
    type(c_ptr), value :: p_d, q_d, r_d, s_d, u_d1, u_d2
    type(c_ptr), value :: w_d, ni_d, mi_d, z_d, mult_d
    integer(c_int) :: n
    real(c_rp) :: alpha, beta, reduction(3)
#ifdef HAVE_HIP
    call hip_pipecg_vecops(p_d, q_d, r_d,&
         s_d, u_d1, u_d2, w_d, z_d, ni_d, mi_d, alpha, beta, mult_d, reduction,n)
#elif HAVE_CUDA
    call cuda_pipecg_vecops(p_d, q_d, r_d,&
         s_d, u_d1, u_d2, w_d, z_d, ni_d, mi_d, alpha, beta, mult_d, reduction,n)
#else
    call neko_error('No device backend configured')
#endif
  end subroutine device_pipecg_vecops
 
  subroutine device_cg_update_xp(x_d, p_d, u_d_d, alpha, beta, p_cur, p_space, n)
    use, intrinsic :: iso_c_binding
    type(c_ptr), value :: x_d, p_d, u_d_d, alpha, beta
    integer(c_int) :: p_cur, n, p_space
#ifdef HAVE_HIP
    call hip_cg_update_xp(x_d, p_d, u_d_d, alpha, beta, p_cur, p_space, n)
#elif HAVE_CUDA
    call cuda_cg_update_xp(x_d, p_d, u_d_d, alpha, beta, p_cur, p_space, n)
#else
    call neko_error('No device backend configured')
#endif
  end subroutine device_cg_update_xp
 
  subroutine pipecg_device_init(this, n, max_iter, M, rel_tol, abs_tol, monitor)
    class(pipecg_device_t), target, intent(inout) :: this
    class(pc_t), optional, intent(in), target :: M
    integer, intent(in) :: n
    integer, intent(in) :: max_iter
    real(kind=rp), optional, intent(in) :: rel_tol
    real(kind=rp), optional, intent(in) :: abs_tol
    logical, optional, intent(in) :: monitor
    type(c_ptr) :: ptr
    integer(c_size_t) :: u_size
    integer :: i
 
    call this%free()
 
    allocate(this%p(n))
    allocate(this%q(n))
    allocate(this%r(n))
    allocate(this%s(n))
    allocate(this%u(n, device_pipecg_p_space+1))
    allocate(this%u_d(device_pipecg_p_space+1))
    allocate(this%w(n))
    allocate(this%z(n))
    allocate(this%mi(n))
    allocate(this%ni(n))
    allocate(this%alpha(device_pipecg_p_space))
    allocate(this%beta(device_pipecg_p_space))
 
    if (present(m)) then
       this%M => m
    end if
 
    call device_map(this%p, this%p_d, n)
    call device_map(this%q, this%q_d, n)
    call device_map(this%r, this%r_d, n)
    call device_map(this%s, this%s_d, n)
    call device_map(this%w, this%w_d, n)
    call device_map(this%z, this%z_d, n)
    call device_map(this%mi, this%mi_d, n)
    call device_map(this%ni, this%ni_d, n)
    call device_map(this%alpha, this%alpha_d, device_pipecg_p_space)
    call device_map(this%beta, this%beta_d, device_pipecg_p_space)
    do i = 1, device_pipecg_p_space+1
       this%u_d(i) = c_null_ptr
       call device_map(this%u(:,i), this%u_d(i), n)
    end do
    !Did not work with 4 for some reason...
    u_size = 8*(device_pipecg_p_space+1)
    call device_alloc(this%u_d_d, u_size)
    ptr = c_loc(this%u_d)
    call device_memcpy(ptr,this%u_d_d, u_size, &
         host_to_device, sync=.false.)
 
    if (present(rel_tol) .and. present(abs_tol) .and. present(monitor)) then
       call this%ksp_init(max_iter, rel_tol, abs_tol, monitor = monitor)
    else if (present(rel_tol) .and. present(abs_tol)) then
       call this%ksp_init(max_iter, rel_tol, abs_tol)
    else if (present(monitor) .and. present(abs_tol)) then
       call this%ksp_init(max_iter, abs_tol = abs_tol, monitor = monitor)
    else if (present(rel_tol) .and. present(monitor)) then
       call this%ksp_init(max_iter, rel_tol, monitor = monitor)
    else if (present(rel_tol)) then
       call this%ksp_init(max_iter, rel_tol = rel_tol)
    else if (present(abs_tol)) then
       call this%ksp_init(max_iter, abs_tol = abs_tol)
    else if (present(monitor)) then
       call this%ksp_init(max_iter, monitor = monitor)
    else
       call this%ksp_init(max_iter)
    end if
 
    call device_event_create(this%gs_event, 2)
 
  end subroutine pipecg_device_init
 
  subroutine pipecg_device_free(this)
    class(pipecg_device_t), intent(inout) :: this
    integer :: i
 
    call this%ksp_free()
 
    if (allocated(this%p)) then
       deallocate(this%p)
    end if
    if (allocated(this%q)) then
       deallocate(this%q)
    end if
    if (allocated(this%r)) then
       deallocate(this%r)
    end if
    if (allocated(this%s)) then
       deallocate(this%s)
    end if
    if (allocated(this%u)) then
       deallocate(this%u)
    end if
    if (allocated(this%w)) then
       deallocate(this%w)
    end if
    if (allocated(this%z)) then
       deallocate(this%z)
    end if
    if (allocated(this%mi)) then
       deallocate(this%mi)
    end if
    if (allocated(this%ni)) then
       deallocate(this%ni)
    end if
    if (allocated(this%alpha)) then
       deallocate(this%alpha)
    end if
    if (allocated(this%beta)) then
       deallocate(this%beta)
    end if
 
 
    if (c_associated(this%p_d)) then
       call device_free(this%p_d)
    end if
    if (c_associated(this%q_d)) then
       call device_free(this%q_d)
    end if
    if (c_associated(this%r_d)) then
       call device_free(this%r_d)
    end if
    if (c_associated(this%s_d)) then
       call device_free(this%s_d)
    end if
    if (c_associated(this%u_d_d)) then
       call device_free(this%u_d_d)
    end if
    if (c_associated(this%w_d)) then
       call device_free(this%w_d)
    end if
    if (c_associated(this%z_d)) then
       call device_free(this%z_d)
    end if
    if (c_associated(this%mi_d)) then
       call device_free(this%mi_d)
    end if
    if (c_associated(this%ni_d)) then
       call device_free(this%ni_d)
    end if
    if (c_associated(this%alpha_d)) then
       call device_free(this%alpha_d)
    end if
    if (c_associated(this%beta_d)) then
       call device_free(this%beta_d)
    end if
    if (allocated(this%u_d)) then
       do i = 1, device_pipecg_p_space
          if (c_associated(this%u_d(i))) then
             call device_free(this%u_d(i))
          end if
       end do
    end if
 
    nullify(this%M)
 
    if (c_associated(this%gs_event)) then
       call device_event_destroy(this%gs_event)
    end if
 
  end subroutine pipecg_device_free
 
  function pipecg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results)
    class(pipecg_device_t), intent(inout) :: this
    class(ax_t), intent(in) :: ax
    type(field_t), intent(inout) :: x
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: f
    type(coef_t), intent(inout) :: coef
    type(bc_list_t), intent(inout) :: blst
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t) :: ksp_results
    integer, optional, intent(in) :: niter
    integer :: iter, max_iter, ierr, p_cur, p_prev, u_prev
    real(kind=rp) :: rnorm, rtr, reduction(3), norm_fac
    real(kind=rp) :: gamma1, gamma2, delta
    real(kind=rp) :: tmp1, tmp2, tmp3
    type(mpi_request) :: request
    type(mpi_status) :: status
    type(c_ptr) :: f_d
    f_d = device_get_ptr(f)
 
    if (present(niter)) then
       max_iter = niter
    else
       max_iter = this%max_iter
    end if
    norm_fac = 1.0_rp / sqrt(coef%volume)
 
    associate(p => this%p, q => this%q, r => this%r, s => this%s, &
         u => this%u, w => this%w, z => this%z, mi => this%mi, ni => this%ni, &
         alpha => this%alpha, beta => this%beta,  &
         alpha_d => this%alpha_d, beta_d => this%beta_d, &
         p_d => this%p_d, q_d => this%q_d, r_d => this%r_d, &
         s_d => this%s_d, u_d => this%u_d, u_d_d => this%u_d_d, &
         w_d => this%w_d, z_d => this%z_d, mi_d => this%mi_d, ni_d => this%ni_d)
 
      p_prev = device_pipecg_p_space !this%p_space
      u_prev = device_pipecg_p_space + 1 !this%p_space+1
      p_cur = 1
      call device_rzero(x%x_d, n)
      call device_rzero(z_d, n)
      call device_rzero(q_d, n)
      call device_rzero(p_d, n)
      call device_rzero(s_d, n)
      call device_copy(r_d, f_d, n)
      !apply u=M^-1r
      !call device_copy(u_d(u_prev), r_d, n)
      call this%M%solve(u(1,u_prev), r, n)
      call ax%compute(w, u(1,u_prev), coef, x%msh, x%Xh)
      call gs_h%op(w, n, gs_op_add, this%gs_event)
      call device_event_sync(this%gs_event)
      call blst%apply_scalar(w, n)
 
      rtr = device_glsc3(r_d, coef%mult_d, r_d, n)
      rnorm = sqrt(rtr)*norm_fac
      ksp_results%res_start = rnorm
      ksp_results%res_final = rnorm
      ksp_results%iter = 0
      if(abscmp(rnorm, 0.0_rp)) return
 
      gamma1 = 0.0_rp
      tmp1 = 0.0_rp
      tmp2 = 0.0_rp
      tmp3 = 0.0_rp
      tmp1 = device_vlsc3(r_d, coef%mult_d, u_d(u_prev), n)
      tmp2 = device_vlsc3(w_d, coef%mult_d, u_d(u_prev), n)
      tmp3 = device_vlsc3(r_d, coef%mult_d, r_d, n)
      reduction(1) = tmp1
      reduction(2) = tmp2
      reduction(3) = tmp3
 
      call this%monitor_start('PipeCG')
      do iter = 1, max_iter
         call mpi_iallreduce(mpi_in_place, reduction, 3, &
              mpi_real_precision, mpi_sum, neko_comm, request, ierr)
 
         call this%M%solve(mi, w, n)
         call ax%compute(ni, mi, coef, x%msh, x%Xh)
         call gs_h%op(ni, n, gs_op_add, this%gs_event)
         call device_event_sync(this%gs_event)
         call blst%apply(ni, n)
 
         call mpi_wait(request, status, ierr)
         gamma2 = gamma1
         gamma1 = reduction(1)
         delta = reduction(2)
         rtr = reduction(3)
 
         rnorm = sqrt(rtr)*norm_fac
         call this%monitor_iter(iter, rnorm)
         if (rnorm .lt. this%abs_tol) exit
 
 
         if (iter .gt. 1) then
            beta(p_cur) = gamma1 / gamma2
            alpha(p_cur) = gamma1 / (delta - (beta(p_cur) * gamma1/alpha(p_prev)))
         else
            beta(p_cur) = 0.0_rp
            alpha(p_cur) = gamma1/delta
         end if
 
         call device_pipecg_vecops(p_d, q_d, r_d,&
              s_d, u_d(u_prev), u_d(p_cur),&
              w_d, z_d, ni_d,&
              mi_d, alpha(p_cur), beta(p_cur),&
              coef%mult_d, reduction, n)
         if (p_cur .eq. device_pipecg_p_space) then
            call device_memcpy(alpha, alpha_d, p_cur, &
                 host_to_device, sync=.false.)
            call device_memcpy(beta, beta_d, p_cur, &
                 host_to_device, sync=.false.)
            call device_cg_update_xp(x%x_d, p_d, u_d_d, alpha_d, beta_d, p_cur, &
                 device_pipecg_p_space, n)
            p_prev = p_cur
            u_prev = device_pipecg_p_space + 1
            alpha(1) = alpha(p_cur)
            beta(1) = beta(p_cur)
            p_cur = 1
         else
            u_prev = p_cur
            p_prev = p_cur
            p_cur = p_cur + 1
         end if
      end do
 
      if ( p_cur .ne. 1) then
         call device_memcpy(alpha, alpha_d, p_cur, host_to_device, sync=.false.)
         call device_memcpy(beta, beta_d, p_cur, host_to_device, sync=.false.)
         call device_cg_update_xp(x%x_d, p_d, u_d_d, alpha_d, beta_d, p_cur, &
              device_pipecg_p_space, n)
      end if
      call this%monitor_stop()
      ksp_results%res_final = rnorm
      ksp_results%iter = iter
      ksp_results%converged = this%is_converged(iter, rnorm)
 
    end associate
 
  end function pipecg_device_solve
 
  function pipecg_device_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
       n, coef, blstx, blsty, blstz, gs_h, niter) result(ksp_results)
    class(pipecg_device_t), intent(inout) :: this
    class(ax_t), intent(in) :: ax
    type(field_t), intent(inout) :: x
    type(field_t), intent(inout) :: y
    type(field_t), intent(inout) :: z
    integer, intent(in) :: n
    real(kind=rp), dimension(n), intent(in) :: fx
    real(kind=rp), dimension(n), intent(in) :: fy
    real(kind=rp), dimension(n), intent(in) :: fz
    type(coef_t), intent(inout) :: coef
    type(bc_list_t), intent(inout) :: blstx
    type(bc_list_t), intent(inout) :: blsty
    type(bc_list_t), intent(inout) :: blstz
    type(gs_t), intent(inout) :: gs_h
    type(ksp_monitor_t), dimension(3) :: ksp_results
    integer, optional, intent(in) :: niter
 
    ksp_results(1) =  this%solve(ax, x, fx, n, coef, blstx, gs_h, niter)
    ksp_results(2) =  this%solve(ax, y, fy, n, coef, blsty, gs_h, niter)
    ksp_results(3) =  this%solve(ax, z, fz, n, coef, blstz, gs_h, niter)
 
  end function pipecg_device_solve_coupled
 
end module pipecg_device