hecmw_solver_PipeCG.f90

Go to the documentation of this file.

!C
!C*** module hecmw_solver_PipeCG
!C
module hecmw_solver_pipecg
 
  public :: hecmw_solve_pipecg
 
contains
  !C
  !C*** Pipelined CG
  !C
  subroutine hecmw_solve_pipecg( 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 :: u  = 2
    integer(kind=kint), parameter :: v  = 3
    integer(kind=kint), parameter :: q  = 4
    integer(kind=kint), parameter :: p  = 5
    integer(kind=kint), parameter :: z  = 6
    integer(kind=kint), parameter :: l  = 7
    integer(kind=kint), parameter :: m  = 8
    integer(kind=kint), parameter :: s  = 9
    integer(kind=kint), parameter :: wk = 10
 
    integer(kind=kint) :: maxit
 
    real(kind=kreal) :: tol
    integer(kind=kint) :: i
    real(kind=kreal) :: s_time, s1_time, e_time, e1_time
    real(kind=kreal) :: start_time, end_time
    real(kind=kreal) :: bnrm2, dnrm2
    real(kind=kreal) :: alpha, alpha1, beta
    real(kind=kreal) :: gamma, gamma1, delta, denom
    real(kind=kreal) :: cg(3)
    real(kind=kreal) :: t_max, t_min, t_avg, t_sd
    integer(kind=kint) :: estcond
    real(kind=kreal), allocatable :: d(:), e(:)
    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()
 
    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.  i stops at MAXIT, so mod() below is never 0
    if (n_iter_recompute_r < 0) n_iter_recompute_r = maxit + 1
 
    error = 0
    iter = 0
    n_indef_precond = 0
    gamma1 = 1.0d0
    alpha1 = 1.0d0
 
    allocate(ww(ndof*np, 10))
    ww = 0.0d0
 
    call hecmw_solver_scaling_fw(hecmesh, hecmat, tcomm)
    call hecmw_mat_integrate(hecmat)
 
    if (hecmw_mat_get_usejad(hecmat).ne.0) then
      call hecmw_jad_init(hecmat)
    endif
 
    if (estcond /= 0 .and. hecmesh%my_rank == 0) then
      allocate(d(maxit), e(maxit-1))
    endif
 
    call hecmw_precond_setup(hecmat, hecmesh, 1)
 
    call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
    call hecmw_innerproduct_r(hecmesh, ndof, b, b, bnrm2, tcomm)
    if (bnrm2.eq.0.0d0) then
      maxit = 0
      resid = 0.0d0
      x = 0.0d0
    else
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,u), ww(:,wk), tcomm)
      call hecmw_matvec(hecmesh, hecmat, ww(:,u), ww(:,v), tcomm)
    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.0d0
    call hecmw_barrier(hecmesh)
    s1_time = hecmw_wtime()
 
    do i = 1, maxit
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,u), cg(1))
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,v), ww(:,u), cg(2))
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,r), cg(3))
      start_time = hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, cg, 3, hecmw_sum)
      end_time = hecmw_wtime()
      tcomm = tcomm + end_time-start_time
 
      gamma = cg(1)
      delta = cg(2)
      dnrm2 = cg(3)
      resid = dsqrt(dnrm2/bnrm2)
 
      if (i.eq.1 .and. resid.le.tol) exit
 
      if (i.gt.1) then
        iter = i-1
        if (my_rank.eq.0 .and. iterlog.eq.1) write(*,'(i7, 1pe16.6)') iter, resid
        if (resid.le.tol) then
          if (mod(iter,n_iter_recompute_r).eq.0) exit
          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
 
          call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,u), ww(:,wk), tcomm)
          call hecmw_matvec(hecmesh, hecmat, ww(:,u), ww(:,v), tcomm)
          call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,u), cg(1))
          call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,v), ww(:,u), cg(2))
          call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,r), cg(3))
          start_time = hecmw_wtime()
          call hecmw_allreduce_r(hecmesh, cg, 3, hecmw_sum)
          end_time = hecmw_wtime()
          tcomm = tcomm + end_time-start_time
          gamma = cg(1)
          delta = cg(2)
          dnrm2 = cg(3)
          resid = dsqrt(dnrm2/bnrm2)
        endif
      endif
 
      if (gamma.eq.0.0d0) exit
      if (gamma.ne.gamma) then
        error = hecmw_solver_error_diverge_nan
        iter = i
        exit
      elseif (i.gt.1 .and. gamma*gamma1.le.0.0d0) then
        n_indef_precond = n_indef_precond + 1
        if (n_indef_precond.ge.3) then
          error = hecmw_solver_error_diverge_pc
          iter = i
          exit
        endif
      endif
 
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,v), ww(:,m), ww(:,wk), tcomm)
      call hecmw_matvec(hecmesh, hecmat, ww(:,m), ww(:,l), tcomm)
 
      if (i.gt.1) then
        beta = gamma/gamma1
        denom = delta-beta*gamma/alpha1
      else
        beta = 0.0d0
        denom = delta
      endif
 
      if (denom.le.0.0d0) then
        error = hecmw_solver_error_diverge_mat
        iter = i
        exit
      elseif (denom.ne.denom) then
        error = hecmw_solver_error_diverge_nan
        iter = i
        exit
      endif
      alpha = gamma/denom
 
      call hecmw_xpay_r(nndof, beta, ww(:,l), ww(:,z))
      call hecmw_xpay_r(nndof, beta, ww(:,m), ww(:,q))
      call hecmw_xpay_r(nndof, beta, ww(:,v), ww(:,s))
      call hecmw_xpay_r(nndof, beta, ww(:,u), ww(:,p))
 
      call hecmw_axpy_r(nndof, alpha, ww(:,p), x)
      if (mod(i,n_iter_recompute_r).eq.0) then
        call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
        call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,u), ww(:,wk), tcomm)
        call hecmw_matvec(hecmesh, hecmat, ww(:,u), ww(:,v), tcomm)
      else
        call hecmw_axpy_r(nndof, -alpha, ww(:,s), ww(:,r))
        call hecmw_axpy_r(nndof, -alpha, ww(:,q), ww(:,u))
        call hecmw_axpy_r(nndof, -alpha, ww(:,z), ww(:,v))
      endif
 
      iter = i
      if (estcond /= 0 .and. hecmesh%my_rank == 0) then
        if (i.eq.1) then
          d(1) = 1.0d0/alpha
        else
          d(i) = 1.0d0/alpha+beta/alpha1
          e(i-1) = dsqrt(beta)/alpha1
        endif
        if (mod(i,estcond).eq.0) call hecmw_estimate_condition_cg(i, d, e)
      endif
 
      gamma1 = gamma
      alpha1 = alpha
 
      if (i.eq.maxit) then
        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 (my_rank.eq.0 .and. iterlog.eq.1) write(*,'(i7, 1pe16.6)') iter, resid
        if (resid.gt.tol) error = hecmw_solver_error_noconv_maxit
      endif
    enddo
 
    call hecmw_solver_scaling_bk(hecmat)
 
    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)
 
    if (hecmw_mat_get_usejad(hecmat).ne.0) then
      call hecmw_jad_finalize(hecmat)
    endif
 
    if (estcond /= 0 .and. hecmesh%my_rank == 0) then
      if (error.eq.0 .and. iter.gt.0) 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_pipecg
 
end module hecmw_solver_PipeCG