hecmw_solver_Iterative.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
 
module hecmw_solver_iterative
contains
  !
  !C***
  !C*** hecmw_solve_nn
  !C***
  !
  subroutine hecmw_solve_iterative (hecMESH, hecMAT)
 
    use hecmw_util
    use hecmw_solver_cg
    use hecmw_solver_pipecg
    use hecmw_solver_groppcg
    use hecmw_solver_bicgstab
    use hecmw_solver_gmres
    use hecmw_solver_gmresr
    use hecmw_solver_gmresren
    use hecmw_solver_gpbicg
    use hecmw_solver_cr
    use m_hecmw_solve_error
    use hecmw_solver_las
    use hecmw_precond
    use hecmw_matrix_misc
    use hecmw_matrix_dump
    use hecmw_solver_misc
 
    implicit none
 
    type (hecmwST_matrix), target :: hecMAT
    type (hecmwST_local_mesh) :: hecMESH
 
    integer(kind=kint) :: error
    integer(kind=kint) :: ITER, METHOD, PRECOND, NSET, METHOD2
    integer(kind=kint) :: iterPREmax
    integer(kind=kint) :: ITERlog, TIMElog
    real(kind=kreal)   :: resid, sigma_diag, thresh, filter, resid2
    real(kind=kreal)   :: time_setup, time_comm, time_sol, tr
    real(kind=kreal)   :: time_ax, time_precond
 
    integer(kind=kint) :: NREST
    real(kind=kreal)   :: omega
 
    integer(kind=kint) :: auto_sigma_diag
 
    !C PARAMETERs
    iter      = hecmw_mat_get_iter(hecmat)
    method    = hecmw_mat_get_method(hecmat)
    method2   = hecmw_mat_get_method2(hecmat)
    precond   = hecmw_mat_get_precond(hecmat)
    nset      = hecmw_mat_get_nset(hecmat)
    iterpremax= hecmw_mat_get_iterpremax(hecmat)
    nrest     = hecmw_mat_get_nrest(hecmat)
    iterlog   = hecmw_mat_get_iterlog(hecmat)
    timelog   = hecmw_mat_get_timelog(hecmat)
    time_setup= 0.d0
    time_comm = 0.d0
    time_sol  = 0.d0
    resid     = hecmw_mat_get_resid(hecmat)
    sigma_diag= hecmw_mat_get_sigma_diag(hecmat)
    thresh    = hecmw_mat_get_thresh(hecmat)
    filter    = hecmw_mat_get_filter(hecmat)
    if (sigma_diag.lt.0.d0) then
      auto_sigma_diag= 1
      sigma_diag= 1.d0
    else
      auto_sigma_diag= 0
    endif
 
    ! the control file reaches Rarray without passing through the setter, so the
    ! range that hecmw_mat_set_omega enforces has to be applied here as well
    omega     = hecmw_mat_get_omega(hecmat)
    call hecmw_mat_set_omega(hecmat, omega)
 
    !C ERROR CHECK
    call hecmw_solve_check_zerodiag(hecmesh, hecmat) !C-- ZERO DIAGONAL component
 
    !C-- RECYCLE SETTING OF PRECONDITIONER
    call hecmw_mat_recycle_precond_setting(hecmat)
 
    ! exchange diagonal elements of overlap region
    call hecmw_mat_dump(hecmat, hecmesh)
 
    !C ITERATIVE solver
    error=0
    !! Auto Sigma_diag loop
    do
      call hecmw_mat_set_flag_converged(hecmat, 0)
      call hecmw_mat_set_flag_diverged(hecmat, 0)
      if (auto_sigma_diag.eq.1) call hecmw_mat_set_sigma_diag(hecmat, sigma_diag)
 
      call hecmw_matvec_clear_timer()
      call hecmw_precond_clear_timer()
      call hecmw_solve_iterative_printmsg(hecmesh,hecmat,method)
 
      select case(method)
        case (1)  !--CG
          hecmat%symmetric = .true.
          call hecmw_solve_cg( hecmesh, hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (2)  !--BiCGSTAB
          hecmat%symmetric = .false.
          call hecmw_solve_bicgstab( hecmesh,hecmat, iter, resid, error,time_setup, time_sol, time_comm )
        case (3)  !--GMRES
          hecmat%symmetric = .false.
          call hecmw_solve_gmres( hecmesh,hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (4)  !--GPBiCG
          hecmat%symmetric = .false.
          call hecmw_solve_gpbicg( hecmesh,hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (5)  !--GMRESR
          hecmat%symmetric = .false.
          call hecmw_solve_gmresr( hecmesh,hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (6)  !--GMRESREN
          hecmat%symmetric = .false.
          call hecmw_solve_gmresren( hecmesh,hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (7)  !--CR
          hecmat%symmetric = .true.
          call hecmw_solve_cr( hecmesh, hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (8)  !--PipeCG
          hecmat%symmetric = .true.
          call hecmw_solve_pipecg( hecmesh, hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case (9)  !--GroppCG
          hecmat%symmetric = .true.
          call hecmw_solve_groppcg( hecmesh, hecmat, iter, resid, error, time_setup, time_sol, time_comm )
        case default
          error = hecmw_solver_error_incons_pc  !!未定義なMETHOD!!
          call hecmw_solve_error (hecmesh, error)
      end select
 
      if (error==hecmw_solver_error_diverge_pc .or. error==hecmw_solver_error_diverge_mat &
           .or. error==hecmw_solver_error_diverge_nan) then
        call hecmw_mat_set_flag_diverged(hecmat, 1)
        if ((precond>=10 .and. precond<20) .and. auto_sigma_diag==1 .and. sigma_diag<2.d0) then
          sigma_diag = sigma_diag + 0.1
          if (hecmesh%my_rank.eq.0) write(*,*) 'Increasing SIGMA_DIAG to', sigma_diag
          cycle
        elseif ((method==1 .or. method==8 .or. method==9) .and. method2>1) then
          if (auto_sigma_diag.eq.1) sigma_diag = 1.0
          method = method2
          cycle
        endif
      endif
 
      if (auto_sigma_diag.eq.1) call hecmw_mat_set_sigma_diag(hecmat, -1.d0)
      exit
    enddo
 
    if (error.ne.0) then
      call hecmw_solve_error (hecmesh, error)
    endif
 
    resid2=hecmw_rel_resid_l2(hecmesh,hecmat)
    if (hecmesh%my_rank.eq.0 .and. (iterlog.eq.1 .or. timelog.ge.1)) then
      write(*,"(a,1pe12.5)")'### Relative residual =', resid2
    endif
    if (resid2 < hecmw_mat_get_resid(hecmat)) call hecmw_mat_set_flag_converged(hecmat, 1)
 
    call hecmw_mat_dump_solution(hecmat)
 
    time_ax = hecmw_matvec_get_timer()
    time_precond = hecmw_precond_get_timer()
 
    if (hecmesh%my_rank.eq.0 .and. timelog.ge.1) then
      tr= (time_sol-time_comm)/(time_sol+1.d-24)*100.d0
      write (*,'(/a)')          '### summary of linear solver'
      write (*,'(i10,a, 1pe16.6)')      iter, ' iterations  ', resid
      write (*,'(a, 1pe16.6 )') '    set-up time      : ', time_setup
      write (*,'(a, 1pe16.6 )') '    solver time      : ', time_sol
      write (*,'(a, 1pe16.6 )') '    solver/comm time : ', time_comm
      write (*,'(a, 1pe16.6 )') '    solver/matvec    : ', time_ax
      write (*,'(a, 1pe16.6 )') '    solver/precond   : ', time_precond
      if (iter > 0) &
        write (*,'(a, 1pe16.6 )') '    solver/1 iter    : ', time_sol / iter
      write (*,'(a, 1pe16.6/)') '    work ratio (%)   : ', tr
    endif
 
    call hecmw_output_flops(hecmesh, hecmat, iter, time_ax)
 
  end subroutine hecmw_solve_iterative
 
  subroutine hecmw_output_flops(hecMESH, hecMAT, count_Ax, time_Ax)
    use hecmw_util
    use m_hecmw_comm_f
    use hecmw_matrix_misc
    use hecmw_solver_misc
    implicit none
    type (hecmwST_local_mesh) :: hecMESH
    type(hecmwst_matrix) :: hecMAT
    integer(kint) :: N, NP, NDOF, NPU, NPL, NZ
    integer(kint) :: base, i, count_Ax
    real(kreal) :: time_Ax, size_matrix, flop_matrix, size_vector, memory_size, tmp, num
    real(kreal) :: t_max, t_min, t_avg, t_sd
    character(2) :: SI(0:6) = ['  ',' K',' M',' G',' T',' P',' E']
 
    if(hecmw_mat_get_timelog(hecmat) /= 2) return
 
    n = hecmat%N
    np = hecmat%NP
    ndof = hecmat%NDOF
    npu = hecmat%indexU(n)
    npl = hecmat%indexL(n)
    nz = n + npu + npl
 
    size_matrix = kreal*nz*ndof**2 &
                + kint*(npu+npl) &
                + kint*2*(n+1)
    size_vector = kreal*n*ndof &
                + kreal*np*ndof
    memory_size = size_matrix + size_vector
    flop_matrix = 2.0d0*nz*ndof**2
 
    call hecmw_allreduce_r1(hecmesh, memory_size, hecmw_sum)
    call hecmw_allreduce_r1(hecmesh, flop_matrix, hecmw_sum)
 
    base = 1000 ! or 1024
    num = memory_size
    i = int(log(num) / log(dble(base)))
    tmp = 1.0d0/base**i * num
    if(hecmesh%my_rank == 0)then
      write (*,"(a,f11.3,a,a)") "memory amount of coef. matrix: ", tmp, si(i),"B"
    endif
 
    num = count_ax*memory_size/time_ax
    i = int(log(num) / log(dble(base)))
    tmp = 1.0d0/base**i * num
    call hecmw_time_statistics(hecmesh, tmp, t_max, t_min, t_avg, t_sd)
    if(hecmesh%my_rank == 0)then
      write(*,"(a,f11.3,a,a)") "matvec memory band width     : ", tmp, si(i),"B/s"
      write(*,"(a,f11.3)") '  Max     :',t_max
      write(*,"(a,f11.3)") '  Min     :',t_min
      write(*,"(a,f11.3)") '  Avg     :',t_avg
      write(*,"(a,f11.3)") '  Std Dev :',t_sd
    endif
 
    num = count_ax*flop_matrix/time_ax
    i = int(log(num) / log(dble(base)))
    tmp = 1.0d0/base**i * num
    call hecmw_time_statistics(hecmesh, tmp, t_max, t_min, t_avg, t_sd)
    if(hecmesh%my_rank == 0)then
      write(*,"(a,f11.3,a,a)") "matvec FLOPs                 : ", tmp, si(i),"FLOPs"
      write(*,"(a,f11.3)") '  Max     :',t_max
      write(*,"(a,f11.3)") '  Min     :',t_min
      write(*,"(a,f11.3)") '  Avg     :',t_avg
      write(*,"(a,f11.3)") '  Std Dev :',t_sd
    endif
  end subroutine hecmw_output_flops
 
  subroutine hecmw_solve_check_zerodiag (hecMESH, hecMAT)
    use hecmw_util
    use hecmw_matrix_misc
    use m_hecmw_solve_error
    use m_hecmw_comm_f
    implicit none
    type (hecmwST_local_mesh) :: hecMESH
    type (hecmwST_matrix), target :: hecMAT
    integer (kind=kint)::PRECOND,iterPREmax,i,j,error
    precond   = hecmw_mat_get_precond(hecmat)
    iterpremax= hecmw_mat_get_iterpremax(hecmat)
    !C
    !C-- ZERO DIAGONAL component
    error= 0
    do i= 1, hecmat%N
      do j = 1, hecmat%NDOF
        if (dabs(hecmat%D(hecmat%NDOF*hecmat%NDOF*(i-1)+(j-1)*(hecmat%NDOF+1)+1)).eq.0.d0) then
          error=hecmw_solver_error_zero_diag
        end if
      end do
    enddo
 
    call hecmw_allreduce_i1 (hecmesh, error, hecmw_max)
    if (error.ne.0 .and. (precond.lt.10 .and. iterpremax.gt.0)) then
      call hecmw_solve_error (hecmesh, error)
    endif
 
  end subroutine hecmw_solve_check_zerodiag
 
  function hecmw_solve_check_zerorhs (hecMESH, hecMAT)
    use hecmw_util
    use hecmw_matrix_misc
    use m_hecmw_solve_error
    use m_hecmw_comm_f
    implicit none
    type (hecmwst_local_mesh) :: hecmesh
    type (hecmwst_matrix), target :: hecmat
    real(kind=kreal),    dimension(1) :: rhs
    integer (kind=kint)::precond,iterpremax,i,j,error
    logical :: hecmw_solve_check_zerorhs
 
    precond   = hecmw_mat_get_precond(hecmat)
    iterpremax= hecmw_mat_get_iterpremax(hecmat)
    !C
    !C-- ZERO RHS norm
    error= 0
    hecmw_solve_check_zerorhs = .false.
 
    rhs(1)= 0.d0
    do i= 1, hecmat%N
      do j = 1, hecmat%NDOF
        rhs(1)=rhs(1) + hecmat%B(hecmat%NDOF*(i-1)+j)**2
      end do
    enddo
    if (hecmesh%mpc%n_mpc > 0) then
      do i= 1, hecmesh%mpc%n_mpc
        rhs(1)= rhs(1) + hecmesh%mpc%mpc_const(i)**2
      enddo
    endif
    call hecmw_allreduce_r (hecmesh, rhs, 1, hecmw_sum)
 
    if (rhs(1).eq.0.d0) then
      error= hecmw_solver_error_zero_rhs
      call hecmw_solve_error (hecmesh, error)
      hecmat%X(:)=0.d0
      hecmw_solve_check_zerorhs = .true.
    endif
 
  end function hecmw_solve_check_zerorhs
 
  subroutine hecmw_solve_iterative_printmsg (hecMESH, hecMAT, METHOD)
    use hecmw_util
    use hecmw_solver_misc
    use hecmw_matrix_misc
 
    implicit none
    type (hecmwST_local_mesh) :: hecMESH
    type (hecmwST_matrix), target :: hecMAT
    integer(kind=kint) :: METHOD
    integer(kind=kint) :: ITER, PRECOND, NSET, iterPREmax, NREST,NBFGS
    integer(kind=kint) :: ITERlog, TIMElog
 
    character(len=30) :: msg_precond
    character(len=30) :: msg_method
 
    iter      = hecmw_mat_get_iter(hecmat)
    ! METHOD    = hecmw_mat_get_method(hecMAT)
    precond   = hecmw_mat_get_precond(hecmat)
    nset      = hecmw_mat_get_nset(hecmat)
    iterpremax= hecmw_mat_get_iterpremax(hecmat)
    nrest     = hecmw_mat_get_nrest(hecmat)
    nbfgs     = hecmw_mat_get_nbfgs(hecmat)
    iterlog= hecmw_mat_get_iterlog(hecmat)
    timelog= hecmw_mat_get_timelog(hecmat)
 
    select case(method)
      case (1)  !--CG
        msg_method="CG"
      case (2)  !--BiCGSTAB
        msg_method="BiCGSTAB"
      case (3)  !--GMRES
        msg_method="GMRES"
      case (4)  !--GPBiCG
        msg_method="GPBiCG"
      case (5)  
        if (nbfgs==0) then
           msg_method="GMRESR"
        else
           msg_method="SUP-GMRESR"
        endif
      case (6)  
        msg_method="GMRESR-EN"
      case (7)  !--CR
        msg_method="CR"
      case (8)  !--PipeCG
        msg_method="PipeCG"
      case (9)  !--GroppCG
        msg_method="GroppCG"
      case default
        msg_method="Unlabeled"
    end select
    select case(precond)
      case (1,2)
        msg_precond="SSOR"
      case (3)
        msg_precond="DIAG"
      case (5)
        msg_precond="ML"
      case (7)
        msg_precond="DirectMUMPS"
      case (10, 11, 12)
        write(msg_precond,"(a,i0,a)") "ILU(",precond-10,")"
      case (20)
        msg_precond="SAINV"
      case (21)
        msg_precond="RIF"
      case (22)
        msg_precond="SAAMG"
      case default
        msg_precond="Unlabeled"
    end select
    if (hecmesh%my_rank.eq.0 .and. (iterlog.eq.1 .or. timelog.ge.1)) then
      write (*,'(a,i0,a,i0,a,a,a,a,a,i0)') '### ',hecmat%NDOF,'x',hecmat%NDOF,' BLOCK ', &
        &   trim(msg_method),", ",trim(msg_precond),", ", iterpremax
    end if
  end subroutine hecmw_solve_iterative_printmsg
 
end module hecmw_solver_iterative