scale_gmres.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> Module common / GMRES
!!
!! @par Description
!!      A module to provide a iterative solver for system of linear equations using generalized minimal residual method (GMRES)
!!
!! @par Reference
!!  - Y. Saad and M.H. Schultz, 1986:
!!     GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems.
!!     SIAM J. Sci. Stat. Comput., 7:856–869
!!
!! @author Yuta Kawai, Team SCALE
!!
#include "scaleFElib.h"
module scale_gmres
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_prc
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type
  !  
 
  !> Derived type to provide a iterative solver for system of linear equations using GMRES
  type, public :: gmres
    real(rp) :: eps
    real(rp) :: eps0
    integer :: n
    integer :: m
    integer :: m_out
    real(rp), allocatable :: v(:,:)
    real(rp), allocatable :: hj(:)
    real(rp), allocatable :: g(:)
    real(rp), allocatable :: r(:,:)
    real(rp), allocatable :: co(:)
    real(rp), allocatable :: si(:)
    real(rp), allocatable :: y(:)
  contains 
    procedure, public :: init => gmres_init
    procedure, public :: iterate_pre => gmres_iterate_pre
    procedure, public :: iterate_step_j => gmres_iterate_step_j
    procedure, public :: iterate_post => gmres_iterate_post
    procedure, public :: final => gmres_final
  end type
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  
  !-----------------------------------------------------------------------------
 
  ! Private procedure
  !
 
  ! Private variable
  !
 
contains
 
subroutine gmres_init( this, N, m, eps, eps0 )
 
  implicit none
  class(gmres), intent(inout) :: this
  integer, intent(in) :: N
  integer, intent(in) :: m
  real(RP), intent(in) :: eps
  real(RP), intent(in) :: eps0
 
  !--------------------------------------
 
  this%N = n
  this%m = m
  this%EPS = eps
  this%EPS0 = eps0
 
  allocate( this%v(n,m+1) )
  allocate( this%hj(m+1) )
  allocate( this%g(m+1) )
  allocate( this%r(m+1,m) )
  allocate( this%co(m), this%si(m) )
  allocate( this%y(m) )
 
  return
end subroutine gmres_init
 
!OCL SERIAL
subroutine gmres_iterate_pre( this, b, w0, is_converged )
  implicit none
  class(gmres), intent(inout) :: this
  real(RP), intent(in) :: b(this%N)
  real(RP), intent(in) :: w0(this%N)
  logical, intent(out) :: is_converged
  
  integer :: i
  !--------------------------------------
 
  this%g(1) = 0.0_rp
  do i=1, this%N
    this%v(i,1) = b(i) - w0(i)
    this%g(1) = this%g(1) + this%v(i,1)**2
  end do
  this%g(1) = sqrt(this%g(1))
 
  if ( this%g(1) < this%EPS0 * this%N ) then
    is_converged = .true.
    return
  else
    is_converged = .false.
  end if
 
  do i=1, this%N
    this%v(i,1) = this%v(i,1) / this%g(1)
  end do
 
  this%m_out = min(this%m, this%N)
  return
end subroutine gmres_iterate_pre
 
!OCL SERIAL
subroutine gmres_iterate_step_j( this, j, wj, is_converged )
  implicit none
  class(gmres), intent(inout) :: this
  integer, intent(in) :: j
  real(RP), intent(inout) :: wj(this%N)
  logical, intent(out) :: is_converged
  
  integer :: i
  real(RP) :: tmp1, tmp2
  !--------------------------------------
 
  do i=1, j
    this%hj(i) = sum( wj(:) * this%v(:,i) )
    wj(:) = wj(:) - this%hj(i) * this%v(:,i)
  end do
  this%hj(j+1) = sqrt( sum(wj(:)**2) )
 
  ! if ( abs(this%hj(j+1)) < this%EPS0 ) then
  !   this%m_out = j
  !   is_converged = .true.
  !   return
  ! else
    this%v(:,j+1) = wj(:) / this%hj(j+1)
  ! end if
 
  this%r(1,j) = this%hj(1)
  do i=1, j-1
    tmp1 =   this%co(i) * this%r(i,j) + this%si(i) * this%hj(i+1)
    tmp2 = - this%si(i) * this%r(i,j) + this%co(i) * this%hj(i+1)
    this%r(i  ,j) = tmp1
    this%r(i+1,j) = tmp2
  end do
 
  tmp1 = 1.0_rp / sqrt(this%r(j,j)**2 + this%hj(j+1)**2)
  this%co(j) = tmp1 * this%r(j,j)
  this%si(j) = tmp1 * this%hj(j+1)
 
  this%g(j+1) = - this%si(j) * this%g(j)
  this%g(j)   =   this%co(j) * this%g(j)
 
  this%r(j,j) = this%co(j) * this%r(j,j) + this%si(j) * this%hj(j+1)
  this%r(j+1,j) = 0.0_rp
 
  if ( abs(this%g(j+1)) < this%EPS ) then
    this%m_out = j
    is_converged = .true.
  else
    is_converged = .false.
  end if
 
  return
end subroutine gmres_iterate_step_j
 
!OCL SERIAL
subroutine gmres_iterate_post( this, x )
  implicit none
  class(gmres), intent(inout) :: this
  real(RP), intent(inout) :: x(this%N)
  
  integer :: i, j
  !--------------------------------------
 
  do j=this%m_out, 1, -1
    this%y(j) = this%g(j)
    do i=j+1, this%m_out
      this%y(j) = this%y(j) - this%r(j,i) * this%y(i)
    end do
    this%y(j) = this%y(j) / this%r(j,j)
  end do
 
  x(:) = x(:) + matmul(this%v(:,1:this%m_out), this%y(1:this%m_out))
 
  return
end subroutine gmres_iterate_post
 
subroutine gmres_final( this )
  implicit none
  class(gmres), intent(inout) :: this
  !--------------------------------------
 
  if ( allocated(this%v) ) then
    deallocate( this%v )
    deallocate( this%hj )
    deallocate( this%g )
    deallocate( this%r )
    deallocate( this%co, this%si )
    deallocate( this%y )
  end if
 
  return
end subroutine gmres_final
 
end module scale_gmres