cg_mod.f90

Go to the documentation of this file.

module cg_mod

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

  implicit none  
  private 


!       %----------------------------------------%
!       |                                        |
!       |          PUBLIC DATA                   |
!       |                                        |
!       %----------------------------------------%

  !! GENERIC
  public :: cg         !! TESTED
  public :: pcg        !! TESTED

contains


  subroutine cg(x, iter, res, b, A, tol, itmax, verb)

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


    real(RP), dimension(size(x,1)) :: r, p, Ap
    real(RP)   :: alpha, beta, nb2
    real(RP)   :: ps_rr, ps_App, ps_rjrj
    integer    :: nn

    if (verb>1) write(*,*)'cg_mod          : cg'

    iter = 0

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

    call a(r, x)
    ! r   = b-r
    call axpy(-1._rp, r, b)

    res = norm_2(r) * nb2

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

    ! p     = r
    p = r
    ps_rr = scalprod(r, r)
    do iter=1, itmax
       ps_rjrj = ps_rr

       call a(ap, p)

       ps_app = scalprod(ap, p)

       alpha = ps_rr / ps_app
       ! x   = x + alpha * p
       call xpay(x, alpha, p)

       ! r   = r - alpha * Ap
       call xpay(r, -alpha, ap)

       res   = norm_2(r) * nb2

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

       ps_rr = scalprod(r, r)
       beta  = ps_rr/ps_rjrj
       ! p  = r + beta * p
       call axpy(beta, p, r)

    end do

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

  end subroutine cg


  subroutine pcg(x, iter, res, b, A, pc, tol, itmax, verb)

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

    real(RP), dimension(size(x,1)) :: r, p, Ap, z

    real(RP) :: alpha, beta, nb2
    real(RP) :: ps_rz, ps_App, ps_rjzj
    integer  :: nn

    if (verb>1) write(*,*)'cg_mod          : pcg'

    iter = 0

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

    call a(r, x)
    ! r   = b-r
    call axpy(-1._rp, r, b)    
    res = norm_2(r) * nb2

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

    if (res<tol) return

    call pc(z, r)
    ! p = z
    call copy(p, z)
    ps_rz = scalprod(r,z)
    do iter=1, itmax
       ps_rjzj = ps_rz

       call a(ap, p)
       
       ps_app = scalprod(ap,p)

       alpha = ps_rz / ps_app
       ! x   = x + alpha * p
       call xpay(x, alpha, p)
       ! r   = r - alpha * Ap
       call xpay(r, -alpha, ap)

       res   = norm_2(r) * nb2

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

       if (res<tol) return

       call pc(z, r)
       ps_rz = scalprod(r,z)
       beta  = ps_rz/ps_rjzj
       !p     = z + beta * p
       call axpy(beta, p, z)

    end do

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

  end subroutine pcg

end module cg_mod