gmres_mod.f90

Go to the documentation of this file.

module gmres_mod

  use real_type
  use basic_tools
  use abstract_interfaces, only: rntorn
  use algebra_lin, only: invtrisup
  use r1d_mod, only: norm_2, axpy, scalprod, xpay

  implicit none  
  private 
  
  public :: gmres   !! TESTED

contains


  subroutine gmres(x, iter, nbPrd, res, &
       & b, A, tol, itmax, rst, verb)

    real(RP), dimension(:), intent(inout) :: x
    integer               , intent(out)   :: iter, nbPrd
    real(RP)              , intent(out)   :: res
    procedure(rntorn)                     :: A
    real(RP), dimension(:), intent(in)    :: b
    real(RP)              , intent(in)    :: tol
    integer               , intent(in)    :: itmax, rst, verb

    real(RP), dimension(rst+1, size(x,1) ) :: V
    real(RP), dimension(rst+1, rst       ) :: H

    real(RP), dimension(size(x,1)) :: r, w
    real(RP), dimension(rst      ) :: sn, cs, y
    real(RP), dimension(rst +1   ) :: s

    real(RP) :: nb2, nr2, temp
    integer  :: nn, ii, kk

    if (verb>1) write(*,*) 'gmres_mod       : gmres'

    iter  = 0
    nbprd = 0


    nb2 = norm_2(b)
    nn = size(x,1)
    if (nb2/re(nn)<1e-8_rp) nb2=sqrt(re(nn))
    nb2 = 1._rp/nb2

    call a(r, x)
    nbprd = 1
    r   = b-r
    nr2 = norm_2(r)
    res = nr2 * nb2

    if (verb>2) write(*,*)'  iter', iter,' residual', res
    if (res<tol) return

    v  = 0._rp
    h  = 0._rp
    cs = 0._rp
    sn = 0._rp

    do iter=1, itmax

       if (verb>2) write(*,*)'  iter', iter,' residual', res

       s     = 0._rp
       s(1)  = nr2

       v(1,:)  = r / nr2
      
       do ii = 1, rst  ! orthonormal basis using Gram-Schmidt

          call a(w, v(ii,:))
          nbprd = nbprd + 1

          do kk = 1, ii
             h(kk,ii)= scalprod( w, v(kk,:))
             w = w - h(kk,ii)*v(kk,:)
          end do

          h(ii+1,ii) = norm_2(w)

          v(ii+1, : ) = w / h(ii+1,ii)

          ! apply Givens rotation
          !
          do kk = 1, ii-1 
             temp       =  cs(kk)*h(kk,ii) + sn(kk)*h(kk+1,ii)
             h(kk+1,ii) = -sn(kk)*h(kk,ii) + cs(kk)*h(kk+1,ii)
             h(kk,ii)   = temp
          end do

          call grotmat(cs(ii), sn(ii), h(ii,ii), h(ii+1,ii) )

          temp    =  cs(ii)*s(ii)             ! approximate residual norm
          s(ii+1) = -sn(ii)*s(ii)
          s(ii)   = temp

          h(ii,ii)   = cs(ii)*h(ii,ii) + sn(ii)*h(ii+1,ii)
          h(ii+1,ii) = 0._rp

          res  = abs(s(ii+1)) * nb2

          if (verb>2) write(*,*)'  iter', iter,' residual', res

          if ( res < tol ) then

             call invtrisup(y(1:ii), h(1:ii,1:ii), s(1:ii))
             do kk=1, ii
                x    = x + v(kk, : )*y(kk)
             end do

             return

          end if

       end do

       ! update approximation
       !
       call invtrisup(y(1:rst), h(1:rst,1:rst), s(1:rst))
       do kk=1, rst
          x    = x + v(kk, : )*y(kk)
       end do

       ! compute residual
       !
       call  a(r, x)                  
       nbprd = nbprd + 1
       r     = b-r
       nr2   = norm_2(r)
       res   = nr2 * nb2
       if (verb>2) write(*,*)'  iter', iter,' residual', res

       if ( res < tol ) return

    end do

    ! Convergence failure flag
    call warning("gmres_mod: gmres: not converged", verb-1)
    iter = -1

  end subroutine gmres

  subroutine grotmat(cs, sn, a, b)

    real(RP), intent(in)  :: a,b
    real(RP), intent(out) :: cs,sn

    real(RP)              :: tmp

    if (abs(b)<real_tol) then
       cs = 1._rp
       sn = 0._rp

    else if (abs(b)>abs(a)) then
       tmp = a/b
       sn  = 1._rp / sqrt( 1._rp + tmp**2)
       cs  = tmp*sn

    else
       tmp = b/a
       cs  = 1._rp / sqrt( 1._rp + tmp**2)
       sn  = tmp*cs

    end if

  end subroutine grotmat

end module gmres_mod