hecmw_solver_GMRESR.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_GMRESR
!C***
!
module hecmw_solver_gmresr
 
  public :: hecmw_solve_gmresr
 
contains
  !C
  !C*** hecmw_solve_GMRESR
  !C
  subroutine hecmw_solve_gmresr( 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(:)
 
    integer(kind=kint),dimension(:), allocatable :: idxbfgs
    real(kind=kreal), dimension(:)  ,  allocatable :: vecr,workpc,tmpvecbfgs,rho,alpha
    real(kind=kreal), dimension(:,:),  allocatable :: u,c,uin,cin,sbfgs,ybfgs
 
    integer(kind=kint ) :: maxit, nrest,nbfgs
 
    real   (kind=kreal) :: tol
 
    real   (kind=kreal)   zero, one
    parameter( zero = 0.0d+0, one = 1.0d+0 )
 
    integer(kind=kint ) :: nrk,i,k,kk,jj,info,ik,iorth,idx,tmpidx,ibfgs
    integer(kind=kint ) :: irow
    real   (kind=kreal) :: s_time,e_time,s1_time,e1_time
    real   (kind=kreal) :: ldh,ldw,bnrm2,dnrm2,rnorm
    real   (kind=kreal) :: commtime,comptime, coef,coef2,val,vcs,vsn,dtemp,aa,bb,r0,scale,rr
    integer(kind=kint ) :: estcond
    real   (kind=kreal) :: t_max,t_min,t_avg,t_sd
 
 
    call hecmw_barrier(hecmesh)
    s_time= hecmw_wtime()
    !C
    !C-- INIT.
    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 )
    nrest  = hecmw_mat_get_nrest( hecmat )
    nbfgs  = hecmw_mat_get_nbfgs( hecmat )
    estcond = hecmw_mat_get_estcond( hecmat )
 
    error= 0
 
    allocate (vecr(ndof*np))
    allocate (workpc(ndof*np))
    allocate (u(ndof*np,nrest))
    allocate (c(ndof*np,nrest))
    allocate (uin(ndof*np,nrest))
    allocate (cin(ndof*np,nrest))
 
    if(nbfgs>0)then
       allocate (tmpvecbfgs(ndof*np))
       allocate (sbfgs(ndof*np,nbfgs))
       allocate (ybfgs(ndof*np,nbfgs))
       allocate (idxbfgs(nbfgs))
       allocate (rho(nbfgs))
       allocate (alpha(nbfgs))
    endif
 
    do idx = 1, nbfgs
       idxbfgs(idx) = idx
    enddo
 
    commtime= 0.d0
    comptime= 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)
 
    !C===
    !C +----------------------+
    !C | SETUP PRECONDITIONER |
    !C +----------------------+
    !C===
    call hecmw_precond_setup(hecmat, hecmesh, 0)
 
 
    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
    !C===
 
 
    call hecmw_barrier(hecmesh)
    s1_time= hecmw_wtime()
    iter= 0
    ibfgs=0
    outer: do
 
      call hecmw_matresid(hecmesh, hecmat, x, b, vecr, tcomm)
      do i = 1, nrest
        iter= iter + 1
 
        !C Solve M*r = uin(:,1)
        if (ibfgs == 0)then
           call hecmw_precond_apply(hecmesh, hecmat, vecr, uin(:,1), workpc, tcomm)
        else
           call hecmw_copy_r(nndof, vecr, tmpvecbfgs)
           do k = 1,ibfgs
              idx = idxbfgs(k)
              call hecmw_innerproduct_r(hecmesh, ndof, sbfgs(:,idx), ybfgs(:,idx), coef, tcomm)
              rho(k) = 1.0d0 / coef
              call hecmw_innerproduct_r(hecmesh, ndof, sbfgs(:,idx), tmpvecbfgs, coef2, tcomm)
              alpha(k) = rho(k)*coef2
              call hecmw_axpy_r(nndof, -alpha(k), ybfgs(:,idx), tmpvecbfgs)
           enddo
           call hecmw_precond_apply(hecmesh, hecmat, tmpvecbfgs, uin(:,1), workpc, tcomm)
           do k = ibfgs,1,-1
              idx = idxbfgs(k)
              call hecmw_innerproduct_r(hecmesh, ndof, ybfgs(:,idx), uin(:,1), coef, tcomm)
              coef2 = rho(k) * coef
              call hecmw_axpy_r(nndof, alpha(k)-coef2, sbfgs(:,idx), uin(:,1))
           enddo
        endif
        !C cin(:,1) = A*uin(:,1)
        call hecmw_matvec(hecmesh, hecmat, uin(:,1), cin(:,1), tcomm)
 
        do iorth = 1, i-1
           !C c_{i}^T cin_{i}
           call hecmw_innerproduct_r(hecmesh, ndof, c(:,iorth), cin(:,iorth), coef, tcomm)
 
           call hecmw_axpyz_r(nndof, -coef, c(:,iorth), cin(:,iorth), cin(:,iorth+1))
           call hecmw_axpyz_r(nndof, -coef, u(:,iorth), uin(:,iorth), uin(:,iorth+1))
        enddo
        call hecmw_innerproduct_r(hecmesh, ndof, cin(:,i), cin(:,i), coef, tcomm)
        coef = 1.0d0 / dsqrt(coef)
        call hecmw_axpby_r(nndof, coef, 0.0d0, cin(:,i), c(:,i))
        call hecmw_axpby_r(nndof, coef, 0.0d0, uin(:,i), u(:,i))
 
        call hecmw_innerproduct_r(hecmesh, ndof, c(:,i), vecr, coef, tcomm)
        call hecmw_axpy_r(nndof,  coef, u(:,i),    x)
        call hecmw_axpy_r(nndof, -coef, c(:,i), vecr)
 
        if (nbfgs > 0)then
           ibfgs = ibfgs + 1
           if (ibfgs == nbfgs+1)then
             tmpidx = idxbfgs(1)
             do kk = 1, nbfgs-1
                idxbfgs(kk) = idxbfgs(kk+1)
             enddo
             idxbfgs(nbfgs) = tmpidx
             ibfgs = ibfgs - 1
           endif
           call hecmw_axpby_r(nndof, coef, 0.0d0, c(:,i), ybfgs(:,idxbfgs(ibfgs)))
           call hecmw_axpby_r(nndof, coef, 0.0d0, u(:,i), sbfgs(:,idxbfgs(ibfgs)))
        endif
 
 
        call hecmw_innerproduct_r(hecmesh, ndof, vecr, vecr, 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 ( resid.le.tol )   exit outer
        if ( iter.gt.maxit ) then
          error = hecmw_solver_error_noconv_maxit
          exit outer
        end if
      end do
 
    end do outer
 
    call hecmw_solver_scaling_bk(hecmat)
 
    !C
    !C-- INTERFACE data EXCHANGE
    s_time = hecmw_wtime()
    !call hecmw_update_m_R (hecMESH, X, hecMAT%NP, hecMAT%NDOF)
    call hecmw_update_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
 
    e_time = hecmw_wtime()
    tcomm = tcomm + e_time - s_time
 
    !deallocate (H, WW, SS)
    deallocate (vecr)
    deallocate (workpc)
    deallocate (u  )
    deallocate (c  )
    deallocate (uin)
    deallocate (cin)
    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
 
    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_gmresr
 
end module     hecmw_solver_gmresr