hecmw_solver_las_11.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_solver_las_11
 
  private
 
  public :: hecmw_matvec_11
  public :: hecmw_matresid_11
  public :: hecmw_rel_resid_l2_11
 
contains
 
  !C
  !C***
  !C*** hecmw_matvec_11
  !C***
  !C
  subroutine hecmw_matvec_11 (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
    use hecmw_util
    use m_hecmw_comm_f
 
    implicit none
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (hecmwst_matrix), intent(in)     :: hecmat
    real(kind=kreal), intent(in) :: x(:)
    real(kind=kreal), intent(out) :: y(:)
    real(kind=kreal), intent(inout) :: time_ax
    real(kind=kreal), intent(inout), optional :: commtime
 
    real(kind=kreal) :: start_time, end_time
    integer(kind=kint) :: i, j, js, je, in
    real(kind=kreal) :: yv
 
    start_time= hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
    end_time= hecmw_wtime()
    if (present(commtime)) commtime = commtime + end_time - start_time
 
    start_time= hecmw_wtime()
#ifdef _OPENACC
    !$acc kernels
    !$acc loop independent
    do i= 1, hecmat%N
      yv= 0
      js= hecmat%indexA(i) + 1
      je= hecmat%indexA(i+1)
      do j= js, je
        in= hecmat%itemA(j)
        yv= yv + hecmat%A(j) * x(in)
      enddo
      y(i)= yv
    enddo
    !$acc end kernels
#else
    do i= 1, hecmat%N
      yv= hecmat%D(i) * x(i)
      js= hecmat%indexL(i-1) + 1
      je= hecmat%indexL(i  )
      do j= js, je
        in= hecmat%itemL(j)
        yv= yv + hecmat%AL(j) * x(in)
      enddo
      js= hecmat%indexU(i-1) + 1
      je= hecmat%indexU(i  )
      do j= js, je
        in= hecmat%itemU(j)
        yv= yv + hecmat%AU(j) * x(in)
      enddo
      y(i)= yv
    enddo
#endif
    end_time = hecmw_wtime()
    time_ax = time_ax + end_time - start_time
 
  end subroutine hecmw_matvec_11
 
  !C
  !C***
  !C*** hecmw_matresid_11
  !C***
  !C
  subroutine hecmw_matresid_11 (hecMESH, hecMAT, X, B, R, time_Ax, COMMtime)
    use hecmw_util
 
    implicit none
    real(kind=kreal) :: x(:), b(:), r(:)
    type (hecmwst_matrix)     :: hecmat
    type (hecmwst_local_mesh) :: hecmesh
    real(kind=kreal) :: time_ax
    real(kind=kreal), optional :: commtime
 
    integer(kind=kint) :: i
    real(kind=kreal) :: tcomm
 
    tcomm = 0.d0
    call hecmw_matvec_11 (hecmesh, hecmat, x, r, time_ax, tcomm)
    if (present(commtime)) commtime = commtime + tcomm
 
    !$acc kernels
    !$acc loop independent
    do i = 1, hecmat%N
      r(i) = b(i) - r(i)
    enddo
    !$acc end kernels
 
  end subroutine hecmw_matresid_11
 
  !C
  !C***
  !C*** hecmw_rel_resid_L2_11
  !C***
  !C
  function hecmw_rel_resid_l2_11 (hecMESH, hecMAT, time_Ax, COMMtime)
    use hecmw_util
    use hecmw_solver_misc
 
    implicit none
    real(kind=kreal) :: hecmw_rel_resid_l2_11
    type ( hecmwst_local_mesh ), intent(in) :: hecmesh
    type ( hecmwst_matrix     ), intent(in) :: hecmat
    real(kind=kreal) :: time_ax
    real(kind=kreal), optional :: commtime
 
    real(kind=kreal) :: r(hecmat%NDOF*hecmat%NP)
    real(kind=kreal) :: bnorm2, rnorm2
    real(kind=kreal) :: tcomm
 
    tcomm = 0.d0
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, hecmat%B, hecmat%B, bnorm2, tcomm)
    if (bnorm2 == 0.d0) then
      bnorm2 = 1.d0
    endif
    call hecmw_matresid_11(hecmesh, hecmat, hecmat%X, hecmat%B, r, time_ax, tcomm)
    call hecmw_innerproduct_r(hecmesh, hecmat%NDOF, r, r, rnorm2, tcomm)
    if (present(commtime)) commtime = commtime + tcomm
    hecmw_rel_resid_l2_11 = sqrt(rnorm2 / bnorm2)
 
  end function hecmw_rel_resid_l2_11
 
end module hecmw_solver_las_11