modernize hybrd1 example

Author: jacobwilliams

examples/example_hybrd.f90

@@ -19,7 +19,7 @@ program example_hybrd
     real(wp) :: x(n), fvec(n), diag(n), fjac(n, n), r(lr), qtf(n), &
                 wa1(n), wa2(n), wa3(n), wa4(n)
 
-    xtol = dsqrt(dpmpar(1))  ! square root of the machine precision.
+    xtol = sqrt(dpmpar(1))  ! square root of the machine precision.
     maxfev = 2000
     ml = 1
     mu = 1

examples/example_hybrd1.f90

@@ -7,36 +7,32 @@
 program example_hybrd1
 
     use minpack_module, only: wp, hybrd1, dpmpar, enorm
+    use iso_fortran_env, only: nwrite => output_unit
+
     implicit none
-    integer j, n, info, lwa, nwrite
-    real(wp) tol, fnorm
-    real(wp) 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.
@@ -62,9 +58,13 @@ subroutine fcn(n, x, fvec, iflag)
         real(wp), intent(in) :: x(n)
         real(wp), intent(out) :: fvec(n)
 
-        integer k
-        real(wp) one, temp, temp1, temp2, three, two, zero
-        data zero, one, two, three/0.0d0, 1.0d0, 2.0d0, 3.0d0/
+        integer :: k
+        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
 
         do k = 1, n
             temp = (three - two*x(k))*x(k)