hecmw_solver_GMRES.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_GMRES
!C***
!
module hecmw_solver_gmres
 
  public :: hecmw_solve_gmres
 
contains
  !C
  !C*** hecmw_solve_GMRES
  !C
  subroutine hecmw_solve_gmres( 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 ) :: maxit, nrest
 
    real   (kind=kreal) :: tol
 
    real   (kind=kreal), dimension(:),   allocatable :: ss
    real   (kind=kreal), dimension(:,:), allocatable :: h
 
    integer(kind=kint ) :: cs, sn
 
    real   (kind=kreal)   zero, one
    parameter( zero = 0.0d+0, one = 1.0d+0 )
 
    integer(kind=kint ) :: nrk,i,k,kk,jj,info,ik
    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,val,vcs,vsn,dtemp,aa,bb,r0,scale,rr
    integer(kind=kint ) :: estcond
    real   (kind=kreal) :: t_max,t_min,t_avg,t_sd
 
    integer(kind=kint), parameter :: r  = 1
    integer(kind=kint), parameter :: zp = r + 1
    integer(kind=kint), parameter :: zq = r + 2
    integer(kind=kint), parameter :: s  = r + 3
    integer(kind=kint), parameter :: w  = s + 1
    integer(kind=kint), parameter :: y  = w
    integer(kind=kint), parameter :: av = y  + 1
    integer(kind=kint), parameter :: v  = av + 1
 
    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 )
    estcond = hecmw_mat_get_estcond( hecmat )
 
    if (nrest >= ndof*np-1) nrest = ndof*np-2
 
    error= 0
    nrk= nrest + 7
 
    allocate (h(nrk,nrk))
    allocate (ww(ndof*np,nrk))
    allocate (ss(nrk))
 
    commtime= 0.d0
    comptime= 0.d0
 
    ldh= nrest + 2
    ldw= n
 
    !C
    !C-- Store the Givens parameters in matrix H.
    cs= nrest + 1
    sn= cs    + 1
 
    !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)
 
    !C
    !C
    !C +--------------------+
    !C | {r}= {b} - [A]{x0} |
    !C +--------------------+
    !C===
    call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
    !C===
 
    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
 
    outer: do
 
      !C
      !C************************************************ GMRES Iteration
      !C
      i= 0
      !C
      !C +---------------+
      !C | {v1}= {r}/|r| |
      !C +---------------+
      !C===
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
      if (dnrm2 == 0.d0) exit ! converged
 
      rnorm= dsqrt(dnrm2)
      coef= one/rnorm
      call hecmw_axpby_r(nndof, coef, 0.d0, ww(:,r), ww(:,v))
      !C===
 
      !C
      !C +--------------+
      !C | {s}= |r|{e1} |
      !C +--------------+
      !C===
      call hecmw_scale_r(nndof, zero, ww(:,s))
      ww(1 ,s) = rnorm
      !C===
 
      !C************************************************ GMRES(m) restart
      do i = 1, nrest
        iter= iter + 1
 
        !C
        !C +-------------------+
        !C | {w}= [A][Minv]{v} |
        !C +-------------------+
        !C===
        call hecmw_precond_apply(hecmesh, hecmat, ww(:,v+i-1), ww(:,zq), ww(:,zp), tcomm)
 
        call hecmw_matvec(hecmesh, hecmat, ww(:,zq), ww(:,w), tcomm)
        !C===
 
        !C
        !C +------------------------------+
        !C | ORTH. BASIS by GRAMM-SCHMIDT |
        !C +------------------------------+
        !C   Construct the I-th column of the upper Hessenberg matrix H
        !C   using the Gram-Schmidt process on V and W.
        !C===
        do k= 1, i
          call hecmw_innerproduct_r(hecmesh, ndof, ww(:,w), ww(:,v+k-1), val, tcomm)
 
          call hecmw_axpy_r(nndof, -val, ww(:,v+k-1), ww(:,w))
          h(k,i)= val
        enddo
 
        call hecmw_innerproduct_r(hecmesh, ndof, ww(:,w), ww(:,w), val, tcomm)
        if (val == 0.d0) exit ! converged
 
        h(i+1,i)= dsqrt(val)
        coef= one / h(i+1,i)
        call hecmw_axpby_r(nndof, coef, 0.d0, ww(:,w), ww(:,v+i+1-1))
        !C===
 
        !C
        !C +-----------------+
        !C | GIVENS ROTARION |
        !C +-----------------+
        !C===
 
        !C
        !C-- Plane Rotation
        do k = 1, i-1
          vcs= h(k,cs)
          vsn= h(k,sn)
          dtemp   = vcs*h(k  ,i) + vsn*h(k+1,i)
          h(k+1,i)= vcs*h(k+1,i) - vsn*h(k  ,i)
          h(k  ,i)= dtemp
        enddo
 
        !C
        !C-- Construct Givens Plane Rotation
        aa = h(i  ,i)
        bb = h(i+1,i)
        r0= bb
        if (dabs(aa).gt.dabs(bb)) r0= aa
        scale= dabs(aa) + dabs(bb)
 
        if (scale.ne.0.d0) then
          rr= scale * dsqrt((aa/scale)**2+(bb/scale)**2)
          rr= dsign(1.d0,r0)*rr
          h(i,cs)= aa/rr
          h(i,sn)= bb/rr
        else
          h(i,cs)= 1.d0
          h(i,sn)= 0.d0
          rr     = 0.d0
        endif
 
        !C
        !C-- Plane Rotation
        vcs= h(i,cs)
        vsn= h(i,sn)
        dtemp    = vcs*h(i  ,i) + vsn*h(i+1,i)
        h(i+1,i)= vcs*h(i+1,i) - vsn*h(i  ,i)
        h(i  ,i)= dtemp
 
        dtemp    = vcs*ww(i  ,s) + vsn*ww(i+1,s)
        ww(i+1,s)= vcs*ww(i+1,s) - vsn*ww(i  ,s)
        ww(i  ,s)= dtemp
 
        resid = dabs( ww(i+1,s))/dsqrt(bnrm2)
 
        if (my_rank.eq.0 .and. iterlog.eq.1)                              &
          &    write (*, '(2i8, 1pe16.6)') iter,i+1, resid
 
        if (estcond /= 0 .and. hecmesh%my_rank == 0) then
          if (mod(iter,estcond) == 0) call hecmw_estimate_condition_gmres(i, h)
        endif
 
        if ( resid.le.tol ) then
          !C-- [H]{y}= {s_tld}
          do ik= 1, i
            ss(ik)= ww(ik,s)
          enddo
          irow= i
          ww(irow,y)= ss(irow) / h(irow,irow)
 
          do kk= irow-1, 1, -1
            do jj= irow, kk+1, -1
              ss(kk)= ss(kk) - h(kk,jj)*ww(jj,y)
            enddo
            ww(kk,y)= ss(kk) / h(kk,kk)
          enddo
 
          !C-- {x}= {x} + {y}{V}
          call hecmw_scale_r(nndof, 0.d0, ww(:,av))
 
          jj= irow
          do jj= 1, irow
            call hecmw_axpy_r(nndof, ww(jj,y), ww(:,v+jj-1), ww(:,av))
          enddo
 
          call hecmw_precond_apply(hecmesh, hecmat, ww(:,av), ww(:,zq), ww(:,zp), tcomm)
 
          call hecmw_axpy_r(nndof, 1.d0, ww(:,zq), x)
 
          exit outer
        endif
 
        if ( iter.gt.maxit ) then
          error = hecmw_solver_error_noconv_maxit
          exit outer
        end if
      end do
      !C===
 
      !C
      !C +------------------+
      !C | CURRENT SOLUTION |
      !C +------------------+
      !C===
 
      !C-- [H]{y}= {s_tld}
      do ik= 1, nrest
        ss(ik)= ww(ik,s)
      enddo
      irow= nrest
      ww(irow,y)= ss(irow) / h(irow,irow)
 
      do kk= irow-1, 1, -1
        do jj= irow, kk+1, -1
          ss(kk)= ss(kk) - h(kk,jj)*ww(jj,y)
        enddo
        ww(kk,y)= ss(kk) / h(kk,kk)
      enddo
 
      !C-- {x}= {x} + {y}{V}
      call hecmw_scale_r(nndof, 0.d0, ww(:,av))
 
      jj= irow
      do jj= 1, irow
        call hecmw_axpy_r(nndof, ww(jj,y), ww(:,v+jj-1), ww(:,av))
      enddo
 
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,av), ww(:,zq), ww(:,zp), tcomm)
 
      call hecmw_axpy_r(nndof, 1.d0, ww(:,zq), x)
 
      !C
      !C-- Compute residual vector R, find norm, then check for tolerance.
      call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
 
      call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
 
      ww(i+1,s)= dsqrt(dnrm2/bnrm2)
      resid    = ww( i+1,s )
 
      if ( resid.le.tol )   exit outer
      if ( iter .gt.maxit ) then
        error = hecmw_solver_error_noconv_maxit
        exit outer
      end if
      !C
      !C-- RESTART
    end do outer
 
    !C
    !C-- iteration FAILED
 
    if (error == hecmw_solver_error_noconv_maxit) then
      info = iter
 
      !C-- [H]{y}= {s_tld}
      do ik= 1, i
        ss(ik)= ww(ik,s)
      enddo
      irow= i
      ww(irow,y)= ss(irow) / h(irow,irow)
 
      do kk= irow-1, 1, -1
        do jj= irow, kk+1, -1
          ss(kk)= ss(kk) - h(kk,jj)*ww(jj,y)
        enddo
        ww(kk,y)= ss(kk) / h(kk,kk)
      enddo
 
      !C-- {x}= {x} + {y}{V}
      call hecmw_scale_r(nndof, 0.d0, ww(:,av))
 
      jj= irow
      do jj= 1, irow
        call hecmw_axpy_r(nndof, ww(jj,y), ww(:,v+jj-1), ww(:,av))
      enddo
 
      call hecmw_precond_apply(hecmesh, hecmat, ww(:,av), ww(:,zq), ww(:,zp), tcomm)
 
      call hecmw_axpy_r(nndof, 1.d0, ww(:,zq), x)
    end if
 
    call hecmw_solver_scaling_bk(hecmat)
 
    if (estcond /= 0 .and. hecmesh%my_rank == 0) then
      call hecmw_estimate_condition_gmres(i, h)
    endif
    !C
    !C-- INTERFACE data EXCHANGE
    s_time = hecmw_wtime()
    call hecmw_update_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
    e_time = hecmw_wtime()
    tcomm = tcomm + e_time - s_time
 
    deallocate (h, ww, ss)
    !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_gmres
 
end module     hecmw_solver_gmres