rootstest.f90

Go to the documentation of this file.

!#############################################################################
!#                                                                           #
!# fosite - 3D hydrodynamical simulation program                             #
!# module: rootstest.f90                                                     #
!#                                                                           #
!# Copyright (C) 2016                                                        #
!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de>                       #
!#                                                                           #
!# 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 of the License, 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, GOOD TITLE or        #
!# NON INFRINGEMENT.  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., 675 Mass Ave, Cambridge, MA 02139, USA.                 #
!#                                                                           #
!#############################################################################

!----------------------------------------------------------------------------!
!----------------------------------------------------------------------------!
PROGRAM rootstest
  USE roots
#include "tap.h"
  IMPLICIT NONE
  !--------------------------------------------------------------------------!
  INTERFACE
    PURE SUBROUTINE func(x,fx,plist)
      IMPLICIT NONE
      REAL, INTENT(IN)  :: x
      REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
      REAL, INTENT(OUT) :: fx
    END SUBROUTINE func
  END INTERFACE
  INTERFACE
    PURE SUBROUTINE funcd(x,fx,dfx,plist)
      IMPLICIT NONE
      REAL, INTENT(IN)  :: x
      REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
      REAL, INTENT(OUT) :: fx,dfx
    END SUBROUTINE funcd
  END INTERFACE
  COMMON /test_num/ k
  !--------------------------------------------------------------------------!
  INTEGER, PARAMETER :: num_tests = 13
  INTEGER, PARAMETER :: num_methods = 8
  CHARACTER(LEN=16), PARAMETER, DIMENSION(NUM_METHODS) :: method_name = (/ &
                                                       "       Bisection", &
                                                       "    Regula Falsi", &
                                                       "         Pegasus", &
                                                       "            King", &
                                                       " Anderson Bjoerk", &
                                                       "          Ridder", &
                                                       "    Brent Dekker", &
                                                       "          Newton" /)
  REAL, PARAMETER, DIMENSION(2,NUM_TESTS) :: bounds = reshape((/ 0.0, 1.2, &
                                                                 0.4, 1.6, &
                                                                -0.5, 1.9, &
                                                                -0.5, 0.7, &
                                                                -1.4, 1.0, &
                                                                -0.8, 1.6, &
                                                                -0.5, 1.9, &
                                                               0.001, 1.201, &
                                                                -0.9, 1.5, &
                                                                 0.4, 1.0, &
                                                                -1.2, 0.0, &
                                                                 1.0, 3.4, &
                                                                 0.0, 5.0 /), &
                                             (/2,num_tests/))
  DOUBLE PRECISION, PARAMETER, DIMENSION(NUM_TESTS) :: ref_roots = &
     (/ 3.9942229171096819451d-01, 8.0413309750366432374d-01, 9.0340766319186021294d-01, &
      7.7014241346192677110d-02, 2.5920449372984746773d-01, 5.3674166257799978186d-01, &
      4.4754176206055907112d-01, 1.1111111111111111111d-01, 5.0000003403025908310d-01, &
      6.7980892150470050192d-01, -3.5938136638046273022d-01, 1.6487212707001281468d-00, &
      1.0000000000000000000d-00 &
     /)
  !--------------------------------------------------------------------------!
  REAL              :: root, xm, dx_rel, dx_acc, plist(1)
  INTEGER           :: i,k, iter, error
  CHARACTER(LEN=64) :: tap_message
  LOGICAL           :: verbose_results = .false.
  !--------------------------------------------------------------------------!
  ! some information
  WRITE (*,"(A)") "========================================================="
  WRITE (*,"(A)") "          << Testing root finding algorithms >>          "
  WRITE (*,"(A)") "========================================================="
  
  
  tap_plan((num_methods-1)*num_tests)

  DO k=1,num_tests
    ! plist is an array of real parameters that can be passed to the function func
    ! so far just a dummy, but it can be used if necessary;
    plist(1) = 1.0
    ! reset error code; if GetRoot returns error .ne. 0 something bad happend
    error = 0
    ! initial guess: arithmetic mean of boundary values
    xm = 0.5*(bounds(1,k)+bounds(2,k))
    print '(A,I2,A,2(ES10.2))', "Test #",k,"  Search Interval: ",bounds(1,k),bounds(2,k)
    print '(A)', "---------------------------------------------------------"
    IF (verbose_results) print '(A20,ES27.19)', "reference result", ref_roots(k)
    DO i=1,num_methods
      SELECT CASE(i)
      CASE(1)
        CALL getroot_bisection(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(2)
        ! don't use Regula Falsi because of slow convergence in some cases
        CALL getroot_regulafalsi(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(3)
        CALL getroot_pegasus(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(4)
        CALL getroot_king(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(5)
        CALL getroot_andersonbjoerk(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(6)
        CALL getroot_ridder(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(7)
        CALL getroot_brentdekker(func,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      CASE(8)
        CALL getroot_newton(funcd,bounds(1,k),bounds(2,k),root,error,plist,xm,iter)
      END SELECT

      ! compute the relative  error using reference root
      dx_rel = abs(1.0-root/ref_roots(k))
      IF (verbose_results) CALL printresults(method_name(i),root,ref_roots(k),error,iter)
      ! exclude Regular Falsi from check
      IF (i.NE.2) THEN
      
        dx_acc = default_accuracy
        IF (k.EQ.12) THEN
          ! this test never yields the required accuracy, see \cite engeln2011
          dx_acc = 1.0e-5
        END IF
        WRITE (tap_message,'(A16,A11,ES9.2,A3,ES9.2)') method_name(i), ": dx_rel = ", dx_rel, " < ", dx_acc
        tap_check(dx_rel.LE.dx_acc,tap_message)
        
      END IF
    END DO
    print '(A)', "---------------------------------------------------------"
  END DO

  tap_done

CONTAINS


  SUBROUTINE printresults(method,root,ref_root,error,iter)
    IMPLICIT NONE
    CHARACTER(LEN=*), INTENT(IN) :: method
    REAL, INTENT(IN) :: root
    DOUBLE PRECISION, INTENT(IN) :: ref_root
    INTEGER, INTENT(IN) :: error,iter
    IF (error.NE.0) THEN
        print *,trim(geterrormessage(error))
    ELSE
        print '(A20,ES23.15,ES9.1,I4)',trim(method),root,abs(1.0-root/ref_root),iter
    END IF
  END SUBROUTINE printresults


!   TAP_PLAN(7)
! 
!   
! 
!   ! Check if the random numbers are in (0,1)
!   ! and if the average is near 0.5
!   TAP_CHECK_CLOSE(r/imax,0.5,1.E-4,"Average close to 0.5.")
!   TAP_CHECK_GE(rmin,0.,"All are bigger (or equal) than 0.")
!   TAP_CHECK_LE(root,1.,"All are smaller (or equal) than 1.")
!   TAP_CHECK_CLOSE(rmin,0.,1.E-4,"Lower limit is close to 0.")
!   TAP_CHECK_CLOSE(rmax,1.,1.E-4,"Upper limit is close to 1.")
! 
!   TAP_CHECK(x.EQ.x0,"SuperKiss64")
! 

END PROGRAM rootstest


! function definitions
PURE SUBROUTINE func(x,fx,plist)
IMPLICIT NONE
  COMMON /test_num/ k
  REAL, INTENT(IN)  :: x
  REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
  REAL, INTENT(OUT) :: fx
  INTEGER :: k
  SELECT CASE(k)
  CASE(1)
    fx = x*x*(x*x/3. + sqrt(2.)*sin(x)) - sqrt(3.)/18.
  CASE(2)
    fx = 11*x**11 - 1.0
  CASE(3)
    fx = 35*x**35 - 1.0
  CASE(4)
    fx = 2*(x*exp(-9.)-exp(-9*x)) + 1.
  CASE(5)
    fx = x*x - (1.-x)**9
  CASE(6)
    fx = (x-1.)*exp(-9*x) + x**9
  CASE(7)
    fx = x*x + sin(x/9.) - 0.25
  CASE(8)
    fx = 0.125*(9.-1./x)
  CASE(9)
    fx = tan(x) - x - 0.0463025
  CASE(10)
    fx = x*(x+sin(sqrt(75.)*x)) - 0.2
  CASE(11)
    fx = x**9 + 1e-4
  CASE(12)
    fx = log(x) + 0.5*x*x/exp(1.) -2*x/sqrt(exp(1.)) + 1.
  CASE DEFAULT
    fx = x - 1.0
  END SELECT
END SUBROUTINE func


PURE SUBROUTINE funcd(x,fx,dfx,plist)
  IMPLICIT NONE
  COMMON /test_num/ k
  REAL, INTENT(IN)  :: x
  REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
  REAL, INTENT(OUT) :: fx,dfx
  INTEGER :: k
  SELECT CASE(k)
  CASE(1)
    fx = x*x*(x*x/3. + sqrt(2.)*sin(x)) - sqrt(3.)/18.
    dfx= 2*x*(sqrt(2.)*sin(x)+x*x/3.)+x*x*(sqrt(2.)*cos(x)+(2*x)/3.)
  CASE(2)
    fx = 11*x**11 - 1.0
    dfx = 121*x**10
  CASE(3)
    fx = 35*x**35 - 1.0
    dfx= 35*35*x**34 - 1.0
  CASE(4)
    fx = 2*(x*exp(-9.)-exp(-9*x)) + 1.
    dfx= 2*(9*exp(-9*x)+exp(-9.))
  CASE(5)
    fx = x*x - (1.-x)**9
    dfx= 2*x + 9*(1.-x)**8
  CASE(6)
    fx = (x-1.)*exp(-9*x) + x**9
    dfx= (10.-9*x)*exp(-9*x) + 9*x**8
  CASE(7)
    fx = x*x + sin(x/9.) - 0.25
    dfx= 2*x + cos(x/9.)/9.
  CASE(8)
    fx = 0.125*(9.-1./x)
    dfx= 0.125/x**2
  CASE(9)
    fx = tan(x) - x - 0.0463025
    dfx= 1./cos(x)**2 - 1.
  CASE(10)
    fx = x*(x+sin(sqrt(75.)*x)) - 0.2
    dfx= sin(sqrt(75.)*x) + x*(2. + sqrt(75.)*cos(sqrt(75.)*x))
  CASE(11)
    fx = x**9 + 1e-4
    dfx= 9*x**8
  CASE(12)
    fx = log(x) + 0.5*x*x/exp(1.) - 2*x/sqrt(exp(1.)) + 1.
    dfx= 1./x + x/exp(1.) - 2./sqrt(exp(1.))
  CASE DEFAULT
    fx = x - 1.0
    dfx= 1.
  END SELECT
END SUBROUTINE funcd