hecmw_solver_CR.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_CR
!C***
!
module hecmw_solver_cr
 
  public :: hecmw_solve_cr
 
contains
  !C
  !C*** hecmw_solve_CR
  !C
  subroutine hecmw_solve_cr(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
 
    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, rtar, rtar_old, aptap
    real   (kind=kreal)::alpha, beta, dnrm2, dnrm2_true
    real   (kind=kreal)::resid_rec, resid_true
    real   (kind=kreal)::t_max, t_min, t_avg, t_sd
    logical :: true_checked
 
    ! PCR vector definitions
    integer(kind=kint), parameter :: r = 1    ! {z} = [Minv]{rt}
    integer(kind=kint), parameter :: p = 2    ! search direction {p}
    integer(kind=kint), parameter :: ap = 3   ! [A]{p}
    integer(kind=kint), parameter :: az = 4   ! [A]{z}
    integer(kind=kint), parameter :: wk = 4   ! preconditioner / true-residual workspace; shares AZ column (lifetimes disjoint)
    integer(kind=kint), parameter :: map = 5  ! [Minv][A]{p}
    integer(kind=kint), parameter :: rt = 6   ! recursive residual {rt} = {b} - [A]{x}
 
    integer(kind=kint), parameter :: n_iter_check_true_r = 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 )
 
    error = 0
    rtar_old = 0.0d0
 
    allocate (ww(ndof*np, 6))
    ww = 0.d0
 
    !C
    !C-- SCALING
    call hecmw_solver_scaling_fw(hecmesh, hecmat, tcomm)
 
    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, 1)
 
    !C===
    !C +-----------------------+
    !C | {rt}= {b} - [A]{x}    |
    !C +-----------------------+
    !C===
    call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,rt), 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
    else
      !C-- Check the true initial residual before constructing PCR vectors
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,rt), ww(:,rt), dnrm2, tcomm)
      resid = dsqrt(dnrm2/bnrm2)
      if (resid.le.tol) then
        iter = 0
        maxit = 0
      else
        !C-- {z}= [Minv]{rt}
        call hecmw_precond_apply(hecmesh, hecmat, ww(:,rt), ww(:,r), ww(:,wk), tcomm)
 
        !C-- {p}= {z}
        call hecmw_copy_r(nndof, ww(:,r), ww(:,p))
 
        !C-- {az}= [A]{z}, {ap}= {az}
        call hecmw_matvec(hecmesh, hecmat, ww(:,r), ww(:,az), tcomm)
        call hecmw_copy_r(nndof, ww(:,az), ww(:,ap))
 
        !C-- rTAr_old= {z}^T[A]{z}
        call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,az), rtar_old, tcomm)
        if (rtar_old /= rtar_old) then
          error = hecmw_solver_error_diverge_nan
          iter = 0
          maxit = 0
        elseif (rtar_old <= 0.0d0) then
          error = hecmw_solver_error_diverge_mat
          iter = 0
          maxit = 0
        endif
      endif
    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*************************************************************** Iteration begins
    !C
    do iter = 1, maxit
 
      !C===
      !C +-------------------------+
      !C | {map}= [Minv]{ap}    |
      !C +-------------------------+
      !C===
      ! AZ and WK share one column because [A]{z} is dead while the preconditioner workspace is needed.
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,ap), ww(:,map), ww(:,wk), tcomm)
 
      !C===
      !C +--------------------------------+
      !C | ApTAp= {ap}^T[Minv]{ap}    |
      !C +--------------------------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,ap), ww(:,map), aptap, tcomm)
 
      if (aptap /= aptap) then
        error = hecmw_solver_error_diverge_nan
        exit
      elseif (aptap <= 0.0d0) then
        ! A symmetric positive definite preconditioner must give {ap}^T[Minv]{ap} > 0 for nonzero {ap}.
        error = hecmw_solver_error_diverge_pc
        exit
      endif
      alpha = rtar_old / aptap
 
      !C===
      !C +--------------------------------+
      !C | {x}= {x} + ALPHA*{p}        |
      !C | {rt}= {rt} - ALPHA*{ap}     |
      !C | {z}= {z} - ALPHA*{map}      |
      !C +--------------------------------+
      !C===
      call hecmw_axpy_r(nndof,  alpha, ww(:,p), x)
      call hecmw_axpy_r(nndof, -alpha, ww(:,ap), ww(:,rt))
      call hecmw_axpy_r(nndof, -alpha, ww(:,map), ww(:,r))
 
      !C===
      !C +------------------------------------+
      !C | RESID_REC= ||{rt}|| / ||{b}||     |
      !C +------------------------------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,rt), ww(:,rt), dnrm2, tcomm)
      resid_rec = dsqrt(dnrm2/bnrm2)
 
      !C===
      !C +------------------------------------------------+
      !C | Explicit residual RESID_TRUE for monitoring   |
      !C +------------------------------------------------+
      !C===
      ! Compute the explicit residual only for convergence monitoring.
      ! Do not replace the recursive PCR residual or restart the search directions.
      ! It is checked when the recursive residual reaches the tolerance, on a
      ! periodic monitoring iteration, or at the maximum iteration count.
      true_checked = .false.
      if (resid_rec.le.tol .or. mod(iter,n_iter_check_true_r)==0 .or. iter.eq.maxit) then
        ! WK reuses the AZ column, which is free until [A]{z} is formed below.
        call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,wk), tcomm)
        call hecmw_innerproduct_r(hecmesh, ndof, ww(:,wk), ww(:,wk), dnrm2_true, tcomm)
        resid_true = dsqrt(dnrm2_true/bnrm2)
        true_checked = .true.
      endif
 
      ! Return/display value: the explicit residual on a monitoring iteration,
      ! the recursive residual otherwise.
      if (true_checked) then
        resid = resid_true
      else
        resid = resid_rec
      endif
 
      !C##### Iteration history output
      if (my_rank.eq.0.and.iterlog.eq.1) write (*,'(i7, 1pe16.6)') iter, resid
      !C#####
 
      if (true_checked) then
        if (resid_true.le.tol) exit
        if (iter.eq.maxit) then
          error = hecmw_solver_error_noconv_maxit
          exit
        endif
      endif
 
      !C===
      !C +--------------------+
      !C | {az}= [A]{z}      |
      !C +--------------------+
      !C===
      call hecmw_matvec(hecmesh, hecmat, ww(:,r), ww(:,az), tcomm)
 
      !C===
      !C +-----------------------------+
      !C | rTAr= {z}^T[A]{z}          |
      !C | BETA= rTAr / rTAr_old      |
      !C +-----------------------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,az), rtar, tcomm)
 
      if (rtar /= rtar) then
        error = hecmw_solver_error_diverge_nan
        exit
      elseif (rtar <= 0.0d0) then
        error = hecmw_solver_error_diverge_mat
        exit
      endif
 
      beta = rtar / rtar_old
      rtar_old = rtar
 
      !C===
      !C +------------------------------+
      !C | {p}= {z} + BETA*{p}        |
      !C | {ap}= {az} + BETA*{ap}     |
      !C +------------------------------+
      !C===
      call hecmw_xpay_r(nndof, beta, ww(:,r), ww(:,p))
      call hecmw_xpay_r(nndof, beta, ww(:,az), ww(:,ap))
 
    enddo
    if (maxit.eq.0) iter = 0
    !C
    !C*************************************************************** Iteration ends
    !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_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_cr
end module hecmw_solver_cr