hecmw_solver_GroppCG.f90

Go to the documentation of this file.

!C
!C*** module hecmw_solver_GroppCG
!C
module hecmw_solver_groppcg
 
  public :: hecmw_solve_groppcg
 
contains
  !C
  !C*** GroppCG
  !C
  subroutine hecmw_solve_groppcg( 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 :: s  = 6
    integer(kind=kint), parameter :: wk = 7
 
    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
    real(kind=kreal) :: cg(2)
    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), parameter :: n_iter_recompute_r = 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)
 
    error = 0
    iter = 0
    resid = 0.0d0
    n_indef_precond = 0
    alpha1 = 0.0d0
    beta = 0.0d0
 
    allocate(ww(ndof*np, 7))
    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
      x = 0.0d0
    else
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,u), ww(:,wk), tcomm)
      call hecmw_copy_r(nndof, ww(:,u), ww(:,p))
      call hecmw_matvec(hecmesh, hecmat, ww(:,p), ww(:,s), tcomm)
 
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,u), cg(1))
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,r), cg(2))
      start_time = hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, cg, 2, hecmw_sum)
      end_time = hecmw_wtime()
      tcomm = tcomm + end_time-start_time
      gamma = cg(1)
      dnrm2 = cg(2)
      resid = dsqrt(dnrm2/bnrm2)
      if (resid.le.tol) maxit = 0
    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
      if (gamma.eq.0.0d0) then
        iter = i
        error = hecmw_solver_error_diverge_pc
        exit
      elseif (gamma.ne.gamma) then
        iter = i
        error = hecmw_solver_error_diverge_nan
        exit
      endif
 
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,p), ww(:,s), delta, tcomm)
      if (delta.le.0.0d0) then
        iter = i
        error = hecmw_solver_error_diverge_mat
        exit
      elseif (delta.ne.delta) then
        iter = i
        error = hecmw_solver_error_diverge_nan
        exit
      endif
 
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,s), ww(:,q), ww(:,wk), tcomm)
 
      alpha = gamma/delta
 
      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
 
      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)
      else
        call hecmw_axpy_r(nndof, -alpha, ww(:,s), ww(:,r))
        call hecmw_axpy_r(nndof, -alpha, ww(:,q), ww(:,u))
      endif
 
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,u), cg(1))
      call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,r), ww(:,r), cg(2))
      start_time = hecmw_wtime()
      call hecmw_allreduce_r(hecmesh, cg, 2, hecmw_sum)
      end_time = hecmw_wtime()
      tcomm = tcomm + end_time-start_time
      gamma1 = cg(1)
      dnrm2 = cg(2)
      resid = dsqrt(dnrm2/bnrm2)
      iter = i
 
      if (resid.le.tol) then
        if (mod(i,n_iter_recompute_r).ne.0) 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)
        endif
        if (my_rank.eq.0 .and. iterlog.eq.1) write(*,'(i7, 1pe16.6)') iter, resid
        if (resid.le.tol) exit
 
        call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,u), ww(:,wk), tcomm)
        call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,u), gamma1, tcomm)
      else
        if (my_rank.eq.0 .and. iterlog.eq.1) write(*,'(i7, 1pe16.6)') iter, resid
      endif
 
      if (i.eq.maxit) then
        if (mod(i,n_iter_recompute_r).ne.0) 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)
        endif
        if (resid.gt.tol) error = hecmw_solver_error_noconv_maxit
        exit
      endif
 
      if (gamma1.eq.0.0d0) then
        error = hecmw_solver_error_diverge_pc
        exit
      elseif (gamma1.ne.gamma1) then
        error = hecmw_solver_error_diverge_nan
        exit
      elseif (gamma1*gamma.le.0.0d0) then
        n_indef_precond = n_indef_precond + 1
        if (n_indef_precond.ge.3) then
          error = hecmw_solver_error_diverge_pc
          exit
        endif
      endif
 
      call hecmw_matvec(hecmesh, hecmat, ww(:,u), ww(:,v), tcomm)
 
      beta = gamma1/gamma
      call hecmw_xpay_r(nndof, beta, ww(:,u), ww(:,p))
      call hecmw_xpay_r(nndof, beta, ww(:,v), ww(:,s))
 
      gamma = gamma1
      alpha1 = alpha
    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_groppcg
 
end module hecmw_solver_GroppCG