diff --git a/fvn_misc/fvn_misc.f90 b/fvn_misc/fvn_misc.f90 index 9c043e7..b193be3 100644 --- a/fvn_misc/fvn_misc.f90 +++ b/fvn_misc/fvn_misc.f90 @@ -7,6 +7,11 @@ interface fvn_muller module procedure fvn_z_muller end interface fvn_muller +public +private :: lmdif,lmpar,qrfac,qrsolv,enorm,fdjac2 +! These are made private to not interfere with a +! possibly linked minpack + contains ! ! Muller @@ -320,31 +325,94 @@ icloop: do ! This uses MINPACK Routines (http://www.netlib.org/minpack) ! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au ! +! The purpose of fvn_lm 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. +! +! call fvn_lm(fcn,m,n,a,info,tol) +! +! fcn : fcn is the user-supplied subroutine which calculates +! the functions. fcn must follow the following interface that must +! be declared in the calling subroutine : +! +! interface +! subroutine fcn(m,n,a,fvec,iflag) +! use fvn_common +! integer(ip_kind), intent(in) :: m +! integer(ip_kind), intent(in) :: n +! real(dp_kind), dimension(:), intent(in) :: a +! real(dp_kind), dimension(:), intent(inout) :: fvec +! integer(ip_kind), intent(inout) :: iflag +! end subroutine +! end interface +! +! This is the function which calculate the differences for which which square sum +! will be minimized outputing this difference in vector fvec. +! Parameters of fcn are as follows : +! m : positive integer input variable set to the number of functions +! (number of measurement points) +! n : positive integer input variable set to the number of variables +! (number of parameters in the function to fit) +! a : vector of length n containing parameters for which fcn should +! perform the calculation +! fvec : vector of length m containing the resulting evaluation +! iflag : integer normally not used, can be used to exit the +! the algorithm by setting it to a negative value +! +! m : positive integer input variable set to the number of functions +! (number of measurement points) +! n : positive integer input variable set to the number of variables +! (number of parameters in the function to fit) +! a : vector of length n, on input must contains an initial guess (or unity vector) +! and on output the solution vector +! info : is an output positive integer +! 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. +! tol : is an optional positive value. 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. If not provided default value is : +! sqrt(epsilon(0.0d0)) +! +subroutine fvn_lm(fcn,m,n,a,info,tol) + integer(ip_kind), intent(in) :: m + integer(ip_kind), intent(in) :: n + real(dp_kind), dimension(:), intent(inout) :: a + integer(ip_kind), intent(out) :: info + real(dp_kind), intent(in), optional :: tol -subroutine fvn_lmdif(zef,m,n,a,info,tol) - integer(4), intent(in) :: m - integer(4), intent(in) :: n - real(8), dimension(:), intent(inout) :: a - integer(4), intent(out) :: info - real(8), intent(in), optional :: tol - - real(8) :: rtol - real(8), dimension(:), allocatable :: fvec - integer(4), dimension(:), allocatable :: iwa + real(dp_kind) :: rtol + real(dp_kind), dimension(:), allocatable :: fvec + integer(ip_kind), dimension(:), allocatable :: iwa interface - subroutine zef(m,n,a,fvec,iflag) - integer(4), intent(in) :: m - integer(4), intent(in) :: n - real(8), dimension(:), intent(in) :: a - real(8), dimension(:), intent(inout) :: fvec - integer(4), intent(inout) :: iflag + subroutine fcn(m,n,a,fvec,iflag) + use fvn_common + integer(ip_kind), intent(in) :: m + integer(ip_kind), intent(in) :: n + real(dp_kind), dimension(:), intent(in) :: a + real(dp_kind), dimension(:), intent(inout) :: fvec + integer(ip_kind), intent(inout) :: iflag end subroutine end interface - integer(4) :: maxfev, mode, nfev, nprint - real(8) :: epsfcn, ftol, gtol, xtol, fjac(m,n) - real(8), parameter :: factor = 100._8, zero = 0.0_8 + integer(ip_kind) :: maxfev, mode, nfev, nprint + real(dp_kind) :: epsfcn, ftol, gtol, xtol, fjac(m,n) + real(dp_kind), parameter :: factor = 100._8, zero = 0.0_8 allocate(fvec(m),iwa(n)) @@ -367,7 +435,7 @@ subroutine fvn_lmdif(zef,m,n,a,info,tol) mode = 1 nprint = 0 - call lmdif(zef, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn, & + call lmdif(fcn, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn, & mode, factor, nprint, info, nfev, fjac, iwa) if (info == 8) info = 4 @@ -376,113 +444,14 @@ subroutine fvn_lmdif(zef,m,n,a,info,tol) end subroutine -! MINPACK routines which are used by both LMDIF & LMDER -! 25 October 2001: -! Changed INTENT of iflag in several places to IN OUT. -! Changed INTENT of fvec to IN OUT in user routine FCN. -! Removed arguments diag and qtv from LMDIF & LMDER. -! Replaced several DO loops with array operations. -! amiller @ bigpond.net.au - - -! ********** - -! subroutine lmdif1 - -! 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. - -! the subroutine statement is - -! subroutine lmdif1(fcn, m, n, x, fvec, tol, info, iwa) - -! where - -! fcn is the name of the user-supplied subroutine which calculates -! the functions. fcn must be declared in an external statement in the -! user calling program, and should be written as follows. - -! subroutine fcn(m, n, x, fvec, iflag) -! integer m, n, iflag -! REAL (dp) x(n), fvec(m) -! ---------- -! calculate the functions at x and return this vector in fvec. -! ---------- -! return -! end - -! the value of iflag should not be changed by fcn unless -! the user wants to terminate execution of lmdif1. -! In this case set iflag to a negative integer. - -! m is a positive integer input variable set to the number of functions. - -! n is a positive integer input variable set to the number of variables. -! n must not exceed m. - -! x is 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. - -! fvec is an output array of length m which contains -! the functions evaluated at the output x. - -! tol is 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. - -! info is 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. - -! iwa is an integer work array of length n. - -! wa is a work array of length lwa. - -! lwa is a positive integer input variable not less than m*n+5*n+m. - -! subprograms called - -! user-supplied ...... fcn - -! minpack-supplied ... lmdif - -! argonne national laboratory. minpack project. march 1980. -! burton s. garbow, kenneth e. hillstrom, jorge j. more - -! ********** - +! lmdif +! MINPACK Subroutine (http://www.netlib.org/minpack) +! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au SUBROUTINE lmdif(fcn, m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, & mode, factor, nprint, info, nfev, fjac, ipvt) ! N.B. Arguments LDFJAC, DIAG, QTF, WA1, WA2, WA3 & WA4 have been removed. -INTEGER, PARAMETER :: dp = 8 +INTEGER, PARAMETER :: dp = dp_kind INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: x(:) @@ -504,10 +473,11 @@ INTEGER, INTENT(OUT) :: ipvt(:) INTERFACE SUBROUTINE fcn(m, n, x, fvec, iflag) - INTEGER(4), INTENT(IN) :: m, n - REAL (8), INTENT(IN) :: x(:) - REAL (8), INTENT(IN OUT) :: fvec(:) - INTEGER(4), INTENT(IN OUT) :: iflag + use fvn_common + INTEGER(ip_kind), INTENT(IN) :: m, n + REAL (dp_kind), INTENT(IN) :: x(:) + REAL (dp_kind), INTENT(IN OUT) :: fvec(:) + INTEGER(ip_kind), INTENT(IN OUT) :: iflag END SUBROUTINE fcn END INTERFACE @@ -912,7 +882,9 @@ END SUBROUTINE lmdif ! ********** - +! lmpar +! MINPACK Subroutine (http://www.netlib.org/minpack) +! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag) ! Code converted using TO_F90 by Alan Miller @@ -920,7 +892,7 @@ SUBROUTINE lmpar(n, r, ipvt, diag, qtb, delta, par, x, sdiag) ! N.B. Arguments LDR, WA1 & WA2 have been removed. -INTEGER, PARAMETER :: dp = 8 +INTEGER, PARAMETER :: dp = dp_kind INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: r(:,:) INTEGER, INTENT(IN) :: ipvt(:) @@ -1158,14 +1130,16 @@ RETURN END SUBROUTINE lmpar - +! qrfac +! MINPACK Subroutine (http://www.netlib.org/minpack) +! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au SUBROUTINE qrfac(m, n, a, pivot, ipvt, rdiag, acnorm) ! Code converted using TO_F90 by Alan Miller ! Date: 1999-12-09 Time: 12:46:17 ! N.B. Arguments LDA, LIPVT & WA have been removed. -INTEGER, PARAMETER :: dp = 8 +INTEGER, PARAMETER :: dp = dp_kind INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: a(:,:) @@ -1320,11 +1294,13 @@ RETURN END SUBROUTINE qrfac - +! qrsolv +! MINPACK Subroutine (http://www.netlib.org/minpack) +! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au SUBROUTINE qrsolv(n, r, ipvt, diag, qtb, x, sdiag) ! N.B. Arguments LDR & WA have been removed. -INTEGER, PARAMETER :: dp = 8 +INTEGER, PARAMETER :: dp = dp_kind INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: r(:,:) INTEGER, INTENT(IN) :: ipvt(:) @@ -1507,12 +1483,14 @@ RETURN END SUBROUTINE qrsolv - +! enorm +! MINPACK Subroutine (http://www.netlib.org/minpack) +! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au FUNCTION enorm(n,x) RESULT(fn_val) ! Code converted using TO_F90 by Alan Miller ! Date: 1999-12-09 Time: 12:45:34 -INTEGER, PARAMETER :: dp = 8 +INTEGER, PARAMETER :: dp = dp_kind INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN) :: x(:) REAL (dp) :: fn_val @@ -1615,14 +1593,16 @@ GO TO 120 END FUNCTION enorm - +! fdjac2 +! MINPACK Subroutine (http://www.netlib.org/minpack) +! Converted to fortran90 by Alan Miller amiller @ bigpond.net.au SUBROUTINE fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn) ! Code converted using TO_F90 by Alan Miller ! Date: 1999-12-09 Time: 12:45:44 ! N.B. Arguments LDFJAC & WA have been removed. -INTEGER, PARAMETER :: dp = 8 +INTEGER, PARAMETER :: dp = dp_kind INTEGER, INTENT(IN) :: m INTEGER, INTENT(IN) :: n REAL (dp), INTENT(IN OUT) :: x(n) @@ -1633,10 +1613,11 @@ REAL (dp), INTENT(IN) :: epsfcn INTERFACE SUBROUTINE fcn(m, n, x, fvec, iflag) - INTEGER(4), INTENT(IN) :: m, n - REAL (8), INTENT(IN) :: x(:) - REAL (8), INTENT(IN OUT) :: fvec(:) - INTEGER(4), INTENT(IN OUT) :: iflag + use fvn_common + INTEGER(ip_kind), INTENT(IN) :: m, n + REAL (dp_kind), INTENT(IN) :: x(:) + REAL (dp_kind), INTENT(IN OUT) :: fvec(:) + INTEGER(ip_kind), INTENT(IN OUT) :: iflag END SUBROUTINE fcn END INTERFACE