root_solver.F90

Go to the documentation of this file.

!! Copyright (C) 2005-2006 Heiko Appel
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module root_solver_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                               &
    root_solver_t,                        &
    root_solver_init,                     &
    root_solver_read,                     &
    droot_solver_run
 
  integer, public, parameter ::           &
    ROOT_NEWTON      =  3
 
  type root_solver_t
    private
    integer :: solver_type
    integer :: dim
    integer :: maxiter
    integer :: usediter
    real(real64)   :: abs_tolerance
    real(real64)   :: rel_tolerance
  end type root_solver_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine root_solver_init(rs, namespace, dimensionality, solver_type, maxiter, &
    rel_tolerance, abs_tolerance)
    type(root_solver_t), intent(out) :: rs
    type(namespace_t),   intent(in)  :: namespace
    integer,             intent(in)  :: dimensionality
    integer, optional,   intent(in)  :: solver_type, maxiter
    real(real64), optional,     intent(in)  :: rel_tolerance, abs_tolerance
 
    ! no push_sub, called too often
 
    ! Set dimension
    rs%dim = dimensionality
 
    ! Read config parameters for the solver
    call root_solver_read(rs, namespace)
 
    if (present(solver_type))     rs%solver_type     = solver_type
    if (present(maxiter))         rs%maxiter         = maxiter
    if (present(rel_tolerance))   rs%rel_tolerance   = rel_tolerance
    if (present(abs_tolerance))   rs%abs_tolerance   = abs_tolerance
 
  end subroutine root_solver_init
 
 
  ! ---------------------------------------------------------
  subroutine root_solver_read(rs, namespace)
    type(root_solver_t),   intent(inout) :: rs
    type(namespace_t),     intent(in)  :: namespace
 
    push_sub(root_solver_read)
 
    !%Variable RootSolver
    !%Type integer
    !%Default root_newton
    !%Section Math::RootSolver
    !%Description
    !% Specifies what kind of root solver will be used.
    !%Option root_newton 3
    !% Newton-Raphson solver with backtracking.
    !%End
    call parse_variable(namespace, 'RootSolver', root_newton, rs%solver_type)
    if (.not. varinfo_valid_option('RootSolver', rs%solver_type)) then
      call messages_input_error(namespace, 'RootSolver')
    end if
 
    !%Variable RootSolverMaxIter
    !%Type integer
    !%Default 500
    !%Section Math::RootSolver
    !%Description
    !% In case of an iterative root solver, this variable determines the maximum number
    !% of iteration steps.
    !%End
    call parse_variable(namespace, 'RootSolverMaxIter', 500, rs%maxiter)
 
    !%Variable RootSolverAbsTolerance
    !%Type float
    !%Default 1e-10
    !%Section Math::RootSolver
    !%Description
    !% Relative tolerance for the root-finding process.
    !%End
    call parse_variable(namespace, 'RootSolverAbsTolerance', 1e-10_real64, rs%abs_tolerance)
 
    pop_sub(root_solver_read)
  end subroutine root_solver_read
 
  ! ---------------------------------------------------------
  subroutine droot_solver_run(rs, func, root, success, startval)
    type(root_solver_t),    intent(in)  :: rs
    real(real64),           intent(out) :: root(:)
    logical,                intent(out) :: success
    real(real64), optional, intent(in)  :: startval(:)
    interface
      subroutine func(z, f, jf)
        import real64
        implicit none
        real(real64), intent(in)  :: z(:)
        real(real64), intent(out) :: f(:), jf(:, :)
      end subroutine func
    end interface
 
    ! no push_sub, called too often
 
    ! Initializations
    root = m_zero
    success = .false.
 
    select case (rs%solver_type)
    case (root_newton)
      call droot_newton(rs, func, root, startval, success)
    case default
      write(message(1), '(a,i4,a)') "Error in droot_solver_run: '", rs%solver_type, "' is not a valid root solver"
      call messages_fatal(1)
    end select
 
  end subroutine
 
  ! ---------------------------------------------------------
  subroutine droot_newton(rs, func, root, startval, success)
    type(root_solver_t),    intent(in)  :: rs
    real(real64),           intent(out) :: root(:)
    real(real64), optional, intent(in)  :: startval(:)
    logical,                intent(out) :: success
    interface
      subroutine func(z, f, jf)
        import real64
        implicit none
        real(real64), intent(in)  :: z(:)
        real(real64), intent(out) :: f(:), jf(:, :)
      end subroutine func
    end interface
 
    integer :: iter
    real(real64)   :: err, prev_err
    real(real64), allocatable :: f(:), jf(:, :), delta(:, :), rhs(:, :), root_trial(:)
 
    real(real64) :: mm, prev_m, alpha, err_trial
 
    ! no push_sub, called too often
 
    safe_allocate(f(1:rs%dim))
    safe_allocate(root_trial(1:rs%dim))
    safe_allocate(jf(1:rs%dim, 1:rs%dim))
    safe_allocate(delta(1:rs%dim, 1))
    safe_allocate(rhs(1:rs%dim, 1))
 
    if (present(startval)) then
      root = startval
    else
      root = m_zero
    end if
 
    call func(root, f, jf)
    err = norm2(f)
    prev_err = err
    mm = m_one
    prev_m = m_one
 
    success = .true.
    iter = 0
    do while(err > rs%abs_tolerance)
      rhs(1:rs%dim, 1) = -mm * f(1:rs%dim)
 
      call lalg_linsyssolve(rs%dim, 1, jf, rhs, delta)
 
      ! Backtracking damping
      alpha = 1.0_real64
      do
        root_trial = root + alpha * delta(:,1)
        call func(root_trial, f, jf)
        err_trial = norm2(f)
 
        if (err_trial < err) exit   ! accept
        alpha = 0.5_real64 * alpha
        if (alpha < 1e-12_real64) then
          success = .false.
          exit
        end if
      end do
      if (.not. success) exit ! Damping failed
 
      root = root_trial
      prev_err = err
      err = err_trial
 
      iter = iter + 1
      if (iter > rs%maxiter) then
        success = .false.
        exit
      end if
 
      ! Let us assume we have a linear convergence due to a multiplicity m of the solution,
      ! Then  we have |err_i+1| = (m-1)/m |err_i|
      ! Or (m-1)/m = |err_i+1|/|err_i| = A, where A is the asymptotic error constant
      ! for i sufficiently large.
      !
      ! Then we can estimate m from the first two iterations and use it to accelerate the convergence
      !
      ! See https://web.archive.org/web/20190524083302/http://mathfaculty.fullerton.edu/mathews/n2003/NewtonAccelerateMod.html
      ! as well as
      ! https://en.wikipedia.org/wiki/Newton%27s_method#Slow_convergence_for_roots_of_multiplicity_greater_than_1
      ! for some more details
      if (iter == 1) then
        prev_m = floor(m_one / (m_one - err/prev_err))
      end if
      if (iter == 2) then
        mm = floor(m_one / (m_one - err/prev_err))
        ! If we get the same multiplicity, we use it for the next iterations
        if (abs(mm - prev_m) > 0.5_real64) mm = m_one
        ! If the multiplicity is going to infinity, we do not want to use this value
        if ( mm > 10._real64 .or. mm < 1.0_real64) mm = m_one
      end if
 
    end do
 
    safe_deallocate_a(f)
    safe_deallocate_a(jf)
    safe_deallocate_a(delta)
    safe_deallocate_a(rhs)
 
  end subroutine droot_newton
 
end module root_solver_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: