Merge pull request #6 from fortran-fans/add_hybrd_interface

Add hybrd interface and example.

Author: zoziha

.github/workflows/fpm.yml

@@ -58,4 +58,5 @@ jobs:
         fpm run --example example_lmder1
         fpm run --example example_lmdif1
         fpm run --example example_primes
+        fpm run --example example_hybrd
         fpm run --example example_hybrd1

examples/CMakeLists.txt

@@ -1,5 +1,8 @@
 include_directories(${PROJECT_BINARY_DIR}/src)
 
+add_executable(example_hybrd example_hybrd.f90)
+target_link_libraries(example_hybrd minpack)
+
 add_executable(example_hybrd1 example_hybrd1.f90)
 target_link_libraries(example_hybrd1 minpack)
 

examples/example_hybrd.f90

@@ -0,0 +1,101 @@
+!> The problem is to determine the values of x(1), x(2), ..., x(9)
+!>  which solve the system of tridiagonal equations
+!>     (3-2*x(1))*x(1)           -2*x(2)                   = -1
+!>             -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
+!>                                 -x(8) + (3-2*x(9))*x(9) = -1
+program example_hybrd
+
+    use minpack, only: hybrd, enorm, dpmpar
+    implicit none
+    integer j, n, maxfev, ml, mu, mode, nprint, info, nfev, ldfjac, lr, nwrite
+    double precision xtol, epsfcn, factor, fnorm
+    double precision x(9), fvec(9), diag(9), fjac(9, 9), r(45), qtf(9), &
+        wa1(9), wa2(9), wa3(9), wa4(9)
+
+    !> Logical output unit is assumed to be number 6.
+    data nwrite/6/
+
+    n = 9
+
+    !> The following starting values provide a rough solution.
+    do j = 1, 9
+        x(j) = -1.0d0
+    end do
+
+    ldfjac = 9
+    lr = 45
+
+    !> Set xtol to the square root of the machine precision.
+    !>  unless high precision solutions are required,
+    !>  this is the recommended setting.
+    xtol = dsqrt(dpmpar(1))
+
+    maxfev = 2000
+    ml = 1
+    mu = 1
+    epsfcn = 0.0d0
+    mode = 2
+    do j = 1, 9
+        diag(j) = 1.0d0
+    end do
+    factor = 1.0d2
+    nprint = 0
+
+    call hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, &
+               mode, factor, nprint, info, nfev, fjac, ldfjac, &
+               r, lr, qtf, wa1, wa2, wa3, wa4)
+    fnorm = enorm(n, fvec)
+    write (nwrite, 1000) fnorm, nfev, info, (x(j), j=1, n)
+
+1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
+           5x, "NUMBER OF FUNCTION EVALUATIONS", i10// &
+           5x, "EXIT PARAMETER", 16x, i10// &
+           5x, "FINAL APPROXIMATE SOLUTION"//(5x, 3d15.7))
+
+    !> Results obtained with different compilers or machines
+    !>  may be slightly different.
+    !>
+    !>> FINAL L2 NORM OF THE RESIDUALS  0.1192636D-07
+    !>>
+    !>> NUMBER OF FUNCTION EVALUATIONS        14
+    !>>
+    !>> EXIT PARAMETER                         1
+    !>>
+    !>> FINAL APPROXIMATE SOLUTION
+    !>>
+    !>>  -0.5706545D+00 -0.6816283D+00 -0.7017325D+00
+    !>>  -0.7042129D+00 -0.7013690D+00 -0.6918656D+00
+    !>>  -0.6657920D+00 -0.5960342D+00 -0.4164121D+00
+
+contains
+
+    !> Subroutine fcn for hybrd example.
+    subroutine fcn(n, x, fvec, iflag)
+
+        implicit none
+        integer n, iflag
+        double precision x(n), fvec(n)
+
+        integer k
+        double precision one, temp, temp1, temp2, three, two, zero
+        data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
+
+        if (iflag /= 0) go to 5
+
+        !! Insert print statements here when nprint is positive.
+
+        return
+5       continue
+        do k = 1, n
+            temp = (three - two*x(k))*x(k)
+            temp1 = zero
+            if (k /= 1) temp1 = x(k - 1)
+            temp2 = zero
+            if (k /= n) temp2 = x(k + 1)
+            fvec(k) = temp - temp1 - two*temp2 + one
+        end do
+        return
+
+    end subroutine fcn
+
+end program example_hybrd

examples/example_hybrd1.f90

@@ -2,17 +2,17 @@
 !>  which solve the system of tridiagonal equations.
 !>
 !>  (3-2*x(1))*x(1)           -2*x(2)                   = -1
-!>  -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
-!>  -x(8) + (3-2*x(9))*x(9) = -1
+!>          -x(i-1) + (3-2*x(i))*x(i)         -2*x(i+1) = -1, i=2-8
+!>                              -x(8) + (3-2*x(9))*x(9) = -1
 program example_hybrd1
 
-    use minpack, only: hybrd1
+    use minpack, only: hybrd1, dpmpar, enorm
     implicit none
     integer j, n, info, lwa, nwrite
     double precision tol, fnorm
     double precision x(9), fvec(9), wa(180)
-    double precision enorm, dpmpar
 
+    !> Logical output unit is assumed to be number 6.
     data nwrite/6/
 
     n = 9
@@ -23,6 +23,10 @@ program example_hybrd1
     end do
 
     lwa = 180
+
+    !> Set tol to the square root of the machine precision.
+    !>  unless high precision solutions are required,
+    !>  this is the recommended setting.
     tol = dsqrt(dpmpar(1))
 
     call hybrd1(fcn, n, x, fvec, tol, info, wa, lwa)
@@ -35,7 +39,7 @@ program example_hybrd1
            (5x, 3d15.7))
 
     !> Results obtained with different compilers or machines
-    !> may be slightly different.
+    !>  may be slightly different.
     !>
     !>> FINAL L2 NORM OF THE RESIDUALS  0.1192636D-07
     !>>
@@ -49,7 +53,7 @@ program example_hybrd1
 
 contains
 
-    !> subroutine fcn for hybrd1 example.
+    !> Subroutine fcn for hybrd1 example.
     subroutine fcn(n, x, fvec, iflag)
 
         implicit none

fpm.toml

@@ -17,6 +17,11 @@ auto-examples = false
 [install]
 library = false
 
+[[ example ]]
+name = "example_hybrd"
+source-dir = "examples"
+main = "example_hybrd.f90"
+
 [[ example ]]
 name = "example_hybrd1"
 source-dir = "examples"

src/minpack.f90

@@ -12,6 +12,26 @@ double precision function enorm(n,x)
     double precision x(n)
     end function
 
+    !> The purpose of `hybrd` is to find a zero of a system of N non-
+    !>  linear functions in N variables by a modification of the Powell
+    !>  hybrid method.  The user must provide a subroutine which calcu-
+    !>  lates the functions.  The Jacobian is then calculated by a for-
+    !>  ward-difference approximation.
+    subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag, &
+                         mode,factor,nprint,info,nfev,fjac,ldfjac, &
+                         r,lr,qtf,wa1,wa2,wa3,wa4)
+        integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr
+        double precision xtol,epsfcn,factor
+        double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),qtf(n), &
+                         wa1(n),wa2(n),wa3(n),wa4(n)
+        interface
+            subroutine fcn(n,x,fvec,iflag)
+                integer n,iflag
+                double precision x(n),fvec(n)
+            end subroutine fcn
+        end interface
+    end subroutine hybrd
+
     !> The purpose of `hybrd1` is to find a zero of a system of
     !>  n nonlinear functions in n variables by a modification
     !>  of the powell hybrid method.