Permalink
Browse files

Add example_primes.f90

  • Loading branch information...
1 parent 56b802f commit d52ca03df7aa63bc94d1d77a0f9f55b2b1f7ea70 @certik committed Feb 25, 2012
Showing with 98 additions and 0 deletions.
  1. +3 −0 examples/CMakeLists.txt
  2. +95 −0 examples/example_primes.f90
@@ -5,3 +5,6 @@ target_link_libraries(example_lmder1 minpack)
add_executable(example_lmdif1 example_lmdif1.f90)
target_link_libraries(example_lmdif1 minpack)
+
+add_executable(example_primes example_primes.f90)
+target_link_libraries(example_primes minpack)
@@ -0,0 +1,95 @@
+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.
+
+use minpack, only: lmdif1
+use types, only: dp
+implicit none
+private
+public find_fit
+
+contains
+
+subroutine find_fit(data_x, data_y, expr, pars)
+! Fits the (data_x, data_y) arrays with the function expr(x, pars).
+! The user can provide any nonlinear function 'expr' depending on any number of
+! 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(dp) :: tol, fvec(size(data_x))
+integer :: iwa(size(pars)), info, m, n
+real(dp), allocatable :: wa(:)
+
+tol = sqrt(epsilon(1._dp))
+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"
+
+contains
+
+subroutine fcn(m, n, x, fvec, iflag)
+integer, intent(in) :: m, n, 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
+
+end subroutine
+
+end module
+
+
+program example_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
+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]
+integer :: i
+pars = [1._dp, 1._dp, 1._dp]
+call find_fit([(real(i, dp), i=1,size(y))], y, expression, pars)
+print *, pars
+
+contains
+
+function expression(x, pars) result(y)
+real(dp), intent(in) :: x(:), pars(:)
+real(dp) :: y(size(x))
+real(dp) :: a, b, c
+a = pars(1)
+b = pars(2)
+c = pars(3)
+y = a*x*log(b + c*x)
+end function
+
+end program

0 comments on commit d52ca03

Please sign in to comment.