hecmw_solver_CG.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_CG
!C***
!C
module hecmw_solver_cg
 
  public :: hecmw_solve_cg
 
contains
  !C
  !C*** CG_nn
  !C
  subroutine hecmw_solve_cg( 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
    use hecmw_estimate_condition
 
    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
 
    integer(kind=kint), parameter ::  r= 1
    integer(kind=kint), parameter ::  z= 2
    integer(kind=kint), parameter ::  q= 2
    integer(kind=kint), parameter ::  p= 3
    integer(kind=kint), parameter :: wk= 4
 
    integer(kind=kint ) :: maxit
 
    ! local variables
    real   (kind=kreal) :: tol
    real   (kind=kreal)::s_time,s1_time,e_time,e1_time, start_time, end_time
    real   (kind=kreal)::bnrm2
    real   (kind=kreal)::rho,rho1,beta,c1,alpha,dnrm2
    real   (kind=kreal)::t_max,t_min,t_avg,t_sd
    integer(kind=kint) ::estcond
    real   (kind=kreal), allocatable::d(:),e(:)
    real   (kind=kreal)::alpha1
    integer(kind=kint) :: n_indef_precond
 
    integer(kind=kint) :: n_iter_recompute_r
    integer(kind=kint), parameter :: n_iter_recompute_r_default= 50
 
    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 )
    estcond = hecmw_mat_get_estcond( 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
    n_indef_precond = 0
    rho1 = 0.0d0
    alpha1 = 0.0d0
    beta = 0.0d0
 
    allocate (ww(ndof*np, 4))
    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)
 
    if (estcond /= 0 .and. hecmesh%my_rank == 0) then
      allocate(d(maxit),e(maxit-1))
    endif
    !C===
    !C +----------------------+
    !C | SETUP PRECONDITIONER |
    !C +----------------------+
    !C===
    call hecmw_precond_setup(hecmat, hecmesh, 1)
 
    !C===
    !C +---------------------+
    !C | {r0}= {b} - [A]{x0} |
    !C +---------------------+
    !C===
    call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
 
    !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************************************************* Conjugate Gradient Iteration start
    !C
    do iter = 1, maxit
 
      !C===
      !C +----------------+
      !C | {z}= [Minv]{r} |
      !C +----------------+
      !C===
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,z), ww(:,wk), tcomm)
 
 
      !C===
      !C +---------------+
      !C | {RHO}= {r}{z} |
      !C +---------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,z), rho, tcomm)
      ! if RHO is NaN or Inf then no converge
      if (rho == 0.d0) then
        ! converged due to RHO==0
        exit
      elseif (iter > 1 .and. rho*rho1 <= 0) then
        n_indef_precond = n_indef_precond + 1
        if (n_indef_precond >= 10) then
          ! diverged due to indefinite preconditioner
          error = hecmw_solver_error_diverge_pc
          exit
        endif
      elseif (rho /= rho) then  ! RHO is NaN
        error = hecmw_solver_error_diverge_nan
        exit
      else
        n_indef_precond = 0
      endif
 
      !C===
      !C +-----------------------------+
      !C | {p} = {z} if      ITER=1    |
      !C | BETA= RHO / RHO1  otherwise |
      !C +-----------------------------+
      !C===
      if ( iter.eq.1 ) then
        call hecmw_copy_r(nndof, ww(:,z), ww(:,p))
      else
        beta = rho / rho1
        call hecmw_xpay_r(nndof, beta, ww(:,z), ww(:,p))
      endif
 
      !C===
      !C +--------------+
      !C | {q}= [A] {p} |
      !C +--------------+
      !C===
      call hecmw_matvec(hecmesh, hecmat, ww(:,p), ww(:,q), tcomm)
 
      !C===
      !C +---------------------+
      !C | ALPHA= RHO / {p}{q} |
      !C +---------------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,p), ww(:,q), c1, tcomm)
      ! if C1 is NaN or Inf, no converge
      if (c1 <= 0) then
        ! diverged due to indefinite or negative definite matrix
        error = hecmw_solver_error_diverge_mat
        exit
      elseif (c1 /= c1) then  ! C1 is NaN
        error = hecmw_solver_error_diverge_nan
        exit
      endif
 
      alpha= rho / c1
 
      !C===
      !C +----------------------+
      !C | {x}= {x} + ALPHA*{p} |
      !C | {r}= {r} - ALPHA*{q} |
      !C +----------------------+
      !C===
      call hecmw_axpy_r(nndof, alpha, ww(:,p), x)
 
      if ( mod(iter,n_iter_recompute_r)==0 ) then
        call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
      else
        call hecmw_axpy_r(nndof, -alpha, ww(:,q), 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 (estcond /= 0 .and. hecmesh%my_rank == 0) then
        if (iter == 1) then
          d(1) = 1.d0 / alpha
        else
          d(iter) = 1.d0 / alpha + beta / alpha1
          e(iter-1) = sqrt(beta) / alpha1
        endif
        alpha1 = alpha
        if (mod(iter,estcond) == 0) call hecmw_estimate_condition_cg(iter, d, e)
      endif
 
      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************************************************* Conjugate Gradient Iteration 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
 
    if (estcond /= 0 .and. error == 0 .and. hecmesh%my_rank == 0) then
      call hecmw_estimate_condition_cg(iter, d, e)
      deallocate(d, e)
    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_cg
 
end module     hecmw_solver_cg