Merge pull request #63 from jacobwilliams/62-modernize-examples

Modernize the examples

Author: jacobwilliams

examples/example_hybrd.f90

@@ -5,52 +5,41 @@
 !>                                 -x(8) + (3-2*x(9))*x(9) = -1
 program example_hybrd
 
-    use minpack_module, 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/
+    use minpack_module, only: wp, hybrd, enorm, dpmpar
+    use iso_fortran_env, only: nwrite => output_unit
 
-    n = 9
-
-    !> The following starting values provide a rough solution.
-    do j = 1, 9
-        x(j) = -1.0d0
-    end do
+    implicit none
 
-    ldfjac = 9
-    lr = 45
+    integer,parameter :: n = 9
+    integer,parameter :: ldfjac = n
+    integer,parameter :: lr = (n*(n+1))/2
 
-    !> 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))
+    integer :: maxfev, ml, mu, mode, nprint, info, nfev
+    real(wp) :: epsfcn, factor, fnorm, xtol
+    real(wp) :: x(n), fvec(n), diag(n), fjac(n, n), r(lr), qtf(n), &
+                wa1(n), wa2(n), wa3(n), wa4(n)
 
+    xtol = sqrt(dpmpar(1))  ! square root of the machine precision.
     maxfev = 2000
     ml = 1
     mu = 1
-    epsfcn = 0.0d0
+    epsfcn = 0.0_wp
     mode = 2
-    do j = 1, 9
-        diag(j) = 1.0d0
-    end do
-    factor = 1.0d2
+    factor = 100.0_wp
     nprint = 0
+    diag = 1.0_wp
+    x = -1.0_wp  !  starting values to provide a rough solution.
 
     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))
+    write (nwrite, '(5x,a,d15.7//5x,a,i10//5x,a,16x,i10//5x,a//(5x,3d15.7))') &
+           "FINAL L2 NORM OF THE RESIDUALS", fnorm, &
+           "NUMBER OF FUNCTION EVALUATIONS", nfev, &
+           "EXIT PARAMETER", info, &
+           "FINAL APPROXIMATE SOLUTION", x
 
     !> Results obtained with different compilers or machines
     !>  may be slightly different.
@@ -75,28 +64,29 @@ subroutine fcn(n, x, fvec, iflag)
         implicit none
         integer, intent(in) :: n
         integer, intent(inout) :: iflag
-        double precision, intent(in) :: x(n)
-        double precision, intent(out) :: 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
+        real(wp), intent(in) :: x(n)
+        real(wp), intent(out) :: fvec(n)
+
+        integer :: k !! counter
+        real(wp) :: temp, temp1, temp2
+
+        real(wp),parameter :: zero = 0.0_wp
+        real(wp),parameter :: one = 1.0_wp
+        real(wp),parameter :: two = 2.0_wp
+        real(wp),parameter :: three = 3.0_wp
+
+        if (iflag == 0) then
+            !! Insert print statements here when nprint is positive.
+        else
+            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
+        end if
 
     end subroutine fcn
 

examples/example_hybrd1.f90

@@ -6,37 +6,33 @@
 !>                              -x(8) + (3-2*x(9))*x(9) = -1
 program example_hybrd1
 
-    use minpack_module, only: hybrd1, dpmpar, enorm
+    use minpack_module, only: wp, hybrd1, dpmpar, enorm
+    use iso_fortran_env, only: nwrite => output_unit
+
     implicit none
-    integer j, n, info, lwa, nwrite
-    double precision tol, fnorm
-    double precision x(9), fvec(9), wa(180)
 
-    !> Logical output unit is assumed to be number 6.
-    data nwrite/6/
+    integer,parameter :: n = 9
+    integer,parameter :: lwa = (n*(3*n+13))/2
 
-    n = 9
+    integer :: j, info
+    real(wp) :: tol, fnorm
+    real(wp) :: x(n), fvec(n), wa(lwa)
 
     !> The following starting values provide a rough solution.
-    do j = 1, 9
-        x(j) = -1.d0
-    end do
-
-    lwa = 180
+    x = -1.0_wp
 
     !> 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))
+    tol = sqrt(dpmpar(1))
 
     call hybrd1(fcn, n, x, fvec, tol, info, wa, lwa)
     fnorm = enorm(n, fvec)
-    write (nwrite, 1000) fnorm, info, (x(j), j=1, n)
 
-1000 format(5x, "FINAL L2 NORM OF THE RESIDUALS", d15.7// &
-           5x, "EXIT PARAMETER", 16x, i10// &
-           5x, "FINAL APPROXIMATE SOLUTION"// &
-           (5x, 3d15.7))
+    write (nwrite, '(5x,a,d15.7//5x,a,16x,i10//5x,a//(5x,3d15.7))') &
+            "FINAL L2 NORM OF THE RESIDUALS", fnorm, &
+            "EXIT PARAMETER", info, &
+            "FINAL APPROXIMATE SOLUTION", x
 
     !> Results obtained with different compilers or machines
     !>  may be slightly different.
@@ -59,12 +55,16 @@ subroutine fcn(n, x, fvec, iflag)
         implicit none
         integer, intent(in) :: n
         integer, intent(inout) :: iflag
-        double precision, intent(in) :: x(n)
-        double precision, intent(out) :: fvec(n)
+        real(wp), intent(in) :: x(n)
+        real(wp), intent(out) :: fvec(n)
+
+        integer :: k
+        real(wp) :: temp, temp1, temp2
 
-        integer k
-        double precision one, temp, temp1, temp2, three, two, zero
-        data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
+        real(wp),parameter :: zero = 0.0_wp
+        real(wp),parameter :: one = 1.0_wp
+        real(wp),parameter :: two = 2.0_wp
+        real(wp),parameter :: three = 3.0_wp
 
         do k = 1, n
             temp = (three - two*x(k))*x(k)

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