hecmw_solver_GPBiCG.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*** module hecmw_solver_GPBiCG
!C***
!
module hecmw_solver_gpbicg
 
  private
  public :: hecmw_solve_gpbicg
 
contains
  !C
  !C*** hecmw_solve_GPBiCG
  !C
  subroutine hecmw_solve_gpbicg( 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
 
    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
 
    real   (kind=kreal), pointer ::  b(:)
    real   (kind=kreal), pointer ::  x(:)
 
    integer(kind=kint ) :: iterlog, timelog
 
    real(kind=kreal), dimension(:,:),  allocatable       ::  ww
 
    real(kind=kreal), dimension(2) :: rr
 
    integer(kind=kint ) :: maxit
    real   (kind=kreal) :: tol
    integer(kind=kint ) :: i,j
    real   (kind=kreal) :: s_time,s1_time,e_time,e1_time
    real   (kind=kreal) :: bnrm2
    real   (kind=kreal) :: rho,rho1,beta,alpha,dnrm2
    real   (kind=kreal) :: qsi,eta,coef1
    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 :: t= 3
    integer(kind=kint), parameter ::tt= 4
    integer(kind=kint), parameter ::t0= 5
    integer(kind=kint), parameter :: p= 6
    integer(kind=kint), parameter ::pt= 7
    integer(kind=kint), parameter :: u= 8
    integer(kind=kint), parameter ::w1= 9
    integer(kind=kint), parameter :: y=10
    integer(kind=kint), parameter :: z=11
    integer(kind=kint), parameter ::wk=12
    integer(kind=kint), parameter ::w2=13
    integer(kind=kint), parameter ::zq=14
 
    integer(kind=kint), parameter :: n_iter_recompute_r= 20
 
    call hecmw_barrier(hecmesh)
    s_time= hecmw_wtime()
    !C
    !C-- INIT.
    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 )
 
    error= 0
    beta = 0.0d0
 
    allocate (ww(ndof*np,14))
    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)
 
    !C===
    !C +----------------------+
    !C | SETUP PRECONDITIONER |
    !C +----------------------+
    !C===
    call hecmw_precond_setup(hecmat, hecmesh, 0)
 
    !C
    !C +----------------------+
    !C | {r}= {b} - [A]{xini} |
    !C +----------------------+
    !C===
    call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
 
    call hecmw_copy_r(nndof, ww(:,r), ww(:,rt))
    !C==
    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
 
    call hecmw_innerproduct_r(hecmesh, ndof, ww(:,rt), ww(:,r), rho, tcomm)
 
    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
    !C===
 
    !C
    !C*************************************************************** ITERATIVE PROC.
    !C
    call hecmw_barrier(hecmesh)
    s1_time= hecmw_wtime()
    do iter= 1, maxit
      !C
      !C +----------------+
      !C | {r}= [Minv]{r} |
      !C +----------------+
      !C===
      call hecmw_copy_r(nndof, ww(:,r), ww(:,wk))
 
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,wk), ww(:,r), ww(:,zq), tcomm)
      !C===
 
      !C
      !C +----------------------------------+
      !C | {p} = {r} + BETA * ( {p} - {u} ) |
      !C +----------------------------------+
      !C===
      if (iter.gt.1) then
        call hecmw_axpy_r(nndof, -1.0d0, ww(:,u), 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 +--------------------------------+
      !C | ALPHA= {r_tld}{r}/{r_tld} A{p} |
      !C +--------------------------------+
      !C===
 
      !C
      !C-- calc. {p_tld}= A{p}
      call hecmw_matvec(hecmesh, hecmat, ww(:,p), ww(:,pt), tcomm)
 
      !C
      !C-- calc. ALPHA
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,rt), ww(:,pt), rho1, tcomm)
 
      alpha= rho / rho1
      !C===
 
      !C
      !C +------------------------------------------+
      !C | {y}= {t} - {r} - ALPHA{w} + ALPHA{p_tld} |
      !C | {t}= {r}                  - ALPHA{p_tld} |
      !C +------------------------------------------+
      !C===
      call hecmw_axpyz_r(nndof, -1.0d0, ww(:,w1), ww(:,pt), ww(:,y))
      call hecmw_xpay_r (nndof,  alpha, ww(:, t), ww(:, y))
      call hecmw_axpy_r (nndof, -1.0d0, ww(:,wk), ww(:, y))
      call hecmw_axpyz_r(nndof, -alpha, ww(:,pt), ww(:,wk), ww(:,t))
      !C===
 
      !C
      !C +-----------------------+
      !C | {t_tld}= [A][Minv]{t} |
      !C +-----------------------+
      !C===
 
      !C
      !C-- calc. {t_tld} and {t0} by [M] inversion
      !C         {W2}   = [Minv]{p_tld}
      !C
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,t), ww(:,tt), ww(:,zq), tcomm)
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,t0), ww(:,w2), ww(:,zq), tcomm)
      call hecmw_copy_r(nndof, ww(:,w2), ww(:,t0))
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,pt), ww(:,w2), ww(:,zq), tcomm)
 
      !C===
 
      !C
      !C-- calc. [A]{t_tld}
      call hecmw_matvec(hecmesh, hecmat, ww(:,tt), ww(:,wk), tcomm)
 
      call hecmw_copy_r(nndof, ww(:,wk), ww(:,tt))
      !C===
 
      !C
      !C +-------------------+
      !C | calc. QSI and ETA |
      !C +-------------------+
      !C===
      !call pol_coef(iter, WW, T, TT, Y, QSI, ETA)
      !call pol_coef_vanilla(iter, WW, T, TT, Y, QSI, ETA)
      call pol_coef_vanilla2(iter, ww, t, tt, y, qsi, eta)
      !C===
 
      !C
      !C +----------------------------------------------------------+
      !C | {u} = QSI [Minv]{pt} + ETA([Minv]{t0}-[Minv]{r}+BETA*{u} |
      !C | {z} = QSI [Minv]{r}  + ETA*{z} - ALPHA*{u}               |
      !C +----------------------------------------------------------+
      !C===
 
      !C
      !C-- compute. {u},{z}
 
      if (iter.gt.1) then
        call hecmw_xpay_r (nndof,   beta, ww(:,t0), ww(:,u))
        call hecmw_axpy_r (nndof, -1.0d0, ww(:, r), ww(:,u))
        call hecmw_axpby_r(nndof, qsi, eta, ww(:,w2), ww(:,u))
        call hecmw_axpby_r(nndof, -alpha, eta, ww(:, u), ww(:,z))
        call hecmw_axpy_r (nndof, qsi, ww(:,r), ww(:,z))
      else
        call hecmw_axpyz_r(nndof, -1.0d0, ww(:,r), ww(:,t0), ww(:,u))
        call hecmw_axpby_r(nndof, qsi, eta, ww(:,w2), ww(:,u))
        call hecmw_axpby_r(nndof, -alpha, eta, ww(:, u), ww(:,z))
        call hecmw_axpy_r (nndof, qsi, ww(:,r), ww(:,z))
      endif
      !C===
 
      !C
      !C +--------------------+
      !C | update {x},{r},{w} |
      !C +--------------------+
      !C===
      call hecmw_axpy_r(nndof, alpha, ww(:,p), x)
      call hecmw_axpy_r(nndof, 1.0d0, ww(:,z), x)
      call hecmw_copy_r(nndof, ww(:,t), ww(:,t0))
      !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, -eta, ww(:,y), ww(:,t), ww(:,r))
        call hecmw_axpy_r (nndof, -qsi, ww(:,tt), ww(:,r))
      endif
 
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,r), rr(1))
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,rt), rr(2))
      s_time= hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, rr, 2, hecmw_sum)
      e_time= hecmw_wtime()
      tcomm =  tcomm + e_time - s_time
      dnrm2 = rr(1)
      coef1 = rr(2)
 
      beta = alpha*coef1 / (qsi*rho)
      call hecmw_axpyz_r(nndof, beta, ww(:,pt), ww(:,tt), ww(:,w1))
 
      resid= dsqrt(dnrm2/bnrm2)
      rho  = coef1
 
      !C##### ITERATION HISTORY
      if (my_rank.eq.0 .and. iterlog.eq.1)                              &
        &    write (*, 1000) iter, resid
      1000   format (i5, 1pe16.6)
      !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
      !C===
    enddo
 
    call hecmw_solver_scaling_bk(hecmat)
    !C
    !C-- INTERFACE data EXCHANGE
 
    s_time = hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
    e_time = hecmw_wtime()
    tcomm = tcomm + e_time - s_time
 
    deallocate (ww)
    !call hecmw_precond_clear(hecMAT)
 
    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
 
  contains
 
    !C
    !C*** pol_coef : computes QSI and ETA in original GPBiCG way
    !C
    subroutine pol_coef(iter, WW, T, TT, Y, QSI, ETA)
      implicit none
      integer(kind=kint), intent(in) :: iter
      real(kind=kreal), intent(in) :: ww(:,:)
      integer(kind=kint), intent(in) :: T, TT, Y
      real(kind=kreal), intent(out) :: qsi, eta
 
      real(kind=kreal), dimension(5) :: cg
      real(kind=kreal), dimension(2) :: eq
      real(kind=kreal) :: delta
 
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,y ), cg(1)) ! myu1
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,tt), ww(:,t ), cg(2)) ! omega2
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,t ), cg(3)) ! omega1
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,tt), ww(:,y ), cg(4)) ! nyu
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,tt), ww(:,tt), cg(5)) ! myu2
      s_time= hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, cg, 5, hecmw_sum)
      e_time= hecmw_wtime()
      tcomm =  tcomm + e_time - s_time
 
      if (iter.eq.1) then
        eq(1)= cg(2)/cg(5) ! omega2 / myu2
        eq(2)= 0.d0
      else
        delta= (cg(5)*cg(1)-cg(4)*cg(4))         ! myu1*myu2 - nyu^2
        eq(1)= (cg(1)*cg(2)-cg(3)*cg(4)) / delta ! (myu1*omega2-nyu*omega1)/delta
        eq(2)= (cg(5)*cg(3)-cg(4)*cg(2)) / delta ! (myu2*omega1-nyu*omega2)/delta
      endif
 
      qsi= eq(1)
      eta= eq(2)
    end subroutine pol_coef
 
    !C
    !C*** pol_coef_vanilla : computes QSI and ETA with vanilla strategy
    !C      see Fujino, Abe, Sugihara and Nakashima(2013) ISBN978-4-621-08741-1
    !C
    subroutine pol_coef_vanilla(iter, WW, T, TT, Y, QSI, ETA)
      implicit none
      integer(kind=kint), intent(in) :: iter
      real(kind=kreal), intent(inout) :: ww(:,:)
      integer(kind=kint), intent(in) :: T, TT, Y
      real(kind=kreal), intent(out) :: qsi, eta
 
      real(kind=kreal), dimension(3) :: cg
      real(kind=kreal) :: gamma1, gamma2
      real(kind=kreal) :: c, c_abs
 
      real(kind=kreal), parameter :: omega = 0.707106781d0
 
      if (iter.eq.1) then
        gamma1 = 0.d0
        gamma2 = 0.d0
      else
        call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,y ), cg(1)) ! myu
        call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,tt), cg(2)) ! nyu
        call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,t ), cg(3)) ! omega
        s_time= hecmw_wtime()
        call hecmw_allreduce_r(hecmesh, cg, 3, hecmw_sum)
        e_time= hecmw_wtime()
        tcomm =  tcomm + e_time - s_time
        gamma1 = cg(3)/cg(1) ! omega / myu
        gamma2 = cg(2)/cg(1) ! nyu / myu
        !!! COMMENTED OUT because no convergence obtained with following updates.
        !         do j= 1, NNDOF
        !           WW(j,T )= WW(j,T ) - gamma1*WW(j,Y)
        !           WW(j,TT)= WW(j,TT) - gamma2*WW(j,Y)
        !         enddo
      endif
 
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t ), ww(:,t ), cg(1)) ! |r|^2
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,tt), ww(:,tt), cg(2)) ! |s|^2
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t ), ww(:,tt), cg(3)) ! r.s
      s_time= hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, cg, 3, hecmw_sum)
      e_time= hecmw_wtime()
      tcomm =  tcomm + e_time - s_time
 
      c = cg(3) / dsqrt(cg(1)*cg(2))
      c_abs = dabs(c)
      if (c_abs > omega) then
        qsi = c * dsqrt(cg(1)/cg(2))
      else
        ! QSI = (c / c_abs) * OMEGA * dsqrt(CG(1)/CG(2))
        if (c >= 0.d0) then
          qsi = omega * dsqrt(cg(1)/cg(2))
        else
          qsi = -omega * dsqrt(cg(1)/cg(2))
        endif
      endif
      eta = gamma1 - qsi*gamma2
    end subroutine pol_coef_vanilla
 
    !C
    !C*** pol_coef_vanilla2 : optimized version of pol_coef_vanilla
    !C
    subroutine pol_coef_vanilla2(iter, WW, T, TT, Y, QSI, ETA)
      implicit none
      integer(kind=kint), intent(in) :: iter
      real(kind=kreal), intent(inout) :: ww(:,:)
      integer(kind=kint), intent(in) :: T, TT, Y
      real(kind=kreal), intent(out) :: qsi, eta
 
      real(kind=kreal), dimension(6) :: cg
      real(kind=kreal) :: gamma1, gamma2
      real(kind=kreal) :: c, c_abs
 
      real(kind=kreal), parameter :: omega = 0.707106781d0
 
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t ), ww(:,t ), cg(1)) ! |r|^2
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,tt), ww(:,tt), cg(2)) ! |s|^2
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t ), ww(:,tt), cg(3)) ! r.s
 
      if (iter.gt.1) then
        call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,y ), cg(4)) ! myu
        call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,tt), cg(5)) ! nyu
        call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,y ), ww(:,t ), cg(6)) ! omega
        s_time= hecmw_wtime()
        call hecmw_allreduce_r(hecmesh, cg, 6, hecmw_sum)
        e_time= hecmw_wtime()
        tcomm =  tcomm + e_time - s_time
        gamma1 = cg(6)/cg(4) ! omega / myu
        gamma2 = cg(5)/cg(4) ! nyu / myu
      else
        s_time= hecmw_wtime()
        call hecmw_allreduce_r(hecmesh, cg, 3, hecmw_sum)
        e_time= hecmw_wtime()
        tcomm =  tcomm + e_time - s_time
        gamma1 = 0.d0
        gamma2 = 0.d0
      endif
 
      c = cg(3) / dsqrt(cg(1)*cg(2))
      c_abs = dabs(c)
      if (c_abs > omega) then
        qsi = c * dsqrt(cg(1)/cg(2))
      else
        if (c >= 0.d0) then
          qsi = omega * dsqrt(cg(1)/cg(2))
        else
          qsi = -omega * dsqrt(cg(1)/cg(2))
        endif
      endif
      eta = gamma1 - qsi*gamma2
    end subroutine pol_coef_vanilla2
 
  end subroutine  hecmw_solve_gpbicg
 
end module     hecmw_solver_gpbicg