refactor primes example

Author: jacobwilliams

examples/example_primes.f90

@@ -1,22 +1,25 @@
-module types
-implicit none
-private
-public dp
-
-integer, parameter :: dp=kind(0d0)
-
-end module
-
 module find_fit_module
 
-! This module contains a general function find_fit() for a nonlinear least
-! squares fitting. The function can fit any nonlinear expression to any data.
+!! This module contains a general function find_fit() for a nonlinear least
+!! squares fitting. The function can fit any nonlinear expression to any data.
+
+use minpack_module, only: wp, lmdif1
 
-use minpack_module, only: lmdif1
-use types, only: dp
 implicit none
+
 private
-public find_fit
+
+abstract interface
+    function expr_f(x, pars) result(y)
+        import :: wp
+        implicit none
+        real(wp), intent(in) :: x(:)
+        real(wp), intent(in) :: pars(:)
+        real(wp) :: y(size(x))
+    end function expr_f
+end interface
+
+public :: find_fit
 
 contains
 
@@ -26,71 +29,68 @@ subroutine find_fit(data_x, data_y, expr, pars)
 ! parameters 'pars' and it must return the evaluated expression on the
 ! array 'x'. The arrays 'data_x' and 'data_y' must have the same
 ! length.
-real(dp), intent(in) :: data_x(:), data_y(:)
-interface
-    function expr(x, pars) result(y)
-    use types, only: dp
-    implicit none
-    real(dp), intent(in) :: x(:), pars(:)
-    real(dp) :: y(size(x))
-    end function
-end interface
-real(dp), intent(inout) :: pars(:)
+real(wp), intent(in) :: data_x(:)
+real(wp), intent(in) :: data_y(:)
+procedure(expr_f) :: expr
+real(wp), intent(inout) :: pars(:)
 
-real(dp) :: tol, fvec(size(data_x))
+real(wp) :: tol, fvec(size(data_x))
 integer :: iwa(size(pars)), info, m, n
-real(dp), allocatable :: wa(:)
+real(wp), allocatable :: wa(:)
 
-tol = sqrt(epsilon(1._dp))
+tol = sqrt(epsilon(1.0_wp))
 m = size(fvec)
 n = size(pars)
 allocate(wa(m*n + 5*n + m))
 call lmdif1(fcn, m, n, pars, fvec, tol, info, iwa, wa, size(wa))
-if (info /= 1) stop "failed to converge"
+if (info /= 1) error stop "failed to converge"
 
 contains
 
-subroutine fcn(m, n, x, fvec, iflag)
-integer, intent(in) :: m, n
-integer, intent(inout) :: iflag
-real(dp), intent(in) :: x(n)
-real(dp), intent(out) :: fvec(m)
-! Suppress compiler warning:
-fvec(1) = iflag
-fvec = data_y - expr(data_x, x)
-end subroutine
+    subroutine fcn(m, n, x, fvec, iflag)
+    integer, intent(in) :: m
+    integer, intent(in) :: n
+    integer, intent(inout) :: iflag
+    real(wp), intent(in) :: x(n)
+    real(wp), intent(out) :: fvec(m)
+    fvec(1) = iflag ! Suppress compiler warning:
+    fvec = data_y - expr(data_x, x)
+    end subroutine fcn
 
-end subroutine
-
-end module
+end subroutine find_fit
 
+end module find_fit_module
 
 program example_primes
 
-! Find a nonlinear fit of the form a*x*log(b + c*x) to a list of primes.
+!! Find a nonlinear fit of the form `a*x*log(b + c*x)` to a list of primes.
 
 use find_fit_module, only: find_fit
-use types, only: dp
+use minpack_module, only: wp
+use iso_fortran_env, only: nwrite => output_unit
+
 implicit none
 
-real(dp) :: pars(3)
-real(dp), parameter :: y(*) = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, &
-    37, 41, 43, 47, 53, 59, 61, 67, 71]
+real(wp) :: pars(3)
+real(wp), parameter :: y(*) = real([2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, &
+                                    37, 41, 43, 47, 53, 59, 61, 67, 71], wp)
 integer :: i
-pars = [1._dp, 1._dp, 1._dp]
-call find_fit([(real(i, dp), i=1,size(y))], y, expression, pars)
-print *, pars
+pars = [1.0_wp, 1.0_wp, 1.0_wp]
+call find_fit([(real(i, wp), i=1,size(y))], y, expression, pars)
+
+write(nwrite, '(1p,5x,a/(5x,3e15.7))') 'pars: ',pars
 
 contains
 
 function expression(x, pars) result(y)
-real(dp), intent(in) :: x(:), pars(:)
-real(dp) :: y(size(x))
-real(dp) :: a, b, c
+real(wp), intent(in) :: x(:)
+real(wp), intent(in) :: pars(:)
+real(wp) :: y(size(x))
+real(wp) :: a, b, c
 a = pars(1)
 b = pars(2)
 c = pars(3)
 y = a*x*log(b + c*x)
-end function
+end function expression
 
-end program
+end program example_primes