hecmw_solver_BiCGSTAB.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
!C
!C***
!C*** module hecmw_solver_BiCGSTAB
!C***
!C
module hecmw_solver_bicgstab
contains
  !C
  !C*** BiCGSTAB_3
  !C
  subroutine hecmw_solve_bicgstab( hecMESH,  hecMAT, ITER, RESID, error, &
      &                                    Tset, Tsol, Tcomm )
 
    use hecmw_util
    use m_hecmw_solve_error
    use m_hecmw_comm_f
    use hecmw_matrix_misc
    use hecmw_solver_misc
    use hecmw_solver_las
    use hecmw_solver_scaling
    use hecmw_precond
    use hecmw_jad_type
 
    implicit none
 
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix) :: hecMAT
    integer(kind=kint ), intent(inout):: ITER, error
    real   (kind=kreal), intent(inout):: resid, tset, tsol, tcomm
 
    integer(kind=kint ) :: N, NP, NDOF, NNDOF
    integer(kind=kint ) :: my_rank
    integer(kind=kint ) :: ITERlog, TIMElog
    real(kind=kreal), pointer :: b(:), x(:)
 
    real(kind=kreal), dimension(:,:), allocatable :: ww
    real(kind=kreal), dimension(2) :: cg
 
    integer(kind=kint ) :: MAXIT
 
    ! local variables
    real   (kind=kreal):: tol
    integer(kind=kint )::i
    real   (kind=kreal)::s_time,s1_time,e_time,e1_time, start_time, end_time
    real   (kind=kreal)::bnrm2,c2
    real   (kind=kreal)::rho,rho1,beta,alpha,dnrm2
    real   (kind=kreal)::omega
    real   (kind=kreal)::t_max,t_min,t_avg,t_sd
 
    integer(kind=kint), parameter :: R = 1
    integer(kind=kint), parameter :: RT= 2
    integer(kind=kint), parameter :: P = 3
    integer(kind=kint), parameter :: PT= 4
    integer(kind=kint), parameter :: S = 5
    integer(kind=kint), parameter :: ST= 1
    integer(kind=kint), parameter :: T = 6
    integer(kind=kint), parameter :: V = 7
    integer(kind=kint), parameter :: WK= 8
 
    integer(kind=kint) :: N_ITER_RECOMPUTE_R
    integer(kind=kint), parameter :: N_ITER_RECOMPUTE_R_DEFAULT= 100
 
    call hecmw_barrier(hecmesh)
    s_time= hecmw_wtime()
 
    !C===
    !C +-------+
    !C | INIT. |
    !C +-------+
    !C===
    n = hecmat%N
    np = hecmat%NP
    ndof = hecmat%NDOF
    nndof = n * ndof
    my_rank = hecmesh%my_rank
    x => hecmat%X
    b => hecmat%B
 
    iterlog = hecmw_mat_get_iterlog( hecmat )
    timelog = hecmw_mat_get_timelog( hecmat )
    maxit   = hecmw_mat_get_iter( hecmat )
    tol     = hecmw_mat_get_resid( hecmat )
 
    n_iter_recompute_r = hecmw_mat_get_recompute_residual( hecmat )
    if (n_iter_recompute_r == 0) n_iter_recompute_r = n_iter_recompute_r_default
    !C----- negative: never recompute periodically.  ITER stops at MAXIT, so mod() below is never 0
    if (n_iter_recompute_r < 0) n_iter_recompute_r = maxit + 1
 
    error = 0
    rho1 = 0.0d0
    alpha = 0.0d0
    beta = 0.0d0
    omega = 0.0d0
 
    allocate (ww(ndof*np, 8))
    ww = 0.d0
 
    !C
    !C-- SCALING
    call hecmw_solver_scaling_fw(hecmesh, hecmat, tcomm)
 
    !C
    !C-- matrix integration for OpenACC
    !C
    !C  @note:
    !C   Combine hecMAT%AL, D, and AU into a single matrix for GPU execution.
    !C   This is a no-op for CPU builds.
    call hecmw_mat_integrate(hecmat)
 
    if (hecmw_mat_get_usejad(hecmat).ne.0) then
      call hecmw_jad_init(hecmat)
    endif
    call hecmw_matvec_setup(hecmesh, hecmat)
 
    !C===
    !C +----------------------+
    !C | SETUP PRECONDITIONER |
    !C +----------------------+
    !C===
    call hecmw_precond_setup(hecmat, hecmesh, 0)
 
    !C===
    !C +---------------------+
    !C | {r0}= {b} - [A]{x0} |
    !C +---------------------+
    !C===
    call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
 
    !C-- set arbitrary {r_tld}
    call hecmw_copy_r(nndof, ww(:,r), ww(:,rt))
 
    !C-- compute ||{b}||
    call hecmw_innerproduct_r(hecmesh, ndof, b, b, bnrm2, tcomm)
    if (bnrm2.eq.0.d0) then
      iter = 0
      maxit = 0
      resid = 0.d0
      x = 0.d0
    endif
 
    e_time = hecmw_wtime()
    if (timelog.eq.2) then
      call hecmw_time_statistics(hecmesh, e_time - s_time, &
        t_max, t_min, t_avg, t_sd)
      if (hecmesh%my_rank.eq.0) then
        write(*,*) 'Time solver setup'
        write(*,*) '  Max     :',t_max
        write(*,*) '  Min     :',t_min
        write(*,*) '  Avg     :',t_avg
        write(*,*) '  Std Dev :',t_sd
      endif
      tset = t_max
    else
      tset = e_time - s_time
    endif
 
    tcomm = 0.d0
    call hecmw_barrier(hecmesh)
    s1_time = hecmw_wtime()
    !C
    !C*************************************************************** iterative procedures start
    !C
    do iter = 1, maxit
 
      !C===
      !C +-----------------+
      !C | RHO= {r}{r_tld} |
      !C +-----------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,rt), rho, tcomm)
 
      !C===
      !C +----------------------------------------+
      !C | BETA= (RHO/RHO1) * (ALPHA/OMEGA)       |
      !C | {p} = {r} + BETA * ( {p} - OMEGA*{v} ) |
      !C +----------------------------------------+
      !C===
      if ( iter.gt.1 ) then
        beta = (rho/rho1) * (alpha/omega)
        call hecmw_axpy_r(nndof, -omega, ww(:,v), ww(:,p))
        call hecmw_xpay_r(nndof,   beta, ww(:,r), ww(:,p))
      else
        call hecmw_copy_r(nndof, ww(:,r), ww(:,p))
      endif
 
      !C===
      !C +--------------------+
      !C | {p_tld}= [Minv]{p} |
      !C +--------------------+
      !C===
      call hecmw_precond_apply(hecmesh, hecmat, ww(:, p), ww(:, pt), ww(:, wk), tcomm)
 
      !C===
      !C +-------------------------+
      !C | {v}= [A] {p_tld} |
      !C +-------------------------+
      !C===
      call hecmw_matvec(hecmesh, hecmat, ww(:,pt), ww(:,v), tcomm)
 
      !C
      !C-- calc. ALPHA
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,rt), ww(:,v), c2, tcomm)
 
      alpha = rho / c2
 
      !C
      !C-- {s}= {r} - ALPHA*{V}
      call hecmw_axpyz_r(nndof, -alpha, ww(:,v), ww(:,r), ww(:,s))
 
      !C===
      !C +--------------------+
      !C | {s_tld}= [Minv]{s} |
      !C +--------------------+
      !C===
      call hecmw_precond_apply(hecmesh, hecmat, ww(:, s), ww(:, st), ww(:, wk), tcomm)
 
      !C===
      !C +-------------------------+
      !C | {t}= [A] {s_tld} |
      !C +-------------------------+
      !C===
      call hecmw_matvec(hecmesh, hecmat, ww(:,st), ww(:,t), tcomm)
 
      !C===
      !C +----------------------------+
      !C | OMEGA= ({t}{s}) / ({t}{t}) |
      !C +----------------------------+
      !C===
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t), ww(:,s), cg(1))
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t), ww(:,t), cg(2))
      s_time= hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, cg, 2, hecmw_sum)
      e_time= hecmw_wtime()
      tcomm = tcomm + e_time - s_time
 
      omega = cg(1) / cg(2)
 
      !C===
      !C +----------------+
      !C | update {x},{r} |
      !C +----------------+
      !C===
      call hecmw_axpy_r(nndof, alpha, ww(:,pt), x)
      call hecmw_axpy_r(nndof, omega, ww(:,st), x)
      !C
      !C--- recompute R sometimes
      if ( mod(iter,n_iter_recompute_r)==0 ) then
        call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
      else
        call hecmw_axpyz_r(nndof, -omega, ww(:,t), ww(:,s), ww(:,r))
      endif
 
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
 
      resid= dsqrt(dnrm2/bnrm2)
 
      !C##### ITERATION HISTORY
      if (my_rank.eq.0.and.iterlog.eq.1) write (*,'(i7, 1pe16.6)') iter, resid
      !C#####
 
      if ( resid.le.tol   ) then
        if ( mod(iter,n_iter_recompute_r)==0 ) exit
        !C----- recompute R to make sure it is really converged
        call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
        call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
        resid= dsqrt(dnrm2/bnrm2)
        if ( resid.le.tol ) exit
      endif
      if ( iter .eq.maxit ) error = hecmw_solver_error_noconv_maxit
 
      rho1 = rho
 
    enddo
    !C
    !C*************************************************************** iterative procedures end
    !C
 
    call hecmw_solver_scaling_bk(hecmat)
    !C
    !C-- INTERFACE data EXCHANGE
    !C
    start_time = hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
    end_time = hecmw_wtime()
    tcomm = tcomm + end_time - start_time
 
    deallocate (ww)
    !call hecmw_precond_clear(hecMAT)
 
    call hecmw_matvec_teardown(hecmat)
    if (hecmw_mat_get_usejad(hecmat).ne.0) then
      call hecmw_jad_finalize(hecmat)
    endif
 
    e1_time = hecmw_wtime()
    if (timelog.eq.2) then
      call hecmw_time_statistics(hecmesh, e1_time - s1_time, &
        t_max, t_min, t_avg, t_sd)
      if (hecmesh%my_rank.eq.0) then
        write(*,*) 'Time solver iterations'
        write(*,*) '  Max     :',t_max
        write(*,*) '  Min     :',t_min
        write(*,*) '  Avg     :',t_avg
        write(*,*) '  Std Dev :',t_sd
      endif
      tsol = t_max
    else
      tsol = e1_time - s1_time
    endif
 
  end subroutine hecmw_solve_bicgstab
end module     hecmw_solver_bicgstab