refactor lmder1 example

Author: jacobwilliams

examples/example_lmder1.f90

@@ -1,93 +1,95 @@
-module testmod_der1
-implicit none
-private
-public fcn, dp
-
-integer, parameter :: dp=kind(0d0)
-
-contains
-
-subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag)
-integer, intent(in) :: m, n, ldfjac
-integer, intent(inout) :: iflag
-real(dp), intent(in) :: x(n)
-real(dp), intent(inout) :: fvec(m), fjac(ldfjac, n)
-
-integer :: i
-real(dp) :: tmp1, tmp2, tmp3, tmp4, y(15)
-y = [1.4D-1, 1.8D-1, 2.2D-1, 2.5D-1, 2.9D-1, 3.2D-1, 3.5D-1, 3.9D-1, &
-        3.7D-1, 5.8D-1, 7.3D-1, 9.6D-1, 1.34D0, 2.1D0, 4.39D0]
-
-if (iflag == 1) then
-    do i = 1, 15
-        tmp1 = i
-        tmp2 = 16 - i
-        tmp3 = tmp1
-        if (i > 8) tmp3 = tmp2
-        fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
-    end do
-else
-    do i = 1, 15
-        tmp1 = i
-        tmp2 = 16 - i
-        tmp3 = tmp1
-        if (i > 8) tmp3 = tmp2
-        tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
-        fjac(i,1) = -1.D0
-        fjac(i,2) = tmp1*tmp2/tmp4
-        fjac(i,3) = tmp1*tmp3/tmp4
-    end do
-end if
-end subroutine
-
-end module
+program example_lmder1
 
+use minpack_module, only: wp, enorm, lmder1, chkder
+use iso_fortran_env, only: nwrite => output_unit
 
-program example_lmder1
-use minpack_module, only: enorm, lmder1, chkder
-use testmod_der1, only: dp, fcn
 implicit none
 
+integer, parameter :: n = 3
+integer, parameter :: m = 15
+integer, parameter :: lwa = 5*n+m
+
 integer :: info
-real(dp) :: tol, x(3), fvec(15), fjac(size(fvec), size(x))
-integer :: ipvt(size(x))
-real(dp), allocatable :: wa(:)
+real(wp) :: tol, x(n), fvec(m), fjac(m,n)
+integer :: ipvt(n)
+real(wp) :: wa(lwa)
 
 ! The following starting values provide a rough fit.
-x = [1._dp, 1._dp, 1._dp]
+x = [1.0_wp, 1.0_wp, 1.0_wp]
 
 call check_deriv()
 
 ! Set tol to the square root of the machine precision. Unless high precision
 ! solutions are required, this is the recommended setting.
-tol = sqrt(epsilon(1._dp))
+tol = sqrt(epsilon(1._wp))
 
-allocate(wa(5*size(x) + size(fvec)))
-call lmder1(fcn, size(fvec), size(x), x, fvec, fjac, size(fjac, 1), tol, &
-    info, ipvt, wa, size(wa))
-print 1000, enorm(size(fvec), fvec), info, x
-1000 format(5x, 'FINAL L2 NORM OF THE RESIDUALS', d15.7 // &
-            5x, 'EXIT PARAMETER', 16x, i10              // &
-            5x, 'FINAL APPROXIMATE SOLUTION'            // &
-            5x, 3d15.7)
+call lmder1(fcn, m, n, x, fvec, fjac, m, tol, info, ipvt, wa, lwa)
+
+write(nwrite, '(5x,a,d15.7//,5x,a,16x,i10//,5x,a//(5x,3d15.7))') &
+        'FINAL L2 NORM OF THE RESIDUALS', enorm(m, fvec), &
+        'EXIT PARAMETER', info, &
+        'FINAL APPROXIMATE SOLUTION', x
 
 contains
 
 subroutine check_deriv()
-integer :: iflag
-real(dp) :: xp(size(x)), fvecp(size(fvec)), err(size(fvec))
-call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, &
-    1, err)
-iflag = 1
-call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag)
-iflag = 2
-call fcn(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), iflag)
-iflag = 1
-call fcn(size(fvec), size(x), xp, fvecp, fjac, size(fjac, 1), iflag)
-call chkder(size(fvec), size(x), x, fvec, fjac, size(fjac, 1), xp, fvecp, &
-    2, err)
-print *, "Derivatives check (1.0 is correct, 0.0 is incorrect):"
-print *, err
-end subroutine
+
+    integer :: iflag
+    real(wp) :: xp(n), fvecp(m), err(m)
+
+    call chkder(m, n, x, fvec, fjac, m, xp, fvecp, 1, err)
+    iflag = 1
+    call fcn(m, n, x, fvec, fjac, m, iflag)
+    iflag = 2
+    call fcn(m, n, x, fvec, fjac, m, iflag)
+    iflag = 1
+    call fcn(m, n, xp, fvecp, fjac, m, iflag)
+    call chkder(m, n, x, fvec, fjac, m, xp, fvecp, 2, err)
+
+    write(nwrite, '(a)') 'Derivatives check (1.0 is correct, 0.0 is incorrect):'
+    write(nwrite,'(1p,(5x,3d15.7))') err
+    if (any(abs(err-1.0_wp)>epsilon(1.0_wp))) error stop 'Derivative check failed'
+
+end subroutine check_deriv
+
+subroutine fcn(m, n, x, fvec, fjac, ldfjac, iflag)
+
+    integer, intent(in) :: m
+    integer, intent(in) :: n
+    real(wp), intent(in) :: x(n)
+    real(wp), intent(inout) :: fvec(m)
+    real(wp), intent(inout) :: fjac(ldfjac, n)
+    integer, intent(in) :: ldfjac
+    integer, intent(inout) :: iflag
+
+    integer :: i
+    real(wp) :: tmp1, tmp2, tmp3, tmp4
+
+    real(wp),parameter :: y(15) = [1.4e-1_wp, 1.8e-1_wp, 2.2e-1_wp, 2.5e-1_wp, 2.9e-1_wp, &
+                                   3.2e-1_wp, 3.5e-1_wp, 3.9e-1_wp, 3.7e-1_wp, 5.8e-1_wp, &
+                                   7.3e-1_wp, 9.6e-1_wp, 1.34e0_wp, 2.1e0_wp,  4.39e0_wp]
+
+    if (iflag == 1) then
+        do i = 1, 15
+            tmp1 = i
+            tmp2 = 16 - i
+            tmp3 = tmp1
+            if (i > 8) tmp3 = tmp2
+            fvec(i) = y(i) - (x(1) + tmp1/(x(2)*tmp2 + x(3)*tmp3))
+        end do
+    else
+        do i = 1, 15
+            tmp1 = i
+            tmp2 = 16 - i
+            tmp3 = tmp1
+            if (i > 8) tmp3 = tmp2
+            tmp4 = (x(2)*tmp2 + x(3)*tmp3)**2
+            fjac(i,1) = -1.0_wp
+            fjac(i,2) = tmp1*tmp2/tmp4
+            fjac(i,3) = tmp1*tmp3/tmp4
+        end do
+    end if
+
+    end subroutine fcn
 
 end program