newton_mod.f90

Go to the documentation of this file.

module newton_mod

  use choral_constants
  use basic_tools
  use real_type
  use abstract_interfaces, only: rntorn
  use r1d_mod, only: norm_2
  use krylov_mod

  implicit none  
  private 

  
  public :: solve, newton

  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type newton

     real(RP)    :: tol = 1e-5_rp, eps = 1e-3_rp
     integer     :: verb=0, itmax=100
     type(krylov):: kry
     
     !! Data returned when solving
     !!
     real(RP) :: res
     integer  :: iter, feval=0
     logical  :: ierr = .false.

  end type newton

  interface newton
     module procedure newton_create
  end interface newton

  interface solve
     module procedure newton_solve
  end interface solve

contains


  function newton_create(tol, itmax, eps, verb) result(n)
    type(newton)                   :: n
    real(RP), intent(in), optional :: tol, eps
    integer , intent(in), optional :: itmax, verb

    if (present(verb )) n%verb = verb
    if (present(tol  )) n%tol  = tol
    if (present(itmax)) n%itMax= itmax
    if (present(eps  )) n%eps  = eps

    n%kry = krylov(type=kry_gmres, tol=1e-4_rp, &
         &   itmax=100, restart=3)

  end function newton_create


  subroutine newton_solve(x, n, F)
    real(RP), dimension(:), intent(inout) :: x
    type(newton)            , intent(inout) :: n
    procedure(rntorn)                       :: F

    real(RP), dimension(size(x, 1)) :: Fx, z, x0
    real(DP) :: cpu
    integer :: ii

    n%ierr = .false.

    if (n%verb>0) then
       write(*,*)'newton_mod: newton_solve: '
       cpu = clock() 
    end if

    ! Initiation
    !
    n%iter = 0
    x0   = x
    call f(fx, x)
    n%Feval   = 1
    n%res  = residual(fx, x)
    if (n%verb>1) write(*,*)'  Iter', n%iter,', residual', n%res
    if (n%res<n%tol) goto 10

    ! Newton loop
    !
    do  ii=1, n%itmax
            
       z = x

       call solve(z, n%kry, fx, df)

       ! Updates
       !
       x  = x - z

       call f(fx, x)
       n%Feval = n%Feval + 1 + n%kry%AEval

       n%res = residual(fx, x)    ! exit loop test
       if (n%verb>1) &
            &    write(*,*)'  Iter', ii,', residual', n%res

       ! Test loop end
       !
       if (n%res<n%tol) goto 10

       x0 = x

    end do

    ! Convergence failure flag
    call warning("newton_mod: newton_solve: not converged", n%verb-1)
    n%ierr = .true.

10  continue
    n%iter = ii
    if (n%verb>0) then
       write(*,*)"  Iter               ",n%iter
       write(*,*)"  Performed x-->F(x) ",n%Feval
       write(*,*)"  Residual           ",n%res

       cpu = clock() - cpu
       write(*,*)'  CPU          ', real(cpu, sp)
    end if

  contains


    ! z = Df(x).y \simeq (F(x + eps*y) - F(x))/eps
    !
    subroutine df(zz, yy)
      
      real(RP), dimension(:), intent(out) :: zz
      real(RP), dimension(:), intent(in)  :: yy

      call f(zz,x + n%eps*yy)
      zz = (zz - fx)/n%eps
      
    end subroutine df
    
  end subroutine newton_solve
  

  
  function residual(Fx, x) result(res)

    real(RP) :: res
    real(RP), dimension(:), intent(in) :: Fx, x

    res = norm_2(x)
    if (res < real_tol) res = 1._rp    
    res = norm_2(fx)/res

  end function residual

end module newton_mod