hecmw_precond.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_precond
  use hecmw_util
  use hecmw_precond_ssor
  use hecmw_precond_diag
  use hecmw_precond_ml
  use hecmw_precond_bilu
  use hecmw_precond_sainv
  use hecmw_precond_rif
  use hecmw_precond_saamg
  use hecmw_matrix_misc
  use hecmw_solver_las
  implicit none
 
  private
  public :: hecmw_precond_setup
  public :: hecmw_precond_clear
  public :: hecmw_precond_apply
  public :: hecmw_precond_clear_timer
  public :: hecmw_precond_get_timer
 
  real(kind=kreal) :: time_precond = 0.d0
 
contains
 
  subroutine hecmw_precond_setup(hecMAT, hecMESH, sym)
    implicit none
    type (hecmwst_matrix), intent(inout) :: hecmat
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    integer(kind=kint) :: sym
 
    if (hecmw_mat_get_iterpremax( hecmat ).le.0) return
 
    select case(hecmw_mat_get_precond( hecmat ))
      case(1,2)
        call hecmw_precond_ssor_setup(hecmat)
      case(3)
        call hecmw_precond_diag_setup(hecmat)
      case(5)
        call hecmw_precond_ml_setup(hecmat, hecmesh, sym)
      case(10,11,12)
        call hecmw_precond_bilu_setup(hecmat)
      case(20)
        call hecmw_precond_sainv_setup(hecmat)
      case(21)
        call hecmw_precond_rif_setup(hecmat)
      case(22)
        call hecmw_precond_saamg_setup(hecmat, hecmesh, sym)
      case default
        write (*,'(/a )')'#### HEC-MW-SOLVER-E-1001'
        write (*,'( a/)')'    inconsistent solver/preconditioning'
        call hecmw_abort( hecmw_comm_get_comm())
    end select
  end subroutine hecmw_precond_setup
 
  subroutine hecmw_precond_clear(hecMAT)
    implicit none
    type (hecmwst_matrix), intent(inout) :: hecmat
 
    if (hecmw_mat_get_iterpremax( hecmat ).le.0) return
 
    select case(hecmw_mat_get_precond( hecmat ))
      case(1,2)
        call hecmw_precond_ssor_clear(hecmat)
      case(3)
        call hecmw_precond_diag_clear(hecmat%NDOF)
      case(5)
        call hecmw_precond_ml_clear(hecmat%NDOF)
      case(10:12)
        call hecmw_precond_bilu_clear(hecmat%NDOF)
      case(20)
        call hecmw_precond_sainv_clear(hecmat%NDOF)
      case(21)
        call hecmw_precond_rif_clear(hecmat%NDOF)
      case(22)
        call hecmw_precond_saamg_clear(hecmat%NDOF)
      case default
    end select
 
  end subroutine hecmw_precond_clear
 
  subroutine hecmw_precond_apply(hecMESH, hecMAT, R, Z, ZP, COMMtime)
    implicit none
    type (hecmwst_local_mesh), intent(inout) :: hecmesh
    type (hecmwst_matrix), intent(inout)     :: hecmat
#ifdef HECMW_ARCH_FX64
    real(kind=kreal), intent(in)    :: r(hecmat%N*hecmat%NDOF)
    real(kind=kreal), intent(inout) :: z(hecmat%NP*hecmat%NDOF)
    real(kind=kreal), intent(inout) :: zp(hecmat%NP*hecmat%NDOF)
#else
    real(kind=kreal), intent(inout) :: r(:)
    real(kind=kreal), intent(inout) :: z(:), zp(:)
#endif
    real(kind=kreal), intent(inout) :: commtime
    integer(kind=kint ) :: i, n, np, nndof, npndof
    integer(kind=kint) :: iterpremax, iterpre
    real(kind=kreal) :: start_time, end_time
 
    start_time = hecmw_wtime()
 
    n = hecmat%N
    np = hecmat%NP
    nndof = n * hecmat%NDOF
    npndof = np * hecmat%NDOF
 
    if (hecmw_mat_get_iterpremax( hecmat ).le.0) then
      !$acc kernels
      !$acc loop independent
      do i= 1, nndof
        z(i)= r(i)
      enddo
      !$acc end kernels
      return
    endif
 
    !C {z}= [Minv]{r}
    !$acc kernels
    !$acc loop independent
    do i= 1, nndof
      zp(i)= r(i)
    enddo
    !$acc end kernels
 
    !$acc kernels
    !$acc loop independent
    do i= nndof+1, npndof
      zp(i) = 0.d0
    enddo
    !$acc end kernels
 
    !$acc kernels
    !$acc loop independent
    do i= 1, npndof
      z(i)= 0.d0
    enddo
    !$acc end kernels
 
    iterpremax = hecmw_mat_get_iterpremax( hecmat )
    do iterpre= 1, iterpremax
 
      select case(hecmw_mat_get_precond( hecmat ))
        case(1,2)
          call hecmw_precond_ssor_apply(zp,hecmat%NDOF)
        case(3)
          call hecmw_precond_diag_apply(zp,hecmat%NDOF)
        case(5)
          call hecmw_precond_ml_apply(zp,hecmat%NDOF)
        case(10:12)
          call hecmw_precond_bilu_apply(zp,hecmat%NDOF)
        case(20)
          call hecmw_precond_sainv_apply(r,zp,hecmat%NDOF)
        case(21)
          call hecmw_precond_rif_apply(zp,hecmat%NDOF)
        case(22)
          call hecmw_precond_saamg_apply(zp,hecmat%NDOF)
        case default
      end select
 
      !C-- additive Schwartz
      !$acc kernels
      !$acc loop independent
      do i= 1, hecmat%N * hecmat%NDOF
        z(i)= z(i) + zp(i)
      enddo
      !$acc end kernels
      if (iterpre.eq.iterpremax) exit
 
      !C--    {ZP} = {R} - [A] {Z}
      call hecmw_matresid (hecmesh, hecmat, z, r, zp, commtime)
    enddo
 
    end_time = hecmw_wtime()
    time_precond = time_precond + end_time - start_time
  end subroutine hecmw_precond_apply
 
  subroutine hecmw_precond_clear_timer
    implicit none
    time_precond = 0.d0
  end subroutine hecmw_precond_clear_timer
 
  function hecmw_precond_get_timer()
    implicit none
    real(kind=kreal) :: hecmw_precond_get_timer
    hecmw_precond_get_timer = time_precond
  end function hecmw_precond_get_timer
 
end module hecmw_precond