added an optional bisect_used output

to check if the method failed and bisection was used afterwards.
LaTeX table updates

Author: jacobwilliams

src/root_module.F90

@@ -389,7 +389,7 @@ end function root_name_to_method
 
     subroutine root_scalar_by_name(method,fun,ax,bx,xzero,fzero,iflag,&
                                    ftol,rtol,atol,maxiter,fax,fbx,&
-                                   bisect_on_failure)
+                                   bisect_on_failure,bisect_used)
 
     implicit none
 
@@ -410,6 +410,7 @@ subroutine root_scalar_by_name(method,fun,ax,bx,xzero,fzero,iflag,&
                                                         !! it will be retried using the bisection method.
                                                         !! (default is False). Note that this can use up
                                                         !! to `maxiter` additional function evaluations.
+    logical,intent(out),optional  :: bisect_used  !! if the bisection method was used after failure
 
     type(root_method) :: r
 
@@ -418,7 +419,7 @@ subroutine root_scalar_by_name(method,fun,ax,bx,xzero,fzero,iflag,&
     if (r%id /= 0) then
         call root_scalar(r,fun,ax,bx,xzero,fzero,iflag,&
                          ftol,rtol,atol,maxiter,fax,fbx,&
-                         bisect_on_failure)
+                         bisect_on_failure,bisect_used)
     else
         iflag = -999    ! invalid method
         return
@@ -433,7 +434,7 @@ end subroutine root_scalar_by_name
 
     subroutine root_scalar_by_type(method,fun,ax,bx,xzero,fzero,iflag,&
                                    ftol,rtol,atol,maxiter,fax,fbx,&
-                                   bisect_on_failure)
+                                   bisect_on_failure,bisect_used)
 
     implicit none
 
@@ -454,6 +455,7 @@ subroutine root_scalar_by_type(method,fun,ax,bx,xzero,fzero,iflag,&
                                                         !! it will be retried using the bisection method.
                                                         !! (default is False). Note that this can use up
                                                         !! to `maxiter` additional function evaluations.
+    logical,intent(out),optional  :: bisect_used  !! if the bisection method was used after failure
 
     class(root_solver),allocatable :: s
 
@@ -486,7 +488,8 @@ subroutine root_scalar_by_type(method,fun,ax,bx,xzero,fzero,iflag,&
     end select
 
     call s%initialize(func_wrapper,ftol,rtol,atol,maxiter)
-    call s%solve(ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure)
+    call s%solve(ax,bx,xzero,fzero,iflag,fax,fbx,&
+                 bisect_on_failure,bisect_used)
 
     contains
 
@@ -505,7 +508,7 @@ end subroutine root_scalar_by_type
 !>
 !  Main wrapper routine for all the methods.
 
-    subroutine solve(me,ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure)
+    subroutine solve(me,ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure,bisect_used)
 
     implicit none
 
@@ -521,10 +524,13 @@ subroutine solve(me,ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure)
                                                           !! it will be retried using the bisection method.
                                                           !! (default is False). Note that this can use up
                                                         !! to `maxiter` additional function evaluations.
+    logical,intent(out),optional     :: bisect_used !! if the bisection method was used after failure
 
     real(wp) :: fa !! `f(ax)` passed to the lower level routine
     real(wp) :: fb !! `f(bx)` passed to the lower level routine
 
+    if (present(bisect_used)) bisect_used = .false.
+
     if (ax==bx) then
         ! ax must be /= bx
         iflag = -4
@@ -563,13 +569,19 @@ subroutine solve(me,ax,bx,xzero,fzero,iflag,fax,fbx,bisect_on_failure)
             ! if it failed, then we have the option to then try bisection
             if (iflag /= 0) then
                 if (present(bisect_on_failure)) then
-                    if (bisect_on_failure) then
-                        ! use the wrapper routine for that with the input class
-                        call root_scalar(root_method_bisection,&
-                                         func_wrapper,ax,bx,xzero,fzero,iflag,&
-                                         me%ftol,me%rtol,me%atol,me%maxiter,fa,fb,&
-                                         bisect_on_failure = .false.)
-                    end if
+                    select type (me)
+                    class is (bisection_solver)
+                        ! can't bisect after failed bisection
+                    class default
+                        if (bisect_on_failure) then
+                            ! use the wrapper routine for that with the input class
+                            call root_scalar(root_method_bisection,&
+                                             func_wrapper,ax,bx,xzero,fzero,iflag,&
+                                             me%ftol,me%rtol,me%atol,me%maxiter,fa,fb,&
+                                             bisect_on_failure = .false.)
+                            if (present(bisect_used)) bisect_used = .true.
+                        end if
+                    end select
                 end if
             end if
 

test/root_tests.f90

@@ -68,7 +68,7 @@ program root_tests
     case(real32)
         atol = 1.0e-5_wp
         rtol = 1.0e-5_wp
-        ftol = 1.0e-5_wp
+        ftol = 1.0e-6_wp
         tol_for_check = 1.0e-4_wp
         kind_str = 'real32'
     case(real64)
@@ -171,6 +171,11 @@ program root_tests
     end do
 
     write(iunit_results,'(a)') '\end{longtable}'
+
+    write(iunit_results,'(a)') ' '
+    write(iunit_results,'(a)') '$^*$ indicates that the method failed and bisection was used.'
+    write(iunit_results,'(a)') ' '
+
     write(iunit_results,'(a)') '\end{landscape}'
     write(iunit_results,'(a)') '\end{document}'
 
@@ -213,6 +218,7 @@ subroutine test()
         character(len=1000) :: line
         integer,dimension(number_of_methods) :: fevals !! func evals for each case
         real(wp),dimension(number_of_methods) :: cases_per_sec !! speed of each case
+        logical,dimension(number_of_methods) :: bisect_used !! if the method failed and bisection was used
         integer :: best_feval, worst_feval
         real(wp) :: best_cases_per_second
         character(len=:),allocatable :: best,failures
@@ -224,9 +230,11 @@ subroutine test()
         character(len=10) :: itmp, nstr
         integer :: n_bounds_cases, bounds_cases_first
         character(len=30) :: root_str
+        logical :: bis
 
         integer,parameter :: n_repeat = 1  !! number of times to repeat each test for timing purposes
         integer,parameter :: maxiter = 1000 !! maximum number of iterations
+        logical,parameter :: bisect_on_failure = .true. !! to fall back to bisection if method fails
 
         write(output_unit,fmt) &
             repeat('-',20),repeat('-',3),repeat('-',4),repeat('-',16),&
@@ -252,13 +260,16 @@ subroutine test()
             call cpu_time(tstart)
             do irepeat = 1, n_repeat
                 ifunc = 0 ! reset func evals counter
-                call root_scalar(methods(imeth),func,ax,bx,xzero,fzero,iflag,bisect_on_failure=.true., &
+                call root_scalar(methods(imeth),func,ax,bx,xzero,fzero,iflag,&
+                                 bisect_on_failure = bisect_on_failure, &
+                                 bisect_used = bis, &
                                  atol = atol, rtol = rtol, ftol = ftol, maxiter = maxiter )
             end do
             call cpu_time(tfinish)
 
             error = xzero-root
-            write(line, dfmt) trim(methods(imeth)),nprob,n,error,xzero,fzero,ifunc,iflag,n_repeat / (tfinish-tstart), trim(tmp)
+            write(line, dfmt) trim(methods(imeth)),nprob,n,error,xzero,fzero,ifunc,iflag,&
+                              n_repeat / (tfinish-tstart), trim(tmp)
 
             root_found = abs(fzero) <= tol_for_check !.and. abs(error) <= tol_for_check
 
@@ -276,6 +287,7 @@ subroutine test()
                 fevals(imeth) = huge(1)
                 cases_per_sec(imeth) = huge(1.0_wp)
             end if
+            bisect_used(imeth) = bis
 
         end do
 
@@ -314,7 +326,10 @@ subroutine test()
         fevals_line = ''
         do i = 1, size(fevals)
             write(itmp, '(I10)') fevals(i); itmp = trim(adjustl(itmp))
-            if (fevals(i) == best_feval) then
+            if (bisect_used(i)) itmp = trim(itmp) // '$^*$'
+            if (bisect_used(i)) then
+                fevals_line = fevals_line // '\cellcolor{red!75} ' // itmp // '&'  ! always color red if it failed
+            else if (fevals(i) == best_feval) then
                 if (count(fevals==best_feval)==1) then
                     ! the sole winner
                     fevals_line = fevals_line // '\cellcolor{green!75} ' // itmp // '&'