From 18a7f677279dc6402f27306cd270b400b5c6f171 Mon Sep 17 00:00:00 2001 From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com> Date: Sat, 23 Apr 2022 22:38:05 +0200 Subject: [PATCH] Allow passing of user data for callback procedure --- examples/example_hybrd.f90 | 3 +- examples/example_hybrd1.f90 | 3 +- examples/example_lmder1.f90 | 3 +- examples/example_lmdif1.f90 | 5 +- examples/example_primes.f90 | 3 +- meson.build | 1 + src/minpack.f90 | 120 ++-- src/minpack_capi.f90 | 314 +++++---- src/minpack_legacy.f90 | 1231 +++++++++++++++++++++++++++++++++++ test/test_hybrd.f90 | 5 +- test/test_hybrj.f90 | 3 +- test/test_lmder.f90 | 5 +- test/test_lmdif.f90 | 5 +- test/test_lmstr.f90 | 5 +- 14 files changed, 1517 insertions(+), 189 deletions(-) create mode 100644 src/minpack_legacy.f90 diff --git a/examples/example_hybrd.f90 b/examples/example_hybrd.f90 index 0db4319..a416143 100644 --- a/examples/example_hybrd.f90 +++ b/examples/example_hybrd.f90 @@ -5,7 +5,8 @@ !> -x(8) + (3-2*x(9))*x(9) = -1 program example_hybrd - use minpack_module, only: wp, hybrd, enorm, dpmpar + use minpack_legacy, only: hybrd + use minpack_module, only: wp, enorm, dpmpar use iso_fortran_env, only: nwrite => output_unit implicit none diff --git a/examples/example_hybrd1.f90 b/examples/example_hybrd1.f90 index bd5beb4..e827f19 100644 --- a/examples/example_hybrd1.f90 +++ b/examples/example_hybrd1.f90 @@ -6,7 +6,8 @@ !> -x(8) + (3-2*x(9))*x(9) = -1 program example_hybrd1 - use minpack_module, only: wp, hybrd1, dpmpar, enorm + use minpack_legacy, only: hybrd1 + use minpack_module, only: wp, dpmpar, enorm use iso_fortran_env, only: nwrite => output_unit implicit none diff --git a/examples/example_lmder1.f90 b/examples/example_lmder1.f90 index cfa1bd4..db61c76 100644 --- a/examples/example_lmder1.f90 +++ b/examples/example_lmder1.f90 @@ -1,6 +1,7 @@ program example_lmder1 -use minpack_module, only: wp, enorm, lmder1, chkder +use minpack_legacy, only: lmder1, chkder +use minpack_module, only: wp, enorm use iso_fortran_env, only: nwrite => output_unit implicit none diff --git a/examples/example_lmdif1.f90 b/examples/example_lmdif1.f90 index a12bc5d..513e04a 100644 --- a/examples/example_lmdif1.f90 +++ b/examples/example_lmdif1.f90 @@ -1,6 +1,7 @@ program example_lmdif1 -use minpack_module, only: wp, enorm, lmdif1 +use minpack_legacy, only: lmdif1 +use minpack_module, only: wp, enorm use iso_fortran_env, only: nwrite => output_unit implicit none @@ -56,4 +57,4 @@ subroutine fcn(m, n, x, fvec, iflag) end subroutine fcn -end program example_lmdif1 \ No newline at end of file +end program example_lmdif1 diff --git a/examples/example_primes.f90 b/examples/example_primes.f90 index a17a76f..6b8b771 100644 --- a/examples/example_primes.f90 +++ b/examples/example_primes.f90 @@ -3,7 +3,8 @@ 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_module, only: wp, lmdif1 +use minpack_legacy, only: lmdif1 +use minpack_module, only: wp implicit none diff --git a/meson.build b/meson.build index 7d5e9e8..ebfdd0a 100644 --- a/meson.build +++ b/meson.build @@ -16,6 +16,7 @@ endif srcs = files( 'src/minpack.f90', + 'src/minpack_legacy.f90', ) if get_option('api') srcs += files('src/minpack_capi.f90') diff --git a/src/minpack.f90 b/src/minpack.f90 index 6e1a240..4bcecbe 100644 --- a/src/minpack.f90 +++ b/src/minpack.f90 @@ -22,7 +22,7 @@ module minpack_module real(wp), parameter, private :: zero = 0.0_wp abstract interface - subroutine func(n, x, fvec, iflag) + subroutine func(n, x, fvec, iflag, user_data) !! user-supplied subroutine for [[hybrd]], [[hybrd1]], and [[fdjac1]] import :: wp implicit none @@ -30,9 +30,10 @@ subroutine func(n, x, fvec, iflag) real(wp), intent(in) :: x(n) !! independent variable vector real(wp), intent(out) :: fvec(n) !! value of function at `x` integer, intent(inout) :: iflag !! set to <0 to terminate execution + class(*), intent(inout), optional :: user_data end subroutine func - subroutine func2(m, n, x, fvec, iflag) + subroutine func2(m, n, x, fvec, iflag, user_data) !! user-supplied subroutine for [[fdjac2]], [[lmdif]], and [[lmdif1]] import :: wp implicit none @@ -43,9 +44,10 @@ subroutine func2(m, n, x, fvec, iflag) integer, intent(inout) :: iflag !! the value of iflag should not be changed unless !! the user wants to terminate execution of lmdif. !! in this case set iflag to a negative integer. + class(*), intent(inout), optional :: user_data end subroutine func2 - subroutine fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag) + subroutine fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag, user_data) !! user-supplied subroutine for [[hybrj]] and [[hybrj1]] import :: wp implicit none @@ -62,9 +64,10 @@ subroutine fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag) !! the value of iflag should not be changed by fcn unless !! the user wants to terminate execution of hybrj. !! in this case set iflag to a negative integer. + class(*), intent(inout), optional :: user_data end subroutine fcn_hybrj - subroutine fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag) + subroutine fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag, user_data) !! user-supplied subroutine for [[lmder]] and [[lmder1]] import :: wp implicit none @@ -82,9 +85,10 @@ subroutine fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag) real(wp), intent(in) :: x(n) !! independent variable vector real(wp), intent(inout) :: fvec(m) !! value of function at `x` real(wp), intent(inout) :: fjac(ldfjac, n) !! jacobian matrix at `x` + class(*), intent(inout), optional :: user_data end subroutine fcn_lmder - subroutine fcn_lmstr(m, n, x, fvec, fjrow, iflag) + subroutine fcn_lmstr(m, n, x, fvec, fjrow, iflag, user_data) import :: wp implicit none integer, intent(in) :: m !! the number of functions. @@ -100,6 +104,7 @@ subroutine fcn_lmstr(m, n, x, fvec, fjrow, iflag) real(wp), intent(in) :: x(n) !! independent variable vector real(wp), intent(inout) :: fvec(m) !! value of function at `x` real(wp), intent(inout) :: fjrow(n) !! jacobian row + class(*), intent(inout), optional :: user_data end subroutine fcn_lmstr end interface @@ -431,7 +436,8 @@ end function enorm ! a banded form, then function evaluations are saved by only ! approximating the nonzero terms. - subroutine fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, Iflag, Ml, Mu, Epsfcn, Wa1, Wa2) + subroutine fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, Iflag, Ml, Mu, Epsfcn, Wa1, Wa2, & + & user_data) implicit none @@ -467,6 +473,7 @@ subroutine fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, Iflag, Ml, Mu, Epsfcn, Wa1, Wa2 real(wp), intent(inout) :: Wa2(n) !! work array of length n. if ml + mu + 1 is at !! least n, then the jacobian is considered dense, and wa2 is !! not referenced. + class(*), intent(inout), optional :: user_data integer :: i, j, k, msum real(wp) :: eps, h, temp @@ -482,7 +489,7 @@ subroutine fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, Iflag, Ml, Mu, Epsfcn, Wa1, Wa2 if (h == zero) h = eps x(j) = Wa2(j) + h end do - call fcn(n, x, Wa1, Iflag) + call fcn(n, x, Wa1, Iflag, user_data) if (Iflag < 0) return do j = k, n, msum x(j) = Wa2(j) @@ -501,7 +508,7 @@ subroutine fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, Iflag, Ml, Mu, Epsfcn, Wa1, Wa2 h = eps*abs(temp) if (h == zero) h = eps x(j) = temp + h - call fcn(n, x, Wa1, Iflag) + call fcn(n, x, Wa1, Iflag, user_data) if (Iflag < 0) return x(j) = temp do i = 1, n @@ -519,7 +526,7 @@ end subroutine fdjac1 ! to the m by n jacobian matrix associated with a specified ! problem of m functions in n variables. - subroutine fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, Iflag, Epsfcn, Wa) + subroutine fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, Iflag, Epsfcn, Wa, user_data) implicit none @@ -546,6 +553,7 @@ subroutine fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, Iflag, Epsfcn, Wa) real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output m by n array which contains the !! approximation to the jacobian matrix evaluated at x. real(wp), intent(inout) :: Wa(m) !! a work array of length m. + class(*), intent(inout), optional :: user_data integer :: i, j real(wp) :: eps, h, temp @@ -556,7 +564,7 @@ subroutine fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, Iflag, Epsfcn, Wa) h = eps*abs(temp) if (h == zero) h = eps x(j) = temp + h - call fcn(m, n, x, Wa, Iflag) + call fcn(m, n, x, Wa, Iflag, user_data) if (Iflag < 0) return x(j) = temp do i = 1, m @@ -577,7 +585,7 @@ end subroutine fdjac2 subroutine 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) + Wa2, Wa3, Wa4, user_data) implicit none @@ -664,6 +672,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & real(wp), intent(inout) :: wa2(n) !! work array real(wp), intent(inout) :: wa3(n) !! work array real(wp), intent(inout) :: wa4(n) !! work array + class(*), intent(inout), optional :: user_data integer :: i, iflag, iter, j, jm1, l, msum, ncfail, ncsuc, nslow1, nslow2 integer :: iwa(1) @@ -693,7 +702,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & ! and calculate its norm. iflag = 1 - call fcn(n, x, Fvec, iflag) + call fcn(n, x, Fvec, iflag, user_data) Nfev = 1 if (iflag < 0) goto 300 fnorm = enorm(n, Fvec) @@ -718,7 +727,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & ! calculate the jacobian matrix. iflag = 2 - call fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, iflag, Ml, Mu, Epsfcn, Wa1, Wa2) + call fdjac1(fcn, n, x, Fvec, Fjac, Ldfjac, iflag, Ml, Mu, Epsfcn, Wa1, Wa2, user_data) Nfev = Nfev + msum if (iflag < 0) goto 300 @@ -798,7 +807,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & 200 if (Nprint > 0) then iflag = 0 - if (mod(iter - 1, Nprint) == 0) call fcn(n, x, Fvec, iflag) + if (mod(iter - 1, Nprint) == 0) call fcn(n, x, Fvec, iflag, user_data) if (iflag < 0) goto 300 end if @@ -822,7 +831,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & ! evaluate the function at x + p and calculate its norm. iflag = 1 - call fcn(n, Wa2, Wa4, iflag) + call fcn(n, Wa2, Wa4, iflag, user_data) Nfev = Nfev + 1 if (iflag >= 0) then fnorm1 = enorm(n, Wa4) @@ -939,7 +948,7 @@ subroutine hybrd(fcn, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & 300 if (iflag < 0) Info = iflag iflag = 0 - if (Nprint > 0) call fcn(n, x, Fvec, iflag) + if (Nprint > 0) call fcn(n, x, Fvec, iflag, user_data) end subroutine hybrd !***************************************************************************************** @@ -954,7 +963,7 @@ end subroutine hybrd ! the jacobian is then calculated by a forward-difference ! approximation. - subroutine hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa) + subroutine hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa, user_data) implicit none @@ -985,6 +994,7 @@ subroutine hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa) integer, intent(in) :: Lwa !! a positive integer input variable not less than !! (n*(3*n+13))/2. real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + class(*), intent(inout), optional :: user_data integer :: index, j, lr, maxfev, ml, mode, mu, nfev, nprint real(wp) :: epsfcn, xtol @@ -1011,7 +1021,8 @@ subroutine hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa) index = 6*n + lr call hybrd(fcn, n, x, Fvec, xtol, maxfev, ml, mu, epsfcn, Wa(1), mode, & factor, nprint, Info, nfev, Wa(index + 1), n, Wa(6*n + 1), lr, & - Wa(n + 1), Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1)) + Wa(n + 1), Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1), & + user_data) if (Info == 5) Info = 4 end if @@ -1027,7 +1038,7 @@ end subroutine hybrd1 subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & Factor, Nprint, Info, Nfev, Njev, r, Lr, Qtf, Wa1, Wa2, & - Wa3, Wa4) + Wa3, Wa4, user_data) implicit none @@ -1105,6 +1116,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & real(wp), intent(inout) :: Wa2(n) !! work array of length n. real(wp), intent(inout) :: Wa3(n) !! work array of length n. real(wp), intent(inout) :: Wa4(n) !! work array of length n. + class(*), intent(inout), optional :: user_data integer :: i, iflag, iter, j, jm1, l, ncfail, ncsuc, nslow1, nslow2 integer :: iwa(1) @@ -1135,7 +1147,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & ! and calculate its norm. iflag = 1 - call fcn(n, x, Fvec, Fjac, Ldfjac, iflag) + call fcn(n, x, Fvec, Fjac, Ldfjac, iflag, user_data) Nfev = 1 if (iflag < 0) goto 300 fnorm = enorm(n, Fvec) @@ -1155,7 +1167,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & ! calculate the jacobian matrix. iflag = 2 - call fcn(n, x, Fvec, Fjac, Ldfjac, iflag) + call fcn(n, x, Fvec, Fjac, Ldfjac, iflag, user_data) Njev = Njev + 1 if (iflag < 0) goto 300 @@ -1238,7 +1250,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & 200 if (Nprint > 0) then iflag = 0 if (mod(iter - 1, Nprint) == 0) & - call fcn(n, x, Fvec, Fjac, Ldfjac, iflag) + call fcn(n, x, Fvec, Fjac, Ldfjac, iflag, user_data) if (iflag < 0) goto 300 end if @@ -1262,7 +1274,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & ! evaluate the function at x + p and calculate its norm. iflag = 1 - call fcn(n, Wa2, Wa4, Fjac, Ldfjac, iflag) + call fcn(n, Wa2, Wa4, Fjac, Ldfjac, iflag, user_data) Nfev = Nfev + 1 if (iflag >= 0) then fnorm1 = enorm(n, Wa4) @@ -1380,7 +1392,7 @@ subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & 300 if (iflag < 0) Info = iflag iflag = 0 - if (Nprint > 0) call fcn(n, x, Fvec, Fjac, Ldfjac, iflag) + if (Nprint > 0) call fcn(n, x, Fvec, Fjac, Ldfjac, iflag, user_data) end subroutine hybrj !***************************************************************************************** @@ -1394,7 +1406,7 @@ end subroutine hybrj ! must provide a subroutine which calculates the functions ! and the jacobian. - subroutine hybrj1(fcn, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa) + subroutine hybrj1(fcn, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa, user_data) implicit none @@ -1431,6 +1443,7 @@ subroutine hybrj1(fcn, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa) !! orthogonal matrix q produced by the qr factorization !! of the final approximate jacobian. real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + class(*), intent(inout), optional :: user_data integer :: j, lr, maxfev, mode, nfev, njev, nprint real(wp) :: xtol @@ -1453,7 +1466,7 @@ subroutine hybrj1(fcn, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa) lr = (n*(n + 1))/2 call hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, xtol, maxfev, Wa(1), mode, & factor, nprint, Info, nfev, njev, Wa(6*n + 1), lr, Wa(n + 1), & - Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1)) + Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1), user_data) if (Info == 5) Info = 4 end if @@ -1469,7 +1482,7 @@ end subroutine hybrj1 subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & Diag, Mode, Factor, Nprint, Info, Nfev, Njev, Ipvt, Qtf, & - Wa1, Wa2, Wa3, Wa4) + Wa1, Wa2, Wa3, Wa4, user_data) implicit none @@ -1575,6 +1588,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & real(wp), intent(inout) :: Wa2(n) !! work array of length n. real(wp), intent(inout) :: Wa3(n) !! work array of length n. real(wp), intent(inout) :: Wa4(m) !! work array of length m. + class(*), intent(inout), optional :: user_data integer :: i, iflag, iter, j, l real(wp) :: actred, delta, dirder, fnorm, fnorm1, gnorm, par, & @@ -1606,7 +1620,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & ! and calculate its norm. iflag = 1 - call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag) + call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag, user_data) Nfev = 1 if (iflag >= 0) then fnorm = enorm(m, Fvec) @@ -1621,7 +1635,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & ! calculate the jacobian matrix. 20 iflag = 2 - call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag) + call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag, user_data) Njev = Njev + 1 if (iflag >= 0) then @@ -1630,7 +1644,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & if (Nprint > 0) then iflag = 0 if (mod(iter - 1, Nprint) == 0) & - call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag) + call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag, user_data) if (iflag < 0) goto 100 end if @@ -1732,7 +1746,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & ! evaluate the function at x + p and calculate its norm. iflag = 1 - call fcn(m, n, Wa2, Wa4, Fjac, Ldfjac, iflag) + call fcn(m, n, Wa2, Wa4, Fjac, Ldfjac, iflag, user_data) Nfev = Nfev + 1 if (iflag >= 0) then fnorm1 = enorm(m, Wa4) @@ -1824,7 +1838,7 @@ subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & 100 if (iflag < 0) Info = iflag iflag = 0 - if (Nprint > 0) call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag) + if (Nprint > 0) call fcn(m, n, x, Fvec, Fjac, Ldfjac, iflag, user_data) end subroutine lmder !***************************************************************************************** @@ -1837,7 +1851,7 @@ end subroutine lmder ! general least-squares solver lmder. the user must provide a ! subroutine which calculates the functions and the jacobian. - subroutine lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) + subroutine lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, user_data) implicit none procedure(fcn_lmder) :: fcn !! user-supplied subroutine which @@ -1897,6 +1911,7 @@ subroutine lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) !! part of fjac contains information generated during !! the computation of r. real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + class(*), intent(inout), optional :: user_data integer :: maxfev, mode, nfev, njev, nprint real(wp) :: ftol, gtol, xtol @@ -1918,7 +1933,7 @@ subroutine lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) nprint = 0 call lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, ftol, xtol, gtol, maxfev, & & Wa(1), mode, factor, nprint, Info, nfev, njev, Ipvt, Wa(n + 1)& - & , Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1)) + & , Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1), user_data) if (Info == 8) Info = 4 end if @@ -1935,7 +1950,7 @@ end subroutine lmder1 subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & Mode, Factor, Nprint, Info, Nfev, Fjac, Ldfjac, Ipvt, & - Qtf, Wa1, Wa2, Wa3, Wa4) + Qtf, Wa1, Wa2, Wa3, Wa4, user_data) implicit none procedure(func2) :: fcn !! the user-supplied subroutine which @@ -2044,6 +2059,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & real(wp), intent(inout) :: Wa2(n) !! work array of length n. real(wp), intent(inout) :: Wa3(n) !! work array of length n. real(wp), intent(inout) :: Wa4(m) !! work array of length n. + class(*), intent(inout), optional :: user_data integer :: i, iflag, iter, j, l real(wp) :: actred, delta, dirder, fnorm, & @@ -2075,7 +2091,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & ! and calculate its norm. iflag = 1 - call fcn(m, n, x, Fvec, iflag) + call fcn(m, n, x, Fvec, iflag, user_data) Nfev = 1 if (iflag >= 0) then fnorm = enorm(m, Fvec) @@ -2090,7 +2106,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & ! calculate the jacobian matrix. 20 iflag = 2 - call fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, iflag, Epsfcn, Wa4) + call fdjac2(fcn, m, n, x, Fvec, Fjac, Ldfjac, iflag, Epsfcn, Wa4, user_data) Nfev = Nfev + n if (iflag >= 0) then @@ -2099,7 +2115,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & if (Nprint > 0) then iflag = 0 if (mod(iter - 1, Nprint) == 0) & - call fcn(m, n, x, Fvec, iflag) + call fcn(m, n, x, Fvec, iflag, user_data) if (iflag < 0) goto 100 end if @@ -2202,7 +2218,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & ! evaluate the function at x + p and calculate its norm. iflag = 1 - call fcn(m, n, Wa2, Wa4, iflag) + call fcn(m, n, Wa2, Wa4, iflag, user_data) Nfev = Nfev + 1 if (iflag >= 0) then fnorm1 = enorm(m, Wa4) @@ -2300,7 +2316,7 @@ subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & 100 if (iflag < 0) Info = iflag iflag = 0 - if (Nprint > 0) call fcn(m, n, x, Fvec, iflag) + if (Nprint > 0) call fcn(m, n, x, Fvec, iflag, user_data) end subroutine lmdif !***************************************************************************************** @@ -2314,7 +2330,7 @@ end subroutine lmdif ! subroutine which calculates the functions. the jacobian is ! then calculated by a forward-difference approximation. - subroutine lmdif1(fcn, m, n, x, Fvec, Tol, Info, Iwa, Wa, Lwa) + subroutine lmdif1(fcn, m, n, x, Fvec, Tol, Info, Iwa, Wa, Lwa, user_data) implicit none procedure(func2) :: fcn !! the user-supplied subroutine which @@ -2356,6 +2372,7 @@ subroutine lmdif1(fcn, m, n, x, Fvec, Tol, Info, Iwa, Wa, Lwa) real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains !! the functions evaluated at the output x. real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + class(*), intent(inout), optional :: user_data integer :: maxfev, mode, mp5n, nfev, nprint real(wp) :: epsfcn, ftol, gtol, xtol @@ -2380,7 +2397,8 @@ subroutine lmdif1(fcn, m, n, x, Fvec, Tol, Info, Iwa, Wa, Lwa) mp5n = m + 5*n call lmdif(fcn, m, n, x, Fvec, ftol, xtol, gtol, maxfev, epsfcn, Wa(1), & mode, factor, nprint, Info, nfev, Wa(mp5n + 1), m, Iwa, & - Wa(n + 1), Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1)) + Wa(n + 1), Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1), & + user_data) if (Info == 8) Info = 4 end if @@ -2625,7 +2643,7 @@ end subroutine lmpar subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & Diag, Mode, Factor, Nprint, Info, Nfev, Njev, Ipvt, Qtf, & - Wa1, Wa2, Wa3, Wa4) + Wa1, Wa2, Wa3, Wa4, user_data) implicit none procedure(fcn_lmstr) :: fcn !! user-supplied subroutine which @@ -2727,6 +2745,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & real(wp), intent(inout) :: Wa2(n) !! work array of length n. real(wp), intent(inout) :: Wa3(n) !! work array of length n. real(wp), intent(inout) :: Wa4(m) !! work array of length m. + class(*), intent(inout), optional :: user_data integer :: i, iflag, iter, j, l real(wp) :: actred, delta, dirder, fnorm, & @@ -2760,7 +2779,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & ! and calculate its norm. iflag = 1 - call fcn(m, n, x, Fvec, Wa3, iflag) + call fcn(m, n, x, Fvec, Wa3, iflag, user_data) Nfev = 1 if (iflag < 0) goto 200 fnorm = enorm(m, Fvec) @@ -2776,7 +2795,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & 100 if (Nprint > 0) then iflag = 0 - if (mod(iter - 1, Nprint) == 0) call fcn(m, n, x, Fvec, Wa3, iflag) + if (mod(iter - 1, Nprint) == 0) call fcn(m, n, x, Fvec, Wa3, iflag, user_data) if (iflag < 0) goto 200 end if @@ -2793,7 +2812,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & end do iflag = 2 do i = 1, m - call fcn(m, n, x, Fvec, Wa3, iflag) + call fcn(m, n, x, Fvec, Wa3, iflag, user_data) if (iflag < 0) goto 200 temp = Fvec(i) call rwupdt(n, Fjac, Ldfjac, Wa3, Qtf, temp, Wa1, Wa2) @@ -2900,7 +2919,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & ! evaluate the function at x + p and calculate its norm. iflag = 1 - call fcn(m, n, Wa2, Wa4, Wa3, iflag) + call fcn(m, n, Wa2, Wa4, Wa3, iflag, user_data) Nfev = Nfev + 1 if (iflag >= 0) then fnorm1 = enorm(m, Wa4) @@ -2996,7 +3015,7 @@ subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & 200 if (iflag < 0) Info = iflag iflag = 0 - if (Nprint > 0) call fcn(m, n, x, Fvec, Wa3, iflag) + if (Nprint > 0) call fcn(m, n, x, Fvec, Wa3, iflag, user_data) end subroutine lmstr !***************************************************************************************** @@ -3010,7 +3029,7 @@ end subroutine lmstr ! lmstr. the user must provide a subroutine which calculates ! the functions and the rows of the jacobian. - subroutine lmstr1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) + subroutine lmstr1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, user_data) implicit none procedure(fcn_lmstr) :: fcn !! user-supplied subroutine which @@ -3068,6 +3087,7 @@ subroutine lmstr1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) !! part of fjac contains information generated during !! the computation of r. real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + class(*), intent(inout), optional :: user_data integer :: maxfev, mode, nfev, njev, nprint real(wp) :: ftol, gtol, xtol @@ -3091,7 +3111,7 @@ subroutine lmstr1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) nprint = 0 call lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, ftol, xtol, gtol, maxfev, & Wa(1), mode, factor, nprint, Info, nfev, njev, Ipvt, Wa(n + 1), & - Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1)) + Wa(2*n + 1), Wa(3*n + 1), Wa(4*n + 1), Wa(5*n + 1), user_data) if (Info == 8) Info = 4 end if diff --git a/src/minpack_capi.f90 b/src/minpack_capi.f90 index 4386c59..58641ba 100644 --- a/src/minpack_capi.f90 +++ b/src/minpack_capi.f90 @@ -1,10 +1,10 @@ module minpack_capi - use, intrinsic :: iso_c_binding, only : c_int, c_double, c_funptr, c_ptr, c_f_procpointer + use, intrinsic :: iso_c_binding, only : c_int, c_double, c_null_ptr, c_ptr use minpack_module implicit none private - public :: minpack_hybrd,minpack_hybrd1, minpack_hybrj, minpack_hybrj1, & + public :: minpack_hybrd, minpack_hybrd1, minpack_hybrj, minpack_hybrj1, & & minpack_lmdif, minpack_lmdif1, minpack_lmder, minpack_lmder1, & & minpack_chkder @@ -73,6 +73,31 @@ subroutine minpack_fcn_lmstr(m, n, x, fvec, fjrow, iflag, udata) bind(c) end subroutine minpack_fcn_lmstr end interface + type :: hybrd_data + procedure(minpack_func), pointer, nopass :: fcn => null() + type(c_ptr) :: udata = c_null_ptr + end type hybrd_data + + type :: hybrj_data + procedure(minpack_fcn_hybrj), pointer, nopass :: fcn => null() + type(c_ptr) :: udata = c_null_ptr + end type hybrj_data + + type :: lmdif_data + procedure(minpack_func2), pointer, nopass :: fcn => null() + type(c_ptr) :: udata = c_null_ptr + end type lmdif_data + + type :: lmder_data + procedure(minpack_fcn_lmder), pointer, nopass :: fcn => null() + type(c_ptr) :: udata = c_null_ptr + end type lmder_data + + type :: lmstr_data + procedure(minpack_fcn_lmstr), pointer, nopass :: fcn => null() + type(c_ptr) :: udata = c_null_ptr + end type lmstr_data + contains function minpack_dpmpar(i) result(par) bind(c) @@ -113,18 +138,13 @@ subroutine minpack_hybrd(fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mo real(c_double), intent(inout) :: wa4(n) type(c_ptr), value :: udata - call hybrd(wrap_fcn, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, & - & factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4) + type(hybrd_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - contains - subroutine wrap_fcn(n, x, fvec, iflag) - integer, intent(in) :: n - real(wp), intent(in) :: x(n) - real(wp), intent(out) :: fvec(n) - integer, intent(inout) :: iflag - - call fcn(n, x, fvec, iflag, udata) - end subroutine wrap_fcn + call hybrd(wrap_func, n, x, fvec, xtol, maxfev, ml, mu, epsfcn, diag, mode, & + & factor, nprint, info, nfev, fjac, ldfjac, r, lr, qtf, wa1, wa2, wa3, wa4, & + & wrapper) end subroutine minpack_hybrd subroutine minpack_hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa, udata) & @@ -139,19 +159,33 @@ subroutine minpack_hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa, udata) & real(c_double), intent(inout) :: Wa(Lwa) type(c_ptr), value :: udata - call hybrd1(wrap_fcn, n, x, fvec, tol, info, Wa, Lwa) - - contains - subroutine wrap_fcn(n, x, fvec, iflag) - integer, intent(in) :: n - real(wp), intent(in) :: x(n) - real(wp), intent(out) :: fvec(n) - integer, intent(inout) :: iflag + type(hybrd_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - call fcn(n, x, fvec, iflag, udata) - end subroutine wrap_fcn + call hybrd1(wrap_func, n, x, fvec, tol, info, Wa, Lwa, wrapper) end subroutine minpack_hybrd1 + subroutine wrap_func(n, x, fvec, iflag, wrapper) + integer, intent(in) :: n + real(wp), intent(in) :: x(n) + real(wp), intent(out) :: fvec(n) + integer, intent(inout) :: iflag + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(hybrd_data) + call wrapper%fcn(n, x, fvec, iflag, wrapper%udata) + class default + iflag = -1 + end select + end subroutine wrap_func + subroutine minpack_hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, & & factor, nprint, info, nfev, njev, r, lr, qtf, wa1, wa2, wa3, wa4, udata) & & bind(c) @@ -179,20 +213,12 @@ subroutine minpack_hybrj(fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode real(c_double), intent(inout) :: wa4(n) type(c_ptr), value :: udata - call hybrj(wrap_fcn, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, & - & factor, nprint, info, nfev, njev, r, lr, qtf, wa1, wa2, wa3, wa4) - - contains - subroutine wrap_fcn(n, x, fvec, fjac, ldfjac, iflag) - integer, intent(in) :: n - real(wp), intent(in) :: x(n) - integer, intent(in) :: ldfjac - real(wp), intent(inout) :: fvec(n) - real(wp), intent(inout) :: fjac(ldfjac, n) - integer, intent(inout) :: iflag + type(hybrj_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - call fcn(n, x, fvec, fjac, ldfjac, iflag, udata) - end subroutine wrap_fcn + call hybrj(wrap_fcn_hybrj, n, x, fvec, fjac, ldfjac, xtol, maxfev, diag, mode, & + & factor, nprint, info, nfev, njev, r, lr, qtf, wa1, wa2, wa3, wa4, wrapper) end subroutine minpack_hybrj subroutine minpack_hybrj1(fcn, n, x, fvec, fjac, ldfjac, tol, info, wa, lwa, udata) & @@ -209,21 +235,35 @@ subroutine minpack_hybrj1(fcn, n, x, fvec, fjac, ldfjac, tol, info, wa, lwa, uda real(c_double), intent(inout) :: wa(lwa) type(c_ptr), value :: udata - call hybrj1(wrap_fcn, n, x, fvec, fjac, ldfjac, tol, info, wa, lwa) + type(hybrj_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - contains - subroutine wrap_fcn(n, x, fvec, fjac, ldfjac, iflag) - integer, intent(in) :: n - real(wp), intent(in) :: x(n) - integer, intent(in) :: ldfjac - real(wp), intent(inout) :: fvec(n) - real(wp), intent(inout) :: fjac(ldfjac, n) - integer, intent(inout) :: iflag - - call fcn(n, x, fvec, fjac, ldfjac, iflag, udata) - end subroutine wrap_fcn + call hybrj1(wrap_fcn_hybrj, n, x, fvec, fjac, ldfjac, tol, info, wa, lwa, wrapper) end subroutine minpack_hybrj1 + subroutine wrap_fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag, wrapper) + integer, intent(in) :: n + real(wp), intent(in) :: x(n) + integer, intent(in) :: ldfjac + real(wp), intent(inout) :: fvec(n) + real(wp), intent(inout) :: fjac(ldfjac, n) + integer, intent(inout) :: iflag + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(hybrj_data) + call wrapper%fcn(n, x, fvec, fjac, ldfjac, iflag, wrapper%udata) + class default + iflag = -1 + end select + end subroutine wrap_fcn_hybrj + subroutine minpack_lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, & & mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4, & & udata) & @@ -254,19 +294,13 @@ subroutine minpack_lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, d real(c_double), intent(inout) :: wa4(m) type(c_ptr), value :: udata - call lmdif(wrap_fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, & - & mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4) + type(lmdif_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - contains - subroutine wrap_fcn(m, n, x, fvec, iflag) - integer, intent(in) :: m - integer, intent(in) :: n - real(wp), intent(in) :: x(n) - real(wp), intent(out) :: fvec(m) - integer, intent(inout) :: iflag - - call fcn(m, n, x, fvec, iflag, udata) - end subroutine wrap_fcn + call lmdif(wrap_func2, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag, & + & mode, factor, nprint, info, nfev, fjac, ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4, & + & wrapper) end subroutine minpack_lmdif subroutine minpack_lmdif1(fcn, m, n, x, fvec, tol, info, iwa, wa, lwa, udata) & @@ -283,20 +317,34 @@ subroutine minpack_lmdif1(fcn, m, n, x, fvec, tol, info, iwa, wa, lwa, udata) & real(c_double), intent(inout) :: wa(lwa) type(c_ptr), value :: udata - call lmdif1(wrap_fcn, m, n, x, fvec, tol, info, iwa, wa, lwa) + type(lmdif_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - contains - subroutine wrap_fcn(m, n, x, fvec, iflag) - integer, intent(in) :: m - integer, intent(in) :: n - real(wp), intent(in) :: x(n) - real(wp), intent(out) :: fvec(m) - integer, intent(inout) :: iflag - - call fcn(m, n, x, fvec, iflag, udata) - end subroutine wrap_fcn + call lmdif1(wrap_func2, m, n, x, fvec, tol, info, iwa, wa, lwa, wrapper) end subroutine minpack_lmdif1 + subroutine wrap_func2(m, n, x, fvec, iflag, wrapper) + integer, intent(in) :: m + integer, intent(in) :: n + real(wp), intent(in) :: x(n) + real(wp), intent(out) :: fvec(m) + integer, intent(inout) :: iflag + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(lmdif_data) + call wrapper%fcn(m, n, x, fvec, iflag, wrapper%udata) + class default + iflag = -1 + end select + end subroutine wrap_func2 + subroutine minpack_lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4, & & udata) & @@ -327,21 +375,13 @@ subroutine minpack_lmder(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, max real(c_double), intent(inout) :: wa4(m) type(c_ptr), value :: udata - call lmder(wrap_fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & - & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4) - - contains - subroutine wrap_fcn(m, n, x, fvec, fjac, ldfjac, iflag) - integer, intent(in) :: m - integer, intent(in) :: n - integer, intent(in) :: ldfjac - integer, intent(inout) :: iflag - real(wp), intent(in) :: x(n) - real(wp), intent(inout) :: fvec(m) - real(wp), intent(inout) :: fjac(ldfjac, n) - - call fcn(m, n, x, fvec, fjac, ldfjac, iflag, udata) - end subroutine wrap_fcn + type(lmder_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata + + call lmder(wrap_fcn_lmder, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & + & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4, & + & wrapper) end subroutine minpack_lmder subroutine minpack_lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, & @@ -361,22 +401,37 @@ subroutine minpack_lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, real(c_double), intent(inout) :: wa(lwa) type(c_ptr), value :: udata - call lmder1(wrap_fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) - - contains - subroutine wrap_fcn(m, n, x, fvec, fjac, ldfjac, iflag) - integer, intent(in) :: m - integer, intent(in) :: n - integer, intent(in) :: ldfjac - integer, intent(inout) :: iflag - real(wp), intent(in) :: x(n) - real(wp), intent(inout) :: fvec(m) - real(wp), intent(inout) :: fjac(ldfjac, n) + type(lmder_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata - call fcn(m, n, x, fvec, fjac, ldfjac, iflag, udata) - end subroutine wrap_fcn + call lmder1(wrap_fcn_lmder, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, & + & wrapper) end subroutine minpack_lmder1 + subroutine wrap_fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag, wrapper) + integer, intent(in) :: m + integer, intent(in) :: n + integer, intent(in) :: ldfjac + integer, intent(inout) :: iflag + real(wp), intent(in) :: x(n) + real(wp), intent(inout) :: fvec(m) + real(wp), intent(inout) :: fjac(ldfjac, n) + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(lmder_data) + call wrapper%fcn(m, n, x, fvec, fjac, ldfjac, iflag, wrapper%udata) + class default + iflag = -1 + end select + end subroutine wrap_fcn_lmder + subroutine minpack_lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4, & & udata) & @@ -407,19 +462,13 @@ subroutine minpack_lmstr(fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, max real(c_double), intent(inout) :: wa4(m) type(c_ptr), value :: udata - call lmstr(wrap_fcn, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & - & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4) - contains - subroutine wrap_fcn(m, n, x, fvec, fjrow, iflag) - integer, intent(in) :: m - integer, intent(in) :: n - integer, intent(inout) :: iflag - real(wp), intent(in) :: x(n) - real(wp), intent(inout) :: fvec(m) - real(wp), intent(inout) :: fjrow(n) - - call fcn(m, n, x, fvec, fjrow, iflag, udata) - end subroutine wrap_fcn + type(lmstr_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata + + call lmstr(wrap_fcn_lmstr, m, n, x, fvec, fjac, ldfjac, ftol, xtol, gtol, maxfev, & + & diag, mode, factor, nprint, info, nfev, njev, ipvt, qtf, wa1, wa2, wa3, wa4, & + & wrapper) end subroutine minpack_lmstr subroutine minpack_lmstr1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info, ipvt, wa, lwa, & @@ -439,20 +488,37 @@ subroutine minpack_lmstr1(fcn, m, n, x, fvec, fjac, ldfjac, tol, info, ipvt, wa, real(c_double), intent(inout) :: wa(lwa) type(c_ptr), value :: udata - call lmstr1(wrap_fcn, m, n, x, fvec, fjac, ldfjac, tol, info, ipvt, wa, lwa) - contains - subroutine wrap_fcn(m, n, x, fvec, fjrow, iflag) - integer, intent(in) :: m - integer, intent(in) :: n - integer, intent(inout) :: iflag - real(wp), intent(in) :: x(n) - real(wp), intent(inout) :: fvec(m) - real(wp), intent(inout) :: fjrow(n) - - call fcn(m, n, x, fvec, fjrow, iflag, udata) - end subroutine wrap_fcn + type(lmstr_data) :: wrapper + wrapper%fcn => fcn + wrapper%udata = udata + + call lmstr1(wrap_fcn_lmstr, m, n, x, fvec, fjac, ldfjac, tol, info, ipvt, wa, lwa, & + & wrapper) end subroutine minpack_lmstr1 + subroutine wrap_fcn_lmstr(m, n, x, fvec, fjrow, iflag, wrapper) + integer, intent(in) :: m + integer, intent(in) :: n + integer, intent(inout) :: iflag + real(wp), intent(in) :: x(n) + real(wp), intent(inout) :: fvec(m) + real(wp), intent(inout) :: fjrow(n) + + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(lmstr_data) + call wrapper%fcn(m, n, x, fvec, fjrow, iflag, wrapper%udata) + class default + iflag = -1 + end select + end subroutine wrap_fcn_lmstr + subroutine minpack_chkder(m, n, x, fvec, fjac, ldfjac, xp, fvecp, mode, err) & & bind(c) integer(c_int), value :: m diff --git a/src/minpack_legacy.f90 b/src/minpack_legacy.f90 new file mode 100644 index 0000000..5cbde8a --- /dev/null +++ b/src/minpack_legacy.f90 @@ -0,0 +1,1231 @@ +!***************************************************************************************** +!> +! Modernized Minpack +! +!### Authors +! * argonne national laboratory. minpack project. march 1980. +! burton s. garbow, kenneth e. hillstrom, jorge j. more. +! * Jacob Williams, Sept 2021, updated to modern standards. + +module minpack_legacy + use minpack_module, only : & + & hybrd_ => hybrd, hybrd1_ => hybrd1, func_ => func, & + & hybrj_ => hybrj, hybrj1_ => hybrj1, fcn_hybrj_ => fcn_hybrj, & + & lmdif_ => lmdif, lmdif1_ => lmdif1, func2_ => func2, & + & lmder_ => lmder, lmder1_ => lmder1, fcn_lmder_ => fcn_lmder, & + & lmstr_ => lmstr, lmstr1_ => lmstr1, fcn_lmstr_ => fcn_lmstr, & + & chkder_ => chkder, wp + implicit none + private + + public :: & + & hybrd, hybrd1, func, & + & hybrj, hybrj1, fcn_hybrj, & + & lmdif, lmdif1, func2, & + & lmder, lmder1, fcn_lmder, & + & lmstr, lmstr1, fcn_lmstr, & + & chkder + + abstract interface + subroutine func(n, x, fvec, iflag) + !! user-supplied subroutine for [[hybrd]], [[hybrd1]], and [[fdjac1]] + import :: wp + implicit none + integer, intent(in) :: n !! the number of variables. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(out) :: fvec(n) !! value of function at `x` + integer, intent(inout) :: iflag !! set to <0 to terminate execution + end subroutine func + + subroutine func2(m, n, x, fvec, iflag) + !! user-supplied subroutine for [[fdjac2]], [[lmdif]], and [[lmdif1]] + import :: wp + implicit none + integer, intent(in) :: m !! the number of functions. + integer, intent(in) :: n !! the number of variables. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(out) :: fvec(m) !! value of function at `x` + integer, intent(inout) :: iflag !! the value of iflag should not be changed unless + !! the user wants to terminate execution of lmdif. + !! in this case set iflag to a negative integer. + end subroutine func2 + + subroutine fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag) + !! user-supplied subroutine for [[hybrj]] and [[hybrj1]] + import :: wp + implicit none + integer, intent(in) :: n !! the number of variables. + real(wp), dimension(n), intent(in) :: x !! independent variable vector + integer, intent(in) :: ldfjac !! leading dimension of the array fjac. + real(wp), dimension(n), intent(inout) :: fvec !! value of function at `x` + real(wp), dimension(ldfjac, n), intent(inout) :: fjac !! jacobian matrix at `x` + integer, intent(inout) :: iflag !! if iflag = 1 calculate the functions at x and + !! return this vector in fvec. do not alter fjac. + !! if iflag = 2 calculate the jacobian at x and + !! return this matrix in fjac. do not alter fvec. + !! + !! the value of iflag should not be changed by fcn unless + !! the user wants to terminate execution of hybrj. + !! in this case set iflag to a negative integer. + end subroutine fcn_hybrj + + subroutine fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag) + !! user-supplied subroutine for [[lmder]] and [[lmder1]] + import :: wp + implicit none + integer, intent(in) :: m !! the number of functions. + integer, intent(in) :: n !! the number of variables. + integer, intent(in) :: ldfjac !! leading dimension of the array fjac. + integer, intent(inout) :: iflag !! if iflag = 1 calculate the functions at x and + !! return this vector in fvec. do not alter fjac. + !! if iflag = 2 calculate the jacobian at x and + !! return this matrix in fjac. do not alter fvec. + !! + !! the value of iflag should not be changed by fcn unless + !! the user wants to terminate execution of lmder. + !! in this case set iflag to a negative integer. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(inout) :: fvec(m) !! value of function at `x` + real(wp), intent(inout) :: fjac(ldfjac, n) !! jacobian matrix at `x` + end subroutine fcn_lmder + + subroutine fcn_lmstr(m, n, x, fvec, fjrow, iflag) + import :: wp + implicit none + integer, intent(in) :: m !! the number of functions. + integer, intent(in) :: n !! the number of variables. + integer, intent(inout) :: iflag !! if iflag = 1 calculate the functions at x and + !! return this vector in fvec. + !! if iflag = i calculate the (i-1)-st row of the + !! jacobian at x and return this vector in fjrow. + !! + !! the value of iflag should not be changed by fcn unless + !! the user wants to terminate execution of lmstr. + !! in this case set iflag to a negative integer. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(inout) :: fvec(m) !! value of function at `x` + real(wp), intent(inout) :: fjrow(n) !! jacobian row + end subroutine fcn_lmstr + + end interface + + + type :: hybrd_data + procedure(func), pointer, nopass :: fcn => null() + end type hybrd_data + + type :: hybrj_data + procedure(fcn_hybrj), pointer, nopass :: fcn => null() + end type hybrj_data + + type :: lmdif_data + procedure(func2), pointer, nopass :: fcn => null() + end type lmdif_data + + type :: lmder_data + procedure(fcn_lmder), pointer, nopass :: fcn => null() + end type lmder_data + + type :: lmstr_data + procedure(fcn_lmstr), pointer, nopass :: fcn => null() + end type lmstr_data + +contains +!***************************************************************************************** + +!***************************************************************************************** +!> +! this subroutine checks the gradients of m nonlinear functions +! in n variables, evaluated at a point x, for consistency with +! the functions themselves. +! +! the subroutine does not perform reliably if cancellation or +! rounding errors cause a severe loss of significance in the +! evaluation of a function. therefore, none of the components +! of x should be unusually small (in particular, zero) or any +! other value which may cause loss of significance. + + subroutine chkder(m, n, x, Fvec, Fjac, Ldfjac, Xp, Fvecp, Mode, Err) + + implicit none + + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. + integer, intent(in) :: Ldfjac !! a positive integer input parameter not less than m + !! which specifies the leading dimension of the array fjac. + integer, intent(in) :: Mode !! an integer input variable set to 1 on the first call + !! and 2 on the second. other values of mode are equivalent + !! to mode = 1. + !! + !! the user must call chkder twice, + !! first with mode = 1 and then with mode = 2. + !! + !! * mode = 1. **on input**, x must contain the point of evaluation. + !! **on output**, xp is set to a neighboring point. + !! + !! * mode = 2. **on input**, fvec must contain the functions and the + !! rows of fjac must contain the gradients + !! of the respective functions each evaluated + !! at x, and fvecp must contain the functions + !! evaluated at xp. + !! **on output**, err contains measures of correctness of + !! the respective gradients. + real(wp), intent(in) :: x(n) !! input array + real(wp), intent(in) :: Fvec(m) !! an array of length m. on input when mode = 2, + !! fvec must contain the functions evaluated at x. + real(wp), intent(in) :: Fjac(Ldfjac, n) !! an m by n array. on input when mode = 2, + !! the rows of fjac must contain the gradients of + !! the respective functions evaluated at x. + real(wp), intent(out) :: Xp(n) !! an array of length n. on output when mode = 1, + !! xp is set to a neighboring point of x. + real(wp), intent(in) :: Fvecp(m) !! an array of length m. on input when mode = 2, + !! fvecp must contain the functions evaluated at xp. + real(wp), intent(out) :: Err(m) !! an array of length m. on output when mode = 2, + !! err contains measures of correctness of the respective + !! gradients. if there is no severe loss of significance, + !! then if err(i) is 1.0 the i-th gradient is correct, + !! while if err(i) is 0.0 the i-th gradient is incorrect. + !! for values of err between 0.0 and 1.0, the categorization + !! is less certain. in general, a value of err(i) greater + !! than 0.5 indicates that the i-th gradient is probably + !! correct, while a value of err(i) less than 0.5 indicates + !! that the i-th gradient is probably incorrect. + + call chkder_(m, n, x, Fvec, Fjac, Ldfjac, Xp, Fvecp, Mode, Err) + end subroutine chkder +!***************************************************************************************** + + +!***************************************************************************************** +!> +! the purpose of hybrd is to find a zero of a system of +! n nonlinear functions in n variables by a modification +! of the powell hybrid method. the user must provide a +! subroutine which calculates the functions. the jacobian is +! then calculated by a forward-difference approximation. + + subroutine 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) + procedure(func) :: fcn !! user-supplied subroutine which calculates the functions + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of functions and variables. + integer, intent(in) :: maxfev !! a positive integer input variable. termination + !! occurs when the number of calls to `fcn` is at least `maxfev` + !! by the end of an iteration. + integer, intent(in) :: ml !! a nonnegative integer input variable which specifies + !! the number of subdiagonals within the band of the + !! jacobian matrix. if the jacobian is not banded, set + !! `ml` to at least `n - 1`. + integer, intent(in) :: mu !! a nonnegative integer input variable which specifies + !! the number of superdiagonals within the band of the + !! jacobian matrix. if the jacobian is not banded, set + !! `mu` to at least` n - 1`. + integer, intent(in) :: mode !! if `mode = 1`, the + !! variables will be scaled internally. if `mode = 2`, + !! the scaling is specified by the input `diag`. other + !! values of `mode` are equivalent to `mode = 1`. + integer, intent(in) :: nprint !! an integer input variable that enables controlled + !! printing of iterates if it is positive. in this case, + !! `fcn` is called with `iflag = 0` at the beginning of the first + !! iteration and every `nprint` iterations thereafter and + !! immediately prior to return, with `x` and `fvec` available + !! for printing. if `nprint` is not positive, no special calls + !! of `fcn` with `iflag = 0` are made. + integer, intent(out) :: info !! an integer output variable. if the user has + !! terminated execution, `info` is set to the (negative) + !! value of `iflag`. see description of `fcn`. otherwise, + !! `info` is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** relative error between two consecutive iterates + !! is at most `xtol`. + !! * ***info = 2*** number of calls to `fcn` has reached or exceeded + !! `maxfev`. + !! * ***info = 3*** `xtol` is too small. no further improvement in + !! the approximate solution `x` is possible. + !! * ***info = 4*** iteration is not making good progress, as + !! measured by the improvement from the last + !! five jacobian evaluations. + !! * ***info = 5*** iteration is not making good progress, as + !! measured by the improvement from the last + !! ten iterations. + integer, intent(out) :: nfev !! output variable set to the number of calls to `fcn`. + integer, intent(in):: ldfjac !! a positive integer input variable not less than `n` + !! which specifies the leading dimension of the array `fjac`. + integer, intent(in) :: lr !! a positive integer input variable not less than `(n*(n+1))/2`. + real(wp), intent(in) :: xtol !! a nonnegative input variable. termination + !! occurs when the relative error between two consecutive + !! iterates is at most `xtol`. + real(wp), intent(in) :: epsfcn !! an input variable used in determining a suitable + !! step length for the forward-difference approximation. this + !! approximation assumes that the relative errors in the + !! functions are of the order of `epsfcn`. if `epsfcn` is less + !! than the machine precision, it is assumed that the relative + !! errors in the functions are of the order of the machine + !! precision. + real(wp), intent(in) :: factor !! a positive input variable used in determining the + !! initial step bound. this bound is set to the product of + !! `factor` and the euclidean norm of `diag*x` if nonzero, or else + !! to `factor` itself. in most cases factor should lie in the + !! interval (.1,100.). 100. is a generally recommended value. + real(wp), intent(inout) :: x(n) !! array of length n. on input `x` must contain + !! an initial estimate of the solution vector. on output `x` + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: fvec(n) !! an output array of length `n` which contains + !! the functions evaluated at the output `x`. + real(wp), intent(inout) :: diag(n) !! an array of length `n`. if `mode = 1` (see + !! below), `diag` is internally set. if `mode = 2`, `diag` + !! must contain positive entries that serve as + !! multiplicative scale factors for the variables. + real(wp), intent(out) :: fjac(ldfjac, n) !! array which contains the + !! orthogonal matrix `q` produced by the QR factorization + !! of the final approximate jacobian. + real(wp), intent(out) :: r(lr) !! an output array which contains the + !! upper triangular matrix produced by the QR factorization + !! of the final approximate jacobian, stored rowwise. + real(wp), intent(out) :: qtf(n) !! an output array of length `n` which contains + !! the vector `(q transpose)*fvec`. + real(wp), intent(inout) :: wa1(n) !! work array + real(wp), intent(inout) :: wa2(n) !! work array + real(wp), intent(inout) :: wa3(n) !! work array + real(wp), intent(inout) :: wa4(n) !! work array + + type(hybrd_data) :: wrapper + wrapper%fcn => fcn + + call hybrd_(wrap_func, n, x, Fvec, Xtol, Maxfev, Ml, Mu, Epsfcn, Diag, Mode, & + Factor, Nprint, Info, Nfev, Fjac, Ldfjac, r, Lr, Qtf, Wa1, & + Wa2, Wa3, Wa4, wrapper) + end subroutine hybrd +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of hybrd1 is to find a zero of a system of +! n nonlinear functions in n variables by a modification +! of the powell hybrid method. this is done by using the +! more general nonlinear equation solver hybrd. the user +! must provide a subroutine which calculates the functions. +! the jacobian is then calculated by a forward-difference +! approximation. + + subroutine hybrd1(fcn, n, x, Fvec, Tol, Info, Wa, Lwa) + + implicit none + + procedure(func) :: fcn !! user-supplied subroutine which calculates the functions + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of functions and variables. + integer, intent(out) :: info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of `iflag`. see description of `fcn`. otherwise, + !! `info` is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** algorithm estimates that the relative error + !! between `x` and the solution is at most `tol`. + !! * ***info = 2*** number of calls to `fcn` has reached or exceeded + !! `200*(n+1)`. + !! * ***info = 3*** `tol` is too small. no further improvement in + !! the approximate solution `x` is possible. + !! * ***info = 4*** iteration is not making good progress. + real(wp), intent(in) :: tol !! a nonnegative input variable. termination occurs + !! when the algorithm estimates that the relative error + !! between `x` and the solution is at most `tol`. + real(wp), dimension(n), intent(inout) :: x !! an array of length `n`. on input `x` must contain + !! an initial estimate of the solution vector. on output `x` + !! contains the final estimate of the solution vector. + real(wp), dimension(n), intent(out) :: fvec !! an output array of length `n` which contains + !! the functions evaluated at the output `x`. + integer, intent(in) :: Lwa !! a positive integer input variable not less than + !! (n*(3*n+13))/2. + real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + + type(hybrd_data) :: wrapper + wrapper%fcn => fcn + + call hybrd1_(wrap_func, n, x, Fvec, Tol, Info, Wa, Lwa, wrapper) + end subroutine hybrd1 +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of hybrj is to find a zero of a system of +! n nonlinear functions in n variables by a modification +! of the powell hybrid method. the user must provide a +! subroutine which calculates the functions and the jacobian. + + subroutine hybrj(fcn, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & + Factor, Nprint, Info, Nfev, Njev, r, Lr, Qtf, Wa1, Wa2, & + Wa3, Wa4) + + implicit none + + procedure(fcn_hybrj) :: fcn !! the user-supplied subroutine which + !! calculates the functions and the jacobian + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of functions and variables. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than n + !! which specifies the leading dimension of the array fjac. + integer, intent(in) :: Maxfev !! a positive integer input variable. termination + !! occurs when the number of calls to fcn with iflag = 1 + !! has reached maxfev. + integer, intent(in) :: Mode !! an integer input variable. if mode = 1, the + !! variables will be scaled internally. if mode = 2, + !! the scaling is specified by the input diag. other + !! values of mode are equivalent to mode = 1. + integer, intent(in) :: Nprint !! an integer input variable that enables controlled + !! printing of iterates if it is positive. in this case, + !! fcn is called with iflag = 0 at the beginning of the first + !! iteration and every nprint iterations thereafter and + !! immediately prior to return, with x and fvec available + !! for printing. fvec and fjac should not be altered. + !! if nprint is not positive, no special calls of fcn + !! with iflag = 0 are made. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** relative error between two consecutive iterates + !! is at most xtol. + !! * ***info = 2*** number of calls to fcn with iflag = 1 has + !! reached maxfev. + !! * ***info = 3*** xtol is too small. no further improvement in + !! the approximate solution x is possible. + !! * ***info = 4*** iteration is not making good progress, as + !! measured by the improvement from the last + !! five jacobian evaluations. + !! * ***info = 5*** iteration is not making good progress, as + !! measured by the improvement from the last + !! ten iterations. + integer, intent(out) :: Nfev !! an integer output variable set to the number of + !! calls to fcn with iflag = 1. + integer, intent(out) :: Njev !! an integer output variable set to the number of + !! calls to fcn with iflag = 2. + integer, intent(in) :: Lr !! a positive integer input variable not less than + !! (n*(n+1))/2. + real(wp), intent(in) :: Xtol !! a nonnegative input variable. termination + !! occurs when the relative error between two consecutive + !! iterates is at most xtol. + real(wp), intent(in) :: Factor !! a positive input variable used in determining the + !! initial step bound. this bound is set to the product of + !! factor and the euclidean norm of diag*x if nonzero, or else + !! to factor itself. in most cases factor should lie in the + !! interval (.1,100.). 100. is a generally recommended value. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(n) !! an output array of length n which contains + !! the functions evaluated at the output x. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output n by n array which contains the + !! orthogonal matrix q produced by the qr factorization + !! of the final approximate jacobian. + real(wp), intent(inout) :: Diag(n) !! an array of length n. if mode = 1 (see + !! below), diag is internally set. if mode = 2, diag + !! must contain positive entries that serve as + !! multiplicative scale factors for the variables. + real(wp), intent(out) :: r(Lr) !! an output array of length lr which contains the + !! upper triangular matrix produced by the qr factorization + !! of the final approximate jacobian, stored rowwise. + real(wp), intent(out) :: Qtf(n) !! an output array of length n which contains + !! the vector (q transpose)*fvec. + real(wp), intent(inout) :: Wa1(n) !! work array of length n. + real(wp), intent(inout) :: Wa2(n) !! work array of length n. + real(wp), intent(inout) :: Wa3(n) !! work array of length n. + real(wp), intent(inout) :: Wa4(n) !! work array of length n. + + type(hybrj_data) :: wrapper + wrapper%fcn => fcn + + call hybrj_(wrap_fcn_hybrj, n, x, Fvec, Fjac, Ldfjac, Xtol, Maxfev, Diag, Mode, & + Factor, Nprint, Info, Nfev, Njev, r, Lr, Qtf, Wa1, Wa2, & + Wa3, Wa4, wrapper) + end subroutine hybrj +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of hybrj1 is to find a zero of a system of +! n nonlinear functions in n variables by a modification +! of the powell hybrid method. this is done by using the +! more general nonlinear equation solver hybrj. the user +! must provide a subroutine which calculates the functions +! and the jacobian. + + subroutine hybrj1(fcn, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa) + + implicit none + + procedure(fcn_hybrj) :: fcn !! the user-supplied subroutine which + !! calculates the functions and the jacobian + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of functions and variables. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than n + !! which specifies the leading dimension of the array fjac. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** algorithm estimates that the relative error + !! between x and the solution is at most tol. + !! * ***info = 2*** number of calls to fcn with iflag = 1 has + !! reached 100*(n+1). + !! * ***info = 3*** tol is too small. no further improvement in + !! the approximate solution x is possible. + !! * ***info = 4*** iteration is not making good progress. + integer, intent(in) :: Lwa !! a positive integer input variable not less than + !! (n*(n+13))/2. + real(wp), intent(in) :: Tol !! a nonnegative input variable. termination occurs + !! when the algorithm estimates that the relative error + !! between x and the solution is at most tol. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(n) !! an output array of length n which contains + !! the functions evaluated at the output x. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output n by n array which contains the + !! orthogonal matrix q produced by the qr factorization + !! of the final approximate jacobian. + real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + + type(hybrj_data) :: wrapper + wrapper%fcn => fcn + + call hybrj1_(wrap_fcn_hybrj, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Wa, Lwa, wrapper) + end subroutine hybrj1 +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of lmder is to minimize the sum of the squares of +! m nonlinear functions in n variables by a modification of +! the levenberg-marquardt algorithm. the user must provide a +! subroutine which calculates the functions and the jacobian. + + subroutine lmder(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & + Diag, Mode, Factor, Nprint, Info, Nfev, Njev, Ipvt, Qtf, & + Wa1, Wa2, Wa3, Wa4) + + implicit none + + procedure(fcn_lmder) :: fcn !! the user-supplied subroutine which + !! calculates the functions and the jacobian + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. n must not exceed m. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than m + !! which specifies the leading dimension of the array fjac. + integer, intent(in) :: Maxfev !! a positive integer input variable. termination + !! occurs when the number of calls to fcn with iflag = 1 + !! has reached maxfev. + integer, intent(in) :: Mode !! an integer input variable. if mode = 1, the + !! variables will be scaled internally. if mode = 2, + !! the scaling is specified by the input diag. other + !! values of mode are equivalent to mode = 1. + integer, intent(in) :: Nprint !! an integer input variable that enables controlled + !! printing of iterates if it is positive. in this case, + !! fcn is called with iflag = 0 at the beginning of the first + !! iteration and every nprint iterations thereafter and + !! immediately prior to return, with x, fvec, and fjac + !! available for printing. fvec and fjac should not be + !! altered. if nprint is not positive, no special calls + !! of fcn with iflag = 0 are made. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** both actual and predicted relative reductions + !! in the sum of squares are at most ftol. + !! * ***info = 2*** relative error between two consecutive iterates + !! is at most xtol. + !! * ***info = 3*** conditions for info = 1 and info = 2 both hold. + !! * ***info = 4*** the cosine of the angle between fvec and any + !! column of the jacobian is at most gtol in + !! absolute value. + !! * ***info = 5*** number of calls to fcn with iflag = 1 has + !! reached maxfev. + !! * ***info = 6*** ftol is too small. no further reduction in + !! the sum of squares is possible. + !! * ***info = 7*** xtol is too small. no further improvement in + !! the approximate solution x is possible. + !! * ***info = 8*** gtol is too small. fvec is orthogonal to the + !! columns of the jacobian to machine precision. + integer, intent(out) :: Nfev !! an integer output variable set to the number of + !! calls to fcn with iflag = 1. + integer, intent(out) :: Njev !! an integer output variable set to the number of + !! calls to fcn with iflag = 2. + integer, intent(out) :: Ipvt(n) !! an integer output array of length n. ipvt + !! defines a permutation matrix p such that jac*p = q*r, + !! where jac is the final calculated jacobian, q is + !! orthogonal (not stored), and r is upper triangular + !! with diagonal elements of nonincreasing magnitude. + !! column j of p is column ipvt(j) of the identity matrix. + real(wp), intent(in) :: Ftol !! a nonnegative input variable. termination + !! occurs when both the actual and predicted relative + !! reductions in the sum of squares are at most ftol. + !! therefore, ftol measures the relative error desired + !! in the sum of squares. + real(wp), intent(in) :: Xtol !! a nonnegative input variable. termination + !! occurs when the relative error between two consecutive + !! iterates is at most xtol. therefore, xtol measures the + !! relative error desired in the approximate solution. + real(wp), intent(in) :: Gtol !! a nonnegative input variable. termination + !! occurs when the cosine of the angle between fvec and + !! any column of the jacobian is at most gtol in absolute + !! value. therefore, gtol measures the orthogonality + !! desired between the function vector and the columns + !! of the jacobian. + real(wp), intent(in) :: Factor !! a positive input variable used in determining the + !! initial step bound. this bound is set to the product of + !! factor and the euclidean norm of diag*x if nonzero, or else + !! to factor itself. in most cases factor should lie in the + !! interval (.1,100.).100. is a generally recommended value. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains + !! the functions evaluated at the output x. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output m by n array. the upper n by n submatrix + !! of fjac contains an upper triangular matrix r with + !! diagonal elements of nonincreasing magnitude such that + !!``` + !! t t t + !! p *(jac *jac)*p = r *r, + !!``` + !! where p is a permutation matrix and jac is the final + !! calculated jacobian. column j of p is column ipvt(j) + !! (see below) of the identity matrix. the lower trapezoidal + !! part of fjac contains information generated during + !! the computation of r. + real(wp), intent(inout) :: Diag(n) !! an array of length n. if mode = 1 (see + !! below), diag is internally set. if mode = 2, diag + !! must contain positive entries that serve as + !! multiplicative scale factors for the variables. + real(wp), intent(out) :: Qtf(n) !! an output array of length n which contains + !! the first n elements of the vector (q transpose)*fvec. + real(wp), intent(inout) :: Wa1(n) !! work array of length n. + real(wp), intent(inout) :: Wa2(n) !! work array of length n. + real(wp), intent(inout) :: Wa3(n) !! work array of length n. + real(wp), intent(inout) :: Wa4(m) !! work array of length m. + + type(lmder_data) :: wrapper + wrapper%fcn => fcn + + call lmder_(wrap_fcn_lmder, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & + Diag, Mode, Factor, Nprint, Info, Nfev, Njev, Ipvt, Qtf, & + Wa1, Wa2, Wa3, Wa4, wrapper) + end subroutine lmder +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of lmder1 is to minimize the sum of the squares of +! m nonlinear functions in n variables by a modification of the +! levenberg-marquardt algorithm. this is done by using the more +! general least-squares solver lmder. the user must provide a +! subroutine which calculates the functions and the jacobian. + + subroutine lmder1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) + implicit none + + procedure(fcn_lmder) :: fcn !! user-supplied subroutine which + !! calculates the functions and the jacobian. + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. n must not exceed m. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than m + !! which specifies the leading dimension of the array fjac. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows. + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** algorithm estimates that the relative error + !! in the sum of squares is at most tol. + !! * ***info = 2*** algorithm estimates that the relative error + !! between x and the solution is at most tol. + !! * ***info = 3*** conditions for info = 1 and info = 2 both hold. + !! * ***info = 4*** fvec is orthogonal to the columns of the + !! jacobian to machine precision. + !! * ***info = 5*** number of calls to fcn with iflag = 1 has + !! reached 100*(n+1). + !! * ***info = 6*** tol is too small. no further reduction in + !! the sum of squares is possible. + !! * ***info = 7*** tol is too small. no further improvement in + !! the approximate solution x is possible. + integer, intent(in) :: Lwa !! a positive integer input variable not less than 5*n+m. + integer, intent(out) :: Ipvt(n) !! an integer output array of length n. ipvt + !! defines a permutation matrix p such that jac*p = q*r, + !! where jac is the final calculated jacobian, q is + !! orthogonal (not stored), and r is upper triangular + !! with diagonal elements of nonincreasing magnitude. + !! column j of p is column ipvt(j) of the identity matrix. + real(wp), intent(in) :: Tol !! a nonnegative input variable. termination occurs + !! when the algorithm estimates either that the relative + !! error in the sum of squares is at most tol or that + !! the relative error between x and the solution is at + !! most tol. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains + !! the functions evaluated at the output x. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output m by n array. the upper n by n submatrix + !! of fjac contains an upper triangular matrix r with + !! diagonal elements of nonincreasing magnitude such that + !!``` + !! t t t + !! p *(jac *jac)*p = r *r, + !!``` + !! where p is a permutation matrix and jac is the final + !! calculated jacobian. column j of p is column ipvt(j) + !! (see below) of the identity matrix. the lower trapezoidal + !! part of fjac contains information generated during + !! the computation of r. + real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + + type(lmder_data) :: wrapper + wrapper%fcn => fcn + + call lmder1_(wrap_fcn_lmder, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, & + & wrapper) + end subroutine lmder1 +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of lmdif is to minimize the sum of the squares of +! m nonlinear functions in n variables by a modification of +! the levenberg-marquardt algorithm. the user must provide a +! subroutine which calculates the functions. the jacobian is +! then calculated by a forward-difference approximation. + + subroutine lmdif(fcn, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & + Mode, Factor, Nprint, Info, Nfev, Fjac, Ldfjac, Ipvt, & + Qtf, Wa1, Wa2, Wa3, Wa4) + implicit none + + procedure(func2) :: fcn !! the user-supplied subroutine which + !! calculates the functions. + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. n must not exceed m. + integer, intent(in) :: Maxfev !! a positive integer input variable. termination + !! occurs when the number of calls to fcn is at least + !! maxfev by the end of an iteration. + integer, intent(in) :: Mode !! an integer input variable. if mode = 1, the + !! variables will be scaled internally. if mode = 2, + !! the scaling is specified by the input diag. other + !! values of mode are equivalent to mode = 1. + integer, intent(in) :: Nprint !! an integer input variable that enables controlled + !! printing of iterates if it is positive. in this case, + !! fcn is called with iflag = 0 at the beginning of the first + !! iteration and every nprint iterations thereafter and + !! immediately prior to return, with x and fvec available + !! for printing. if nprint is not positive, no special calls + !! of fcn with iflag = 0 are made. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** both actual and predicted relative reductions + !! in the sum of squares are at most ftol. + !! * ***info = 2*** relative error between two consecutive iterates + !! is at most xtol. + !! * ***info = 3*** conditions for info = 1 and info = 2 both hold. + !! * ***info = 4*** the cosine of the angle between fvec and any + !! column of the jacobian is at most gtol in + !! absolute value. + !! * ***info = 5*** number of calls to fcn has reached or + !! exceeded maxfev. + !! * ***info = 6*** ftol is too small. no further reduction in + !! the sum of squares is possible. + !! * ***info = 7*** xtol is too small. no further improvement in + !! the approximate solution x is possible. + !! * ***info = 8*** gtol is too small. fvec is orthogonal to the + !! columns of the jacobian to machine precision. + integer, intent(out) :: Nfev !! an integer output variable set to the number of + !! calls to fcn. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than m + !! which specifies the leading dimension of the array fjac. + integer, intent(out) :: Ipvt(n) !! an integer output array of length n. ipvt + !! defines a permutation matrix p such that jac*p = q*r, + !! where jac is the final calculated jacobian, q is + !! orthogonal (not stored), and r is upper triangular + !! with diagonal elements of nonincreasing magnitude. + !! column j of p is column ipvt(j) of the identity matrix. + real(wp), intent(in) :: Ftol !! a nonnegative input variable. termination + !! occurs when both the actual and predicted relative + !! reductions in the sum of squares are at most ftol. + !! therefore, ftol measures the relative error desired + !! in the sum of squares. + real(wp), intent(in) :: Xtol !! a nonnegative input variable. termination + !! occurs when the relative error between two consecutive + !! iterates is at most xtol. therefore, xtol measures the + !! relative error desired in the approximate solution. + real(wp), intent(in) :: Gtol !! a nonnegative input variable. termination + !! occurs when the cosine of the angle between fvec and + !! any column of the jacobian is at most gtol in absolute + !! value. therefore, gtol measures the orthogonality + !! desired between the function vector and the columns + !! of the jacobian. + real(wp), intent(in) :: Epsfcn !! an input variable used in determining a suitable + !! step length for the forward-difference approximation. this + !! approximation assumes that the relative errors in the + !! functions are of the order of epsfcn. if epsfcn is less + !! than the machine precision, it is assumed that the relative + !! errors in the functions are of the order of the machine + !! precision. + real(wp), intent(in) :: Factor !! a positive input variable used in determining the + !! initial step bound. this bound is set to the product of + !! factor and the euclidean norm of diag*x if nonzero, or else + !! to factor itself. in most cases factor should lie in the + !! interval (.1,100.). 100. is a generally recommended value. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains + !! the functions evaluated at the output x. + real(wp), intent(inout) :: Diag(n) !! an array of length n. if mode = 1 (see + !! below), diag is internally set. if mode = 2, diag + !! must contain positive entries that serve as + !! multiplicative scale factors for the variables. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output m by n array. the upper n by n submatrix + !! of fjac contains an upper triangular matrix r with + !! diagonal elements of nonincreasing magnitude such that + !!``` + !! t t t + !! p *(jac *jac)*p = r *r, + !!``` + !! where p is a permutation matrix and jac is the final + !! calculated jacobian. column j of p is column ipvt(j) + !! (see below) of the identity matrix. the lower trapezoidal + !! part of fjac contains information generated during + !! the computation of r. + real(wp), intent(out) :: Qtf(n) !! an output array of length n which contains + !! the first n elements of the vector (q transpose)*fvec. + real(wp), intent(inout) :: Wa1(n) !! work array of length n. + real(wp), intent(inout) :: Wa2(n) !! work array of length n. + real(wp), intent(inout) :: Wa3(n) !! work array of length n. + real(wp), intent(inout) :: Wa4(m) !! work array of length n. + + type(lmdif_data) :: wrapper + wrapper%fcn => fcn + + call lmdif_(wrap_func2, m, n, x, Fvec, Ftol, Xtol, Gtol, Maxfev, Epsfcn, Diag, & + Mode, Factor, Nprint, Info, Nfev, Fjac, Ldfjac, Ipvt, & + Qtf, Wa1, Wa2, Wa3, Wa4, wrapper) + end subroutine lmdif +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of lmdif1 is to minimize the sum of the squares of +! m nonlinear functions in n variables by a modification of the +! levenberg-marquardt algorithm. this is done by using the more +! general least-squares solver lmdif. the user must provide a +! subroutine which calculates the functions. the jacobian is +! then calculated by a forward-difference approximation. + + subroutine lmdif1(fcn, m, n, x, Fvec, Tol, Info, Iwa, Wa, Lwa) + implicit none + + procedure(func2) :: fcn !! the user-supplied subroutine which + !! calculates the functions. + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. n must not exceed m. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** algorithm estimates that the relative error + !! in the sum of squares is at most tol. + !! * ***info = 2*** algorithm estimates that the relative error + !! between x and the solution is at most tol. + !! * ***info = 3*** conditions for info = 1 and info = 2 both hold. + !! * ***info = 4*** fvec is orthogonal to the columns of the + !! jacobian to machine precision. + !! * ***info = 5*** number of calls to fcn has reached or + !! exceeded 200*(n+1). + !! * ***info = 6*** tol is too small. no further reduction in + !! the sum of squares is possible. + !! * ***info = 7*** tol is too small. no further improvement in + !! the approximate solution x is possible. + integer, intent(in) :: Lwa !! a positive integer input variable not less than + !! m*n+5*n+m. + integer, intent(inout) :: Iwa(n) !! an integer work array of length n. + real(wp), intent(in) :: Tol !! a nonnegative input variable. termination occurs + !! when the algorithm estimates either that the relative + !! error in the sum of squares is at most tol or that + !! the relative error between x and the solution is at + !! most tol. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains + !! the functions evaluated at the output x. + real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + + type(lmdif_data) :: wrapper + wrapper%fcn => fcn + + call lmdif1_(wrap_func2, m, n, x, Fvec, Tol, Info, Iwa, Wa, Lwa, wrapper) + end subroutine lmdif1 +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of lmstr is to minimize the sum of the squares of +! m nonlinear functions in n variables by a modification of +! the levenberg-marquardt algorithm which uses minimal storage. +! the user must provide a subroutine which calculates the +! functions and the rows of the jacobian. + + subroutine lmstr(fcn, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & + Diag, Mode, Factor, Nprint, Info, Nfev, Njev, Ipvt, Qtf, & + Wa1, Wa2, Wa3, Wa4) + implicit none + + procedure(fcn_lmstr) :: fcn !! user-supplied subroutine which + !! calculates the functions and the rows of the jacobian. + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. n must not exceed m. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than n + !! which specifies the leading dimension of the array fjac. + integer, intent(in) :: Maxfev !! a positive integer input variable. termination + !! occurs when the number of calls to fcn with iflag = 1 + !! has reached maxfev. + integer, intent(in) :: Mode !! an integer input variable. if mode = 1, the + !! variables will be scaled internally. if mode = 2, + !! the scaling is specified by the input diag. other + !! values of mode are equivalent to mode = 1. + integer, intent(in) :: Nprint !! an integer input variable that enables controlled + !! printing of iterates if it is positive. in this case, + !! fcn is called with iflag = 0 at the beginning of the first + !! iteration and every nprint iterations thereafter and + !! immediately prior to return, with x and fvec available + !! for printing. if nprint is not positive, no special calls + !! of fcn with iflag = 0 are made. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** both actual and predicted relative reductions + !! in the sum of squares are at most ftol. + !! * ***info = 2*** relative error between two consecutive iterates + !! is at most xtol. + !! * ***info = 3*** conditions for info = 1 and info = 2 both hold. + !! * ***info = 4*** the cosine of the angle between fvec and any + !! column of the jacobian is at most gtol in + !! absolute value. + !! * ***info = 5*** number of calls to fcn with iflag = 1 has + !! reached maxfev. + !! * ***info = 6*** ftol is too small. no further reduction in + !! the sum of squares is possible. + !! * ***info = 7*** xtol is too small. no further improvement in + !! the approximate solution x is possible. + !! * ***info = 8*** gtol is too small. fvec is orthogonal to the + !! columns of the jacobian to machine precision. + integer, intent(out) :: Nfev !! an integer output variable set to the number of + !! calls to fcn with iflag = 1. + integer, intent(out) :: Njev !! an integer output variable set to the number of + !! calls to fcn with iflag = 2. + integer, intent(out) :: Ipvt(n) !! an integer output array of length n. ipvt + !! defines a permutation matrix p such that jac*p = q*r, + !! where jac is the final calculated jacobian, q is + !! orthogonal (not stored), and r is upper triangular. + !! column j of p is column ipvt(j) of the identity matrix. + real(wp), intent(in) :: Ftol !! a nonnegative input variable. termination + !! occurs when both the actual and predicted relative + !! reductions in the sum of squares are at most ftol. + !! therefore, ftol measures the relative error desired + !! in the sum of squares. + real(wp), intent(in) :: Xtol !! a nonnegative input variable. termination + !! occurs when the relative error between two consecutive + !! iterates is at most xtol. therefore, xtol measures the + !! relative error desired in the approximate solution. + real(wp), intent(in) :: Gtol !! a nonnegative input variable. termination + !! occurs when the cosine of the angle between fvec and + !! any column of the jacobian is at most gtol in absolute + !! value. therefore, gtol measures the orthogonality + !! desired between the function vector and the columns + !! of the jacobian. + real(wp), intent(in) :: Factor !! a positive input variable used in determining the + !! initial step bound. this bound is set to the product of + !! factor and the euclidean norm of diag*x if nonzero, or else + !! to factor itself. in most cases factor should lie in the + !! interval (.1,100.). 100. is a generally recommended value. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains + !! the functions evaluated at the output x. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output n by n array. the upper triangle of fjac + !! contains an upper triangular matrix r such that + !!``` + !! t t t + !! p *(jac *jac)*p = r *r, + !!``` + !! where p is a permutation matrix and jac is the final + !! calculated jacobian. column j of p is column ipvt(j) + !! (see below) of the identity matrix. the lower triangular + !! part of fjac contains information generated during + !! the computation of r. + real(wp), intent(inout) :: Diag(n) !! an array of length n. if mode = 1 (see + !! below), diag is internally set. if mode = 2, diag + !! must contain positive entries that serve as + !! multiplicative scale factors for the variables. + real(wp), intent(out) :: Qtf(n) !! an output array of length n which contains + !! the first n elements of the vector (q transpose)*fvec. + real(wp), intent(inout) :: Wa1(n) !! work array of length n. + real(wp), intent(inout) :: Wa2(n) !! work array of length n. + real(wp), intent(inout) :: Wa3(n) !! work array of length n. + real(wp), intent(inout) :: Wa4(m) !! work array of length m. + + type(lmstr_data) :: wrapper + wrapper%fcn => fcn + + call lmstr_(wrap_fcn_lmstr, m, n, x, Fvec, Fjac, Ldfjac, Ftol, Xtol, Gtol, Maxfev, & + Diag, Mode, Factor, Nprint, Info, Nfev, Njev, Ipvt, Qtf, & + Wa1, Wa2, Wa3, Wa4, wrapper) + end subroutine lmstr +!***************************************************************************************** + +!***************************************************************************************** +!> +! the purpose of lmstr1 is to minimize the sum of the squares of +! m nonlinear functions in n variables by a modification of +! the levenberg-marquardt algorithm which uses minimal storage. +! this is done by using the more general least-squares solver +! lmstr. the user must provide a subroutine which calculates +! the functions and the rows of the jacobian. + + subroutine lmstr1(fcn, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa) + implicit none + + procedure(fcn_lmstr) :: fcn !! user-supplied subroutine which + !! calculates the functions and the rows of the jacobian. + integer, intent(in) :: m !! a positive integer input variable set to the number + !! of functions. + integer, intent(in) :: n !! a positive integer input variable set to the number + !! of variables. n must not exceed m. + integer, intent(in) :: Ldfjac !! a positive integer input variable not less than n + !! which specifies the leading dimension of the array fjac. + integer, intent(out) :: Info !! an integer output variable. if the user has + !! terminated execution, info is set to the (negative) + !! value of iflag. see description of fcn. otherwise, + !! info is set as follows: + !! + !! * ***info = 0*** improper input parameters. + !! * ***info = 1*** algorithm estimates that the relative error + !! in the sum of squares is at most tol. + !! * ***info = 2*** algorithm estimates that the relative error + !! between x and the solution is at most tol. + !! * ***info = 3*** conditions for info = 1 and info = 2 both hold. + !! * ***info = 4*** fvec is orthogonal to the columns of the + !! jacobian to machine precision. + !! * ***info = 5*** number of calls to fcn with iflag = 1 has + !! reached 100*(n+1). + !! * ***info = 6*** tol is too small. no further reduction in + !! the sum of squares is possible. + !! * ***info = 7*** tol is too small. no further improvement in + !! the approximate solution x is possible. + integer, intent(in) :: Lwa !! a positive integer input variable not less than 5*n+m. + integer, intent(out) :: Ipvt(n) !! an integer output array of length n. ipvt + !! defines a permutation matrix p such that jac*p = q*r, + !! where jac is the final calculated jacobian, q is + !! orthogonal (not stored), and r is upper triangular. + !! column j of p is column ipvt(j) of the identity matrix. + real(wp), intent(in) :: Tol !! a nonnegative input variable. termination occurs + !! when the algorithm estimates either that the relative + !! error in the sum of squares is at most tol or that + !! the relative error between x and the solution is at + !! most tol. + real(wp), intent(inout) :: x(n) !! an array of length n. on input x must contain + !! an initial estimate of the solution vector. on output x + !! contains the final estimate of the solution vector. + real(wp), intent(out) :: Fvec(m) !! an output array of length m which contains + !! the functions evaluated at the output x. + real(wp), intent(out) :: Fjac(Ldfjac, n) !! an output n by n array. the upper triangle of fjac + !! contains an upper triangular matrix r such that + !!``` + !! t t t + !! p *(jac *jac)*p = r *r, + !!``` + !! where p is a permutation matrix and jac is the final + !! calculated jacobian. column j of p is column ipvt(j) + !! (see below) of the identity matrix. the lower triangular + !! part of fjac contains information generated during + !! the computation of r. + real(wp), intent(inout) :: Wa(Lwa) !! a work array of length lwa. + + type(lmstr_data) :: wrapper + wrapper%fcn => fcn + + call lmstr1_(wrap_fcn_lmstr, m, n, x, Fvec, Fjac, Ldfjac, Tol, Info, Ipvt, Wa, Lwa, & + & wrapper) + end subroutine lmstr1 +!***************************************************************************************** + + + subroutine wrap_func(n, x, fvec, iflag, wrapper) + !! user-supplied subroutine for [[hybrd]], [[hybrd1]], and [[fdjac1]] + integer, intent(in) :: n !! the number of variables. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(out) :: fvec(n) !! value of function at `x` + integer, intent(inout) :: iflag !! set to <0 to terminate execution + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(hybrd_data) + call wrapper%fcn(n, x, fvec, iflag) + class default + iflag = -1 + end select + end subroutine wrap_func + + subroutine wrap_func2(m, n, x, fvec, iflag, wrapper) + !! user-supplied subroutine for [[fdjac2]], [[lmdif]], and [[lmdif1]] + integer, intent(in) :: m !! the number of functions. + integer, intent(in) :: n !! the number of variables. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(out) :: fvec(m) !! value of function at `x` + integer, intent(inout) :: iflag !! the value of iflag should not be changed unless + !! the user wants to terminate execution of lmdif. + !! in this case set iflag to a negative integer. + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(lmdif_data) + call wrapper%fcn(m, n, x, fvec, iflag) + class default + iflag = -1 + end select + end subroutine wrap_func2 + + subroutine wrap_fcn_hybrj(n, x, fvec, fjac, ldfjac, iflag, wrapper) + !! user-supplied subroutine for [[hybrj]] and [[hybrj1]] + integer, intent(in) :: n !! the number of variables. + real(wp), dimension(n), intent(in) :: x !! independent variable vector + integer, intent(in) :: ldfjac !! leading dimension of the array fjac. + real(wp), dimension(n), intent(inout) :: fvec !! value of function at `x` + real(wp), dimension(ldfjac, n), intent(inout) :: fjac !! jacobian matrix at `x` + integer, intent(inout) :: iflag !! if iflag = 1 calculate the functions at x and + !! return this vector in fvec. do not alter fjac. + !! if iflag = 2 calculate the jacobian at x and + !! return this matrix in fjac. do not alter fvec. + !! + !! the value of iflag should not be changed by fcn unless + !! the user wants to terminate execution of hybrj. + !! in this case set iflag to a negative integer. + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(hybrj_data) + call wrapper%fcn(n, x, fvec, fjac, ldfjac, iflag) + class default + iflag = -1 + end select + end subroutine wrap_fcn_hybrj + + subroutine wrap_fcn_lmder(m, n, x, fvec, fjac, ldfjac, iflag, wrapper) + !! user-supplied subroutine for [[lmder]] and [[lmder1]] + integer, intent(in) :: m !! the number of functions. + integer, intent(in) :: n !! the number of variables. + integer, intent(in) :: ldfjac !! leading dimension of the array fjac. + integer, intent(inout) :: iflag !! if iflag = 1 calculate the functions at x and + !! return this vector in fvec. do not alter fjac. + !! if iflag = 2 calculate the jacobian at x and + !! return this matrix in fjac. do not alter fvec. + !! + !! the value of iflag should not be changed by fcn unless + !! the user wants to terminate execution of lmder. + !! in this case set iflag to a negative integer. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(inout) :: fvec(m) !! value of function at `x` + real(wp), intent(inout) :: fjac(ldfjac, n) !! jacobian matrix at `x` + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(lmder_data) + call wrapper%fcn(m, n, x, fvec, fjac, ldfjac, iflag) + class default + iflag = -1 + end select + end subroutine wrap_fcn_lmder + + subroutine wrap_fcn_lmstr(m, n, x, fvec, fjrow, iflag, wrapper) + integer, intent(in) :: m !! the number of functions. + integer, intent(in) :: n !! the number of variables. + integer, intent(inout) :: iflag !! if iflag = 1 calculate the functions at x and + !! return this vector in fvec. + !! if iflag = i calculate the (i-1)-st row of the + !! jacobian at x and return this vector in fjrow. + !! + !! the value of iflag should not be changed by fcn unless + !! the user wants to terminate execution of lmstr. + !! in this case set iflag to a negative integer. + real(wp), intent(in) :: x(n) !! independent variable vector + real(wp), intent(inout) :: fvec(m) !! value of function at `x` + real(wp), intent(inout) :: fjrow(n) !! jacobian row + class(*), intent(inout), optional :: wrapper + + if (.not.present(wrapper)) then + iflag = -1 + return + end if + + select type(wrapper) + class is(lmstr_data) + call wrapper%fcn(m, n, x, fvec, fjrow, iflag) + class default + iflag = -1 + end select + end subroutine wrap_fcn_lmstr + +!***************************************************************************************** +end module minpack_legacy +!***************************************************************************************** diff --git a/test/test_hybrd.f90 b/test/test_hybrd.f90 index 9bd829d..d60c5da 100644 --- a/test/test_hybrd.f90 +++ b/test/test_hybrd.f90 @@ -12,7 +12,8 @@ program test_hybrd - use minpack_module, only: wp, dpmpar, enorm, hybrd1 + use minpack_legacy, only: hybrd1 + use minpack_module, only: wp, dpmpar, enorm use iso_fortran_env, only: output_unit implicit none @@ -698,4 +699,4 @@ end subroutine initpt !***************************************************************************************** end program test_hybrd -!***************************************************************************************** \ No newline at end of file +!***************************************************************************************** diff --git a/test/test_hybrj.f90 b/test/test_hybrj.f90 index 2955cee..235e2ee 100644 --- a/test/test_hybrj.f90 +++ b/test/test_hybrj.f90 @@ -12,7 +12,8 @@ program test_hybrj - use minpack_module, only: wp, dpmpar, enorm, hybrj1 + use minpack_legacy, only: hybrj1 + use minpack_module, only: wp, dpmpar, enorm use iso_fortran_env, only: nwrite => output_unit implicit none diff --git a/test/test_lmder.f90 b/test/test_lmder.f90 index 2ed31ac..99e63df 100644 --- a/test/test_lmder.f90 +++ b/test/test_lmder.f90 @@ -12,7 +12,8 @@ program test_lmder - use minpack_module, only: wp, dpmpar, enorm, lmder1 + use minpack_legacy, only: lmder1 + use minpack_module, only: wp, dpmpar, enorm use iso_fortran_env, only: output_unit implicit none @@ -1194,4 +1195,4 @@ end subroutine ssqfcn !***************************************************************************************** end program test_lmder -!***************************************************************************************** \ No newline at end of file +!***************************************************************************************** diff --git a/test/test_lmdif.f90 b/test/test_lmdif.f90 index 978aae5..54c35c1 100644 --- a/test/test_lmdif.f90 +++ b/test/test_lmdif.f90 @@ -12,7 +12,8 @@ program test_lmdif - use minpack_module, only: wp, dpmpar, enorm, lmdif1 + use minpack_legacy, only: lmdif1 + use minpack_module, only: wp, dpmpar, enorm use iso_fortran_env, only: nwrite => output_unit implicit none @@ -753,4 +754,4 @@ end subroutine initpt !***************************************************************************************** end program test_lmdif -!***************************************************************************************** \ No newline at end of file +!***************************************************************************************** diff --git a/test/test_lmstr.f90 b/test/test_lmstr.f90 index 908b798..eec09f4 100644 --- a/test/test_lmstr.f90 +++ b/test/test_lmstr.f90 @@ -12,7 +12,8 @@ program test_lmstr - use minpack_module, only: wp, dpmpar, enorm, lmstr1 + use minpack_legacy, only : lmstr1 + use minpack_module, only: wp, dpmpar, enorm use iso_fortran_env, only: nwrite => output_unit implicit none @@ -1062,4 +1063,4 @@ end subroutine ssqfcn !***************************************************************************************** end program test_lmstr -!***************************************************************************************** \ No newline at end of file +!*****************************************************************************************