Merge pull request #26 from jacobwilliams/25-itp-fix

Fix for ITP method
jacobwilliams · Sep 24, 2023 · a811bfd · a811bfd
2 parents 6ee8d35 + 029c19c
commit a811bfd
Show file tree

Hide file tree

Showing 4 changed files with 54 additions and 32 deletions.
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -1,9 +1,11 @@
 name: CI
-on: [push, pull_request]
+on: [push]
 jobs:
 
   Build:
     runs-on: ${{ matrix.os }}
+    permissions:
+      contents: write
     strategy:
       fail-fast: false
       matrix:
@@ -16,20 +18,20 @@ jobs:
 
     steps:
     - name: Checkout code
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
       with:
         submodules: recursive
 
     - name: Install Python
-      uses: actions/setup-python@v1 # Use pip to install latest CMake, & FORD/Jin2For, etc.
+      uses: actions/setup-python@v4 # Use pip to install latest CMake, & FORD/Jin2For, etc.
       with:
         python-version: ${{ matrix.python-version }}
 
     - name: Setup Graphviz
       uses: ts-graphviz/setup-graphviz@v1
 
     - name: Setup Fortran Package Manager
-      uses: fortran-lang/setup-fpm@v4
+      uses: fortran-lang/setup-fpm@v5
       with:
         github-token: ${{ secrets.GITHUB_TOKEN }}
 
@@ -44,6 +46,9 @@ jobs:
       if: contains( matrix.os, 'ubuntu')
       run: |
         sudo apt-get install lcov
+        sudo add-apt-repository ppa:ubuntu-toolchain-r/test
+        sudo apt-get update
+        sudo apt-get install -y gcc-${{ matrix.gcc_v }} gfortran-${{ matrix.gcc_v }}
         sudo update-alternatives \
            --install /usr/bin/gcc gcc /usr/bin/gcc-${{ matrix.gcc_v }} 100 \
            --slave /usr/bin/gfortran gfortran /usr/bin/gfortran-${{ matrix.gcc_v }} \
@@ -58,14 +63,15 @@ jobs:
     - name: Create coverage report
       run: |
         mkdir -p ${{ env.COV_DIR }}
-        lcov --capture --initial --base-directory . --directory build/gfortran_*/ --output-file ${{ env.COV_DIR }}/coverage.base
-        lcov --capture --base-directory . --directory build/gfortran_*/ --output-file ${{ env.COV_DIR }}/coverage.capture
+        mv ./build/gfortran_*/*/* ${{ env.COV_DIR }}
+        lcov --capture --initial --base-directory . --directory ${{ env.COV_DIR }} --output-file ${{ env.COV_DIR }}/coverage.base
+        lcov --capture           --base-directory . --directory ${{ env.COV_DIR }} --output-file ${{ env.COV_DIR }}/coverage.capture
         lcov --add-tracefile ${{ env.COV_DIR }}/coverage.base --add-tracefile ${{ env.COV_DIR }}/coverage.capture --output-file ${{ env.COV_DIR }}/coverage.info
       env:
         COV_DIR: build/coverage
 
     - name: Upload coverage report
-      uses: codecov/codecov-action@v2
+      uses: codecov/codecov-action@v3
       with:
         files: build/coverage/coverage.info
 
@@ -74,7 +80,7 @@ jobs:
 
     - name: Deploy Documentation
       if: github.ref == 'refs/heads/master'
-      uses: JamesIves/[email protected].0
+      uses: JamesIves/github-pages-deploy-action@v4.4.1
       with:
         branch: gh-pages # The branch the action should deploy to.
         folder: doc  # The folder the action should deploy.
diff --git a/README.md b/README.md
@@ -58,28 +58,28 @@ To generate the documentation using [ford](https://github.com/Fortran-FOSS-Progr
 
 ### Methods
 
-The module contains the following methods (in alphabetical order):
+The module contains the following methods:
 
 Procedure | Description | Reference
 --- | --- | ---
-[`anderson_bjorck`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck.html)           | Anderson-Bjorck method | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
-[`anderson_bjorck_king`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck_king.html) | a [variant](https://link.springer.com/content/pdf/bbm%3A978-3-642-05175-3%2F1.pdf) of `anderson_bjorck` | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
-[`barycentric`](https://jacobwilliams.github.io/roots-fortran/proc/barycentric.html)                   | Barycentric interpolation method | [Mendez & Castillo (2021)](https://www.researchgate.net/publication/352162661_A_highly_efficient_numerical_method_to_solve_non-linear_functions_using_barycentric_interpolation)
-[`bdqrf`](https://jacobwilliams.github.io/roots-fortran/proc/bdqrf.html)                               | Bisected Direct Quadratic Regula Falsi | [Gottlieb & Thompson (2010)](http://www.m-hikari.com/ams/ams-2010/ams-13-16-2010/gottliebAMS13-16-2010.pdf)
 [`bisection`](https://jacobwilliams.github.io/roots-fortran/proc/bisection.html)                       | Classic bisection method | Bolzano (1817)
-[`blendtf`](https://jacobwilliams.github.io/roots-fortran/proc/blendtf.html)                           | Blended method of trisection and false position | [Badr, Almotairi, Ghamry (2021)](https://www.mdpi.com/2227-7390/9/11/1306/htm)
+[`regula_falsi`](https://jacobwilliams.github.io/roots-fortran/proc/regula_falsi.html)                 | Classic regula falsi method | ?
+[`muller`](https://jacobwilliams.github.io/roots-fortran/proc/muller.html)                             | Improved Muller method (for real roots only) | [Muller (1956)](https://www.ams.org/journals/mcom/1956-10-056/S0025-5718-1956-0083822-0/S0025-5718-1956-0083822-0.pdf)
 [`brent`](https://jacobwilliams.github.io/roots-fortran/proc/brent.html)                               | Classic Brent's method (a.k.a. Zeroin) | [Brent (1971)](https://maths-people.anu.edu.au/~brent/pd/rpb005.pdf)
 [`brenth`](https://jacobwilliams.github.io/roots-fortran/proc/brenth.html)                             | SciPy variant of `brent` | [SciPy](https://github.com/scipy/scipy/)
 [`brentq`](https://jacobwilliams.github.io/roots-fortran/proc/brentq.html)                             | SciPy variant of `brent` | [SciPy](https://github.com/scipy/scipy/)
-[`chandrupatla`](https://jacobwilliams.github.io/roots-fortran/proc/chandrupatla.html)                 | Hybrid quadratic/bisection algorithm | [Chandrupatla (1997)](https://dl.acm.org/doi/10.1016/S0965-9978%2896%2900051-8)
 [`illinois`](https://jacobwilliams.github.io/roots-fortran/proc/illinois.html)                         | Illinois method | [Dowell & Jarratt (1971)](https://personal.math.ubc.ca/~loew/mech2/Dowell+Jarratt.pdf)
-[`itp`](https://jacobwilliams.github.io/roots-fortran/proc/itp.html)                                   | Interpolate Truncate and Project method | [Oliveira & Takahashi (2020)](https://dl.acm.org/doi/abs/10.1145/3423597)
-[`muller`](https://jacobwilliams.github.io/roots-fortran/proc/muller.html)                             | Improved Muller method (for real roots only) | [Muller (1956)](https://www.ams.org/journals/mcom/1956-10-056/S0025-5718-1956-0083822-0/S0025-5718-1956-0083822-0.pdf)
 [`pegasus`](https://jacobwilliams.github.io/roots-fortran/proc/pegasus.html)                           | Pegasus method | [Dowell & Jarratt (1972)](https://link.springer.com/article/10.1007/BF01932959)
-[`regula_falsi`](https://jacobwilliams.github.io/roots-fortran/proc/regula_falsi.html)                 | Classic regula falsi method | ?
+[`anderson_bjorck`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck.html)           | Anderson-Bjorck method | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
+[`anderson_bjorck_king`](https://jacobwilliams.github.io/roots-fortran/proc/anderson_bjorck_king.html) | a [variant](https://link.springer.com/content/pdf/bbm%3A978-3-642-05175-3%2F1.pdf) of `anderson_bjorck` | [King (1973)](https://link.springer.com/article/10.1007/BF01933405)
 [`ridders`](https://jacobwilliams.github.io/roots-fortran/proc/ridders.html)                           | Classic Ridders method | [Ridders (1979)](https://cs.fit.edu/~dmitra/SciComp/Resources/RidderMethod.pdf)
 [`toms748`](https://jacobwilliams.github.io/roots-fortran/proc/toms748.html)                           | Algorithm 748 |  [Alefeld, Potra, Shi (1995)](https://dl.acm.org/doi/abs/10.1145/210089.210111)
+[`chandrupatla`](https://jacobwilliams.github.io/roots-fortran/proc/chandrupatla.html)                 | Hybrid quadratic/bisection algorithm | [Chandrupatla (1997)](https://dl.acm.org/doi/10.1016/S0965-9978%2896%2900051-8)
+[`bdqrf`](https://jacobwilliams.github.io/roots-fortran/proc/bdqrf.html)                               | Bisected Direct Quadratic Regula Falsi | [Gottlieb & Thompson (2010)](http://www.m-hikari.com/ams/ams-2010/ams-13-16-2010/gottliebAMS13-16-2010.pdf)
 [`zhang`](https://jacobwilliams.github.io/roots-fortran/proc/zhang.html)                               | Zhang's method (with corrections from Stage) | [Zhang (2011)](https://www.cscjournals.org/download/issuearchive/IJEA/Volume2/IJEA_V2_I1.pdf)
+[`itp`](https://jacobwilliams.github.io/roots-fortran/proc/itp.html)                                   | Interpolate Truncate and Project method | [Oliveira & Takahashi (2020)](https://dl.acm.org/doi/abs/10.1145/3423597)
+[`barycentric`](https://jacobwilliams.github.io/roots-fortran/proc/barycentric.html)                   | Barycentric interpolation method | [Mendez & Castillo (2021)](https://www.researchgate.net/publication/352162661_A_highly_efficient_numerical_method_to_solve_non-linear_functions_using_barycentric_interpolation)
+[`blendtf`](https://jacobwilliams.github.io/roots-fortran/proc/blendtf.html)                           | Blended method of trisection and false position | [Badr, Almotairi, Ghamry (2021)](https://www.mdpi.com/2227-7390/9/11/1306/htm)
 
 In general, all the methods are guaranteed to converge. Some will be more efficient (in terms of number of function evaluations) than others for various problems. The methods can be broadly classified into three groups:
 

diff --git a/src/root_module.F90 b/src/root_module.F90
@@ -1237,7 +1237,7 @@ end subroutine bdqrf
 !  * D. E. Muller, "A Method for Solving Algebraic Equations Using an Automatic Computer",
 !    Mathematical Tables and Other Aids to Computation, 10 (1956), 208-215.
 !    [link](https://www.ams.org/journals/mcom/1956-10-056/S0025-5718-1956-0083822-0/S0025-5718-1956-0083822-0.pdf)
-!  * [Roots.jl](https://github.com/JuliaMath/Roots.jl/blob/master/src/simple.jl), 
+!  * [Roots.jl](https://github.com/JuliaMath/Roots.jl/blob/master/src/simple.jl),
 !    Julia version of standard Muller
 
     subroutine muller (me,ax,bx,fax,fbx,xzero,fzero,iflag)
@@ -2517,24 +2517,25 @@ subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)
     integer :: j !! iteration counter
     logical :: root_found !! convergence in x
     logical :: fail !! if we can't interpolate
+    real(wp) :: factor !! factor to ensure f(ax)<f(bx)
 
     real(wp),parameter :: log2 = log(2.0_wp)
 
     ! initialize:
-    if (fax < fbx) then
-        a  = ax
-        b  = bx
-        ya = fax
-        yb = fbx
+    if (fax < fbx) then  ! a requirement for the algorithm ?
+        factor = 1.0_wp
     else
-        a  = bx
-        b  = ax
-        ya = fbx
-        yb = fax
+        ! flip the sign of the function to ensure fax < fbx
+        ! note that a < b is already enforced by the calling routine.
+        factor = -1.0_wp
     end if
+    a  = ax
+    b  = bx
+    ya = factor*fax
+    yb = factor*fbx
     iflag = 0
     term = (b-a)/(2.0_wp*me%rtol)
-    n12 = ceiling ( log(term) / log2 ) ! ceiling(log2(term))
+    n12 = ceiling ( log(abs(term)) / log2 ) ! ceiling(log2(term))
     nmax = n12 + me%n0
     aprev = huge(1.0_wp) ! initialize to unusual values
     bprev = huge(1.0_wp) !
@@ -2571,7 +2572,7 @@ subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)
             end if
 
             ! updating interval:
-            yitp = me%f(xitp)
+            yitp = factor*me%f(xitp)
             if (me%solution(xitp,yitp,xzero,fzero)) return
             if (yitp > 0.0_wp) then
                 b = xitp
@@ -2591,7 +2592,7 @@ subroutine itp(me,ax,bx,fax,fbx,xzero,fzero,iflag)
             ! [this can happen in the test cases].
             ! So just do a bisection
             xitp = bisect(a,b)
-            yitp = me%f(xitp)
+            yitp = factor*me%f(xitp)
             if (me%solution(xitp,yitp,xzero,fzero)) return
             if (ya*yitp<0.0_wp) then
                 ! root lies between a and xitp

diff --git a/test/root_tests.f90 b/test/root_tests.f90
@@ -1300,12 +1300,27 @@ subroutine problems(x, ax, bx, fx, xroot, cases, num_of_problems, latex)
         if (present(x)) f = x
         if (present(latex)) latex = 'x'
 
+    ! another linear test case (see Issue #25)
+    case(133)
+
+        tmp : block
+            real(wp),parameter :: x1 = 0.0_wp
+            real(wp),parameter :: x2 = 157.08_wp
+            real(wp),parameter :: y1 = 1012.0_wp
+            real(wp),parameter :: y2 = -1114.0_wp
+            a = x1
+            b = x2
+            root = 7.4771853245531515E+01_wp
+            if (present(x)) f = y1 + ((x - x1) / (x2 - x1)) * (y2 - y1)
+            if (present(latex)) latex = '1012 - 2126 x / 157.08'
+        end block tmp
+
     case default
         write(*,*) 'invalid case: ', nprob
         error stop 'invalid case'
     end select
 
-    if (present(num_of_problems)) num_of_problems = 132
+    if (present(num_of_problems)) num_of_problems = 133
 
     ! outputs:
     if (present(ax))    ax = a