module fvn_misc
use fvn_common
implicit none
! Muller
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
!
!
!
! William Daniau 2007
!
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f
!
! it can be used as a replacement for imsl routine dzanly with minor changes
!
!-----------------------------------------------------------------------
!
! purpose - zeros of an analytic complex function
! using the muller method with deflation
!
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier)
!
! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose
! roots are to be found. f must appear in
! an external statement in the calling pro-
! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given
! root agree within eps1. (input)
! note. if either or both of the stopping
! criteria are fulfilled, the root is
! accepted.
! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by
! muller (input)
! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for
! the n new roots which are to be computed.
! if the user does not provide an initial
! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller.
! itmax - the maximum allowable number of iterations
! per root (input)
! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of
! iterations used in finding the j-th root
! when convergence was achieved. if
! convergence was not obtained in itmax
! iterations, infer(j) will be greater than
! itmax (output).
! ier - error parameter (output)
! warning error
! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of
! the (n) new roots.
!
!
! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla-
! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
!
! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax
! the convergence criterion.
!
!-----------------------------------------------------------------------
!
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq,eps1w
double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5
double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n)
integer :: infer(kn+n)
data zero/(0.0d0,0.0d0)/,p1/(0.1d0,0.0d0)/, &
one/(1.0d0,0.0d0)/,four/(4.0d0,0.0d0)/, &
p5/(0.5d0,0.0d0)/, &
rzero/0.0d0/,rten/10.0d0/,rhun/100.0d0/, &
ax/0.1d0/,ickmax/3/,rp01/0.01d0/
ier = 0
if (n .lt. 1) then ! What the hell are doing here then ...
return
end if
!eps1 = rten **(-nsig)
eps1w = min(eps1,rp01)
knp1 = kn+1
knpn = kn+n
knpng = kn+nguess
do i=1,knpn
infer(i) = 0
if (i .gt. knpng) x(i) = zero
end do
l= knp1
ic=0
rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0
ick = 0
xl = x(l)
icloop: do
ic = 0
h = ax
h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are
! xl+h, xl-h, xl
rt = xl+h
call deflated_work(errcode)
if (errcode == 1) then
exit icloop
end if
z0 = fprt
y0 = frt
x0 = rt
rt = xl-h
call deflated_work(errcode)
if (errcode == 1) then
exit icloop
end if
z1 = fprt
y1 = frt
h = xl-rt
d = h/(rt-x0)
rt = xl
call deflated_work(errcode)
if (errcode == 1) then
exit icloop
end if
z2 = fprt
y2 = frt
! begin main algorithm
iloop: do
dd = one + d
t1 = z0*d*d
t2 = z1*dd*dd
xx = z2*dd
t3 = z2*d
bi = t1-t2+xx+t3
den = bi*bi-four*(xx*t1-t3*(t2-xx))
! use denominator of maximum amplitude
t1 = sqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1
tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero
t3 = bi - t1
tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero
den = t2
qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3
! test for zero denominator
if (cdabs(den) .eq. rzero) then
call trans_rt()
call deflated_work(errcode)
if (errcode == 1) then
exit icloop
end if
z2 = fprt
y2 = frt
cycle iloop
end if
d = -xx/den
d = d+d
h = d*h
rt = rt + h
! check convergence of the first kind
if (cdabs(h) .le. eps1w*max(cdabs(rt),ax)) then
if (ic .ne. 0) then
exit icloop
end if
ic = 1
z0 = y1
z1 = y2
z2 = f(rt)
xl = rt
ick = ick+1
if (ick .le. ickmax) then
cycle iloop
end if
! warning error, itmax = maximum
jk = itmax + jk
ier = 33
end if
if (ic .ne. 0) then
cycle icloop
end if
call deflated_work(errcode)
if (errcode == 1) then
exit icloop
end if
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce
! convergence
d = d*p5
h = h*p5
rt = rt-h
call deflated_work(errcode)
if (errcode == 1) then
exit icloop
end if
end do
z0 = z1
z1 = z2
z2 = fprt
y0 = y1
y1 = y2
y2 = frt
end do iloop
end do icloop
x(l) = rt
infer(l) = jk
l = l+1
end do rloop
contains
subroutine trans_rt()
tem = rten*eps1w
if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem
d = (h+tem)*d/h
h = h+tem
end subroutine trans_rt
subroutine deflated_work(errcode)
! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag
flag=1
loop1: do while(flag==1)
errcode=0
jk = jk+1
if (jk .gt. itmax) then
ier=33
errcode=1
return
end if
frt = f(rt)
fprt = frt
if (l /= 1) then
lm1 = l-1
do i=1,lm1
tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt()
cycle loop1
end if
fprt = fprt/tem
end do
end if
flag=0
end do loop1
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1
return
end if
end subroutine deflated_work
end subroutine
!
!
!
! Non linear least square using Levenberg-Marquardt algorithm and
! a finite difference jacobian
!
! 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
real(dp_kind) :: rtol
real(dp_kind), dimension(:), allocatable :: fvec
integer(ip_kind), dimension(:), allocatable :: iwa
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
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))
rtol=sqrt(epsilon(0.d0))
if (present(tol)) rtol=tol
info = 0
! check the input parameters for errors.
if (n <= 0 .or. m < n .or. rtol < zero) return
! call lmdif.
maxfev = 200*(n + 1)
ftol = rtol
xtol = rtol
gtol = zero
epsfcn = zero
mode = 1
nprint = 0
call lmdif(fcn, m, n, a, fvec, ftol, xtol, gtol, maxfev, epsfcn, &
mode, factor, nprint, info, nfev, fjac, iwa)
if (info == 8) info = 4
deallocate(fvec,iwa)
end subroutine
! 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 = dp_kind
INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: x(:)
REAL (dp), INTENT(OUT) :: fvec(:)
REAL (dp), INTENT(IN) :: ftol
REAL (dp), INTENT(IN) :: xtol
REAL (dp), INTENT(IN OUT) :: gtol
INTEGER, INTENT(IN OUT) :: maxfev
REAL (dp), INTENT(IN OUT) :: epsfcn
INTEGER, INTENT(IN) :: mode
REAL (dp), INTENT(IN) :: factor
INTEGER, INTENT(IN) :: nprint
INTEGER, INTENT(OUT) :: info
INTEGER, INTENT(OUT) :: nfev
REAL (dp), INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
INTEGER, INTENT(OUT) :: ipvt(:)
! EXTERNAL fcn
INTERFACE
SUBROUTINE fcn(m, n, x, fvec, 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
! **********
! subroutine lmdif
! 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.
! the subroutine statement is
! 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)
! N.B. 7 of these arguments have been removed in this version.
! 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(:), 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 lmdif.
! 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.
! ftol is 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.
! xtol is 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.
! gtol is 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.
! maxfev is a positive integer input variable. Termination occurs when the
! number of calls to fcn is at least maxfev by the end of an iteration.
! epsfcn is 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.
! diag is 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.
! mode is 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.
! factor is 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.
! nprint is 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.
! 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 both actual and predicted relative reductions
! in the sum of squares are at most ftol.
! info = 2 relative error between two consecutive iterates <= 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.
! nfev is an integer output variable set to the number of calls to fcn.
! fjac is 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.
! ldfjac is a positive integer input variable not less than m
! which specifies the leading dimension of the array fjac.
! ipvt is 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.
! qtf is an output array of length n which contains
! the first n elements of the vector (q transpose)*fvec.
! wa1, wa2, and wa3 are work arrays of length n.
! wa4 is a work array of length m.
! subprograms called
! user-supplied ...... fcn
! minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
! fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
! **********
INTEGER :: i, iflag, iter, j, l
REAL (dp) :: actred, delta, dirder, epsmch, fnorm, fnorm1, gnorm, &
par, pnorm, prered, ratio, sum, temp, temp1, temp2, xnorm
REAL (dp) :: diag(n), qtf(n), wa1(n), wa2(n), wa3(n), wa4(m)
REAL (dp), PARAMETER :: one = 1.0_dp, p1 = 0.1_dp, p5 = 0.5_dp, &
p25 = 0.25_dp, p75 = 0.75_dp, p0001 = 0.0001_dp, &
zero = 0.0_dp
! epsmch is the machine precision.
epsmch = EPSILON(zero)
info = 0
iflag = 0
nfev = 0
! check the input parameters for errors.
IF (n <= 0 .OR. m < n .OR. ftol < zero .OR. xtol < zero .OR. gtol < zero &
.OR. maxfev <= 0 .OR. factor <= zero) GO TO 300
IF (mode /= 2) GO TO 20
DO j = 1, n
IF (diag(j) <= zero) GO TO 300
END DO
! evaluate the function at the starting point and calculate its norm.
20 iflag = 1
CALL fcn(m, n, x, fvec, iflag)
nfev = 1
IF (iflag < 0) GO TO 300
fnorm = enorm(m, fvec)
! initialize levenberg-marquardt parameter and iteration counter.
par = zero
iter = 1
! beginning of the outer loop.
! calculate the jacobian matrix.
30 iflag = 2
CALL fdjac2(fcn, m, n, x, fvec, fjac, iflag, epsfcn)
nfev = nfev + n
IF (iflag < 0) GO TO 300
! If requested, call fcn to enable printing of iterates.
IF (nprint <= 0) GO TO 40
iflag = 0
IF (MOD(iter-1,nprint) == 0) CALL fcn(m, n, x, fvec, iflag)
IF (iflag < 0) GO TO 300
! Compute the qr factorization of the jacobian.
40 CALL qrfac(m, n, fjac, .true., ipvt, wa1, wa2)
! On the first iteration and if mode is 1, scale according
! to the norms of the columns of the initial jacobian.
IF (iter /= 1) GO TO 80
IF (mode == 2) GO TO 60
DO j = 1, n
diag(j) = wa2(j)
IF (wa2(j) == zero) diag(j) = one
END DO
! On the first iteration, calculate the norm of the scaled x
! and initialize the step bound delta.
60 wa3(1:n) = diag(1:n)*x(1:n)
xnorm = enorm(n, wa3)
delta = factor*xnorm
IF (delta == zero) delta = factor
! Form (q transpose)*fvec and store the first n components in qtf.
80 wa4(1:m) = fvec(1:m)
DO j = 1, n
IF (fjac(j,j) == zero) GO TO 120
sum = DOT_PRODUCT( fjac(j:m,j), wa4(j:m) )
temp = -sum/fjac(j,j)
DO i = j, m
wa4(i) = wa4(i) + fjac(i,j)*temp
END DO
120 fjac(j,j) = wa1(j)
qtf(j) = wa4(j)
END DO
! compute the norm of the scaled gradient.
gnorm = zero
IF (fnorm == zero) GO TO 170
DO j = 1, n
l = ipvt(j)
IF (wa2(l) == zero) CYCLE
sum = zero
DO i = 1, j
sum = sum + fjac(i,j)*(qtf(i)/fnorm)
END DO
gnorm = MAX(gnorm, ABS(sum/wa2(l)))
END DO
! test for convergence of the gradient norm.
170 IF (gnorm <= gtol) info = 4
IF (info /= 0) GO TO 300
! rescale if necessary.
IF (mode == 2) GO TO 200
DO j = 1, n
diag(j) = MAX(diag(j), wa2(j))
END DO
! beginning of the inner loop.
! determine the Levenberg-Marquardt parameter.
200 CALL lmpar(n, fjac, ipvt, diag, qtf, delta, par, wa1, wa2)
! store the direction p and x + p. calculate the norm of p.
DO j = 1, n
wa1(j) = -wa1(j)
wa2(j) = x(j) + wa1(j)
wa3(j) = diag(j)*wa1(j)
END DO
pnorm = enorm(n, wa3)
! on the first iteration, adjust the initial step bound.
IF (iter == 1) delta = MIN(delta, pnorm)
! evaluate the function at x + p and calculate its norm.
iflag = 1
CALL fcn(m, n, wa2, wa4, iflag)
nfev = nfev + 1
IF (iflag < 0) GO TO 300
fnorm1 = enorm(m, wa4)
! compute the scaled actual reduction.
actred = -one
IF (p1*fnorm1 < fnorm) actred = one - (fnorm1/fnorm)**2
! Compute the scaled predicted reduction and
! the scaled directional derivative.
DO j = 1, n
wa3(j) = zero
l = ipvt(j)
temp = wa1(l)
DO i = 1, j
wa3(i) = wa3(i) + fjac(i,j)*temp
END DO
END DO
temp1 = enorm(n,wa3)/fnorm
temp2 = (SQRT(par)*pnorm)/fnorm
prered = temp1**2 + temp2**2/p5
dirder = -(temp1**2 + temp2**2)
! compute the ratio of the actual to the predicted reduction.
ratio = zero
IF (prered /= zero) ratio = actred/prered
! update the step bound.
IF (ratio <= p25) THEN
IF (actred >= zero) temp = p5
IF (actred < zero) temp = p5*dirder/(dirder + p5*actred)
IF (p1*fnorm1 >= fnorm .OR. temp < p1) temp = p1
delta = temp*MIN(delta,pnorm/p1)
par = par/temp
ELSE
IF (par /= zero .AND. ratio < p75) GO TO 260
delta = pnorm/p5
par = p5*par
END IF
! test for successful iteration.
260 IF (ratio < p0001) GO TO 290
! successful iteration. update x, fvec, and their norms.
DO j = 1, n
x(j) = wa2(j)
wa2(j) = diag(j)*x(j)
END DO
fvec(1:m) = wa4(1:m)
xnorm = enorm(n, wa2)
fnorm = fnorm1
iter = iter + 1
! tests for convergence.
290 IF (ABS(actred) <= ftol .AND. prered <= ftol .AND. p5*ratio <= one) info = 1
IF (delta <= xtol*xnorm) info = 2
IF (ABS(actred) <= ftol .AND. prered <= ftol &
.AND. p5*ratio <= one .AND. info == 2) info = 3
IF (info /= 0) GO TO 300
! tests for termination and stringent tolerances.
IF (nfev >= maxfev) info = 5
IF (ABS(actred) <= epsmch .AND. prered <= epsmch &
.AND. p5*ratio <= one) info = 6
IF (delta <= epsmch*xnorm) info = 7
IF (gnorm <= epsmch) info = 8
IF (info /= 0) GO TO 300
! end of the inner loop. repeat if iteration unsuccessful.
IF (ratio < p0001) GO TO 200
! end of the outer loop.
GO TO 30
! termination, either normal or user imposed.
300 IF (iflag < 0) info = iflag
iflag = 0
IF (nprint > 0) CALL fcn(m, n, x, fvec, iflag)
RETURN
! last card of subroutine lmdif.
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
! Date: 1999-12-09 Time: 12:46:12
! N.B. Arguments LDR, WA1 & WA2 have been removed.
INTEGER, PARAMETER :: dp = dp_kind
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: r(:,:)
INTEGER, INTENT(IN) :: ipvt(:)
REAL (dp), INTENT(IN) :: diag(:)
REAL (dp), INTENT(IN) :: qtb(:)
REAL (dp), INTENT(IN) :: delta
REAL (dp), INTENT(OUT) :: par
REAL (dp), INTENT(OUT) :: x(:)
REAL (dp), INTENT(OUT) :: sdiag(:)
! **********
! subroutine lmpar
! Given an m by n matrix a, an n by n nonsingular diagonal matrix d,
! an m-vector b, and a positive number delta, the problem is to determine a
! value for the parameter par such that if x solves the system
! a*x = b , sqrt(par)*d*x = 0 ,
! in the least squares sense, and dxnorm is the euclidean
! norm of d*x, then either par is zero and
! (dxnorm-delta) <= 0.1*delta ,
! or par is positive and
! abs(dxnorm-delta) <= 0.1*delta .
! This subroutine completes the solution of the problem if it is provided
! with the necessary information from the r factorization, with column
! qpivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
! q has orthogonal columns, and r is an upper triangular matrix with diagonal
! elements of nonincreasing magnitude, then lmpar expects the full upper
! triangle of r, the permutation matrix p, and the first n components of
! (q transpose)*b.
! On output lmpar also provides an upper triangular matrix s such that
! t t t
! p *(a *a + par*d*d)*p = s *s .
! s is employed within lmpar and may be of separate interest.
! Only a few iterations are generally needed for convergence of the algorithm.
! If, however, the limit of 10 iterations is reached, then the output par
! will contain the best value obtained so far.
! the subroutine statement is
! subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag, wa1,wa2)
! where
! n is a positive integer input variable set to the order of r.
! r is an n by n array. on input the full upper triangle
! must contain the full upper triangle of the matrix r.
! On output the full upper triangle is unaltered, and the
! strict lower triangle contains the strict upper triangle
! (transposed) of the upper triangular matrix s.
! ldr is a positive integer input variable not less than n
! which specifies the leading dimension of the array r.
! ipvt is an integer input array of length n which defines the
! permutation matrix p such that a*p = q*r. column j of p
! is column ipvt(j) of the identity matrix.
! diag is an input array of length n which must contain the
! diagonal elements of the matrix d.
! qtb is an input array of length n which must contain the first
! n elements of the vector (q transpose)*b.
! delta is a positive input variable which specifies an upper
! bound on the euclidean norm of d*x.
! par is a nonnegative variable. on input par contains an
! initial estimate of the levenberg-marquardt parameter.
! on output par contains the final estimate.
! x is an output array of length n which contains the least
! squares solution of the system a*x = b, sqrt(par)*d*x = 0,
! for the output par.
! sdiag is an output array of length n which contains the
! diagonal elements of the upper triangular matrix s.
! wa1 and wa2 are work arrays of length n.
! subprograms called
! minpack-supplied ... dpmpar,enorm,qrsolv
! fortran-supplied ... ABS,MAX,MIN,SQRT
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
! **********
INTEGER :: iter, j, jm1, jp1, k, l, nsing
REAL (dp) :: dxnorm, dwarf, fp, gnorm, parc, parl, paru, sum, temp
REAL (dp) :: wa1(n), wa2(n)
REAL (dp), PARAMETER :: p1 = 0.1_dp, p001 = 0.001_dp, zero = 0.0_dp
! dwarf is the smallest positive magnitude.
dwarf = TINY(zero)
! compute and store in x the gauss-newton direction. if the
! jacobian is rank-deficient, obtain a least squares solution.
nsing = n
DO j = 1, n
wa1(j) = qtb(j)
IF (r(j,j) == zero .AND. nsing == n) nsing = j - 1
IF (nsing < n) wa1(j) = zero
END DO
DO k = 1, nsing
j = nsing - k + 1
wa1(j) = wa1(j)/r(j,j)
temp = wa1(j)
jm1 = j - 1
wa1(1:jm1) = wa1(1:jm1) - r(1:jm1,j)*temp
END DO
DO j = 1, n
l = ipvt(j)
x(l) = wa1(j)
END DO
! initialize the iteration counter.
! evaluate the function at the origin, and test
! for acceptance of the gauss-newton direction.
iter = 0
wa2(1:n) = diag(1:n)*x(1:n)
dxnorm = enorm(n, wa2)
fp = dxnorm - delta
IF (fp <= p1*delta) GO TO 220
! if the jacobian is not rank deficient, the newton
! step provides a lower bound, parl, for the zero of
! the function. Otherwise set this bound to zero.
parl = zero
IF (nsing < n) GO TO 120
DO j = 1, n
l = ipvt(j)
wa1(j) = diag(l)*(wa2(l)/dxnorm)
END DO
DO j = 1, n
sum = DOT_PRODUCT( r(1:j-1,j), wa1(1:j-1) )
wa1(j) = (wa1(j) - sum)/r(j,j)
END DO
temp = enorm(n,wa1)
parl = ((fp/delta)/temp)/temp
! calculate an upper bound, paru, for the zero of the function.
120 DO j = 1, n
sum = DOT_PRODUCT( r(1:j,j), qtb(1:j) )
l = ipvt(j)
wa1(j) = sum/diag(l)
END DO
gnorm = enorm(n,wa1)
paru = gnorm/delta
IF (paru == zero) paru = dwarf/MIN(delta,p1)
! if the input par lies outside of the interval (parl,paru),
! set par to the closer endpoint.
par = MAX(par,parl)
par = MIN(par,paru)
IF (par == zero) par = gnorm/dxnorm
! beginning of an iteration.
150 iter = iter + 1
! evaluate the function at the current value of par.
IF (par == zero) par = MAX(dwarf, p001*paru)
temp = SQRT(par)
wa1(1:n) = temp*diag(1:n)
CALL qrsolv(n, r, ipvt, wa1, qtb, x, sdiag)
wa2(1:n) = diag(1:n)*x(1:n)
dxnorm = enorm(n, wa2)
temp = fp
fp = dxnorm - delta
! if the function is small enough, accept the current value
! of par. also test for the exceptional cases where parl
! is zero or the number of iterations has reached 10.
IF (ABS(fp) <= p1*delta .OR. parl == zero .AND. fp <= temp &
.AND. temp < zero .OR. iter == 10) GO TO 220
! compute the newton correction.
DO j = 1, n
l = ipvt(j)
wa1(j) = diag(l)*(wa2(l)/dxnorm)
END DO
DO j = 1, n
wa1(j) = wa1(j)/sdiag(j)
temp = wa1(j)
jp1 = j + 1
wa1(jp1:n) = wa1(jp1:n) - r(jp1:n,j)*temp
END DO
temp = enorm(n,wa1)
parc = ((fp/delta)/temp)/temp
! depending on the sign of the function, update parl or paru.
IF (fp > zero) parl = MAX(parl,par)
IF (fp < zero) paru = MIN(paru,par)
! compute an improved estimate for par.
par = MAX(parl, par+parc)
! end of an iteration.
GO TO 150
! termination.
220 IF (iter == 0) par = zero
RETURN
! last card of subroutine lmpar.
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 = dp_kind
INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: a(:,:)
LOGICAL, INTENT(IN) :: pivot
INTEGER, INTENT(OUT) :: ipvt(:)
REAL (dp), INTENT(OUT) :: rdiag(:)
REAL (dp), INTENT(OUT) :: acnorm(:)
! **********
! subroutine qrfac
! This subroutine uses Householder transformations with column pivoting
! (optional) to compute a qr factorization of the m by n matrix a.
! That is, qrfac determines an orthogonal matrix q, a permutation matrix p,
! and an upper trapezoidal matrix r with diagonal elements of nonincreasing
! magnitude, such that a*p = q*r. The householder transformation for
! column k, k = 1,2,...,min(m,n), is of the form
! t
! i - (1/u(k))*u*u
! where u has zeros in the first k-1 positions. The form of this
! transformation and the method of pivoting first appeared in the
! corresponding linpack subroutine.
! the subroutine statement is
! subroutine qrfac(m, n, a, lda, pivot, ipvt, lipvt, rdiag, acnorm, wa)
! N.B. 3 of these arguments have been omitted in this version.
! where
! m is a positive integer input variable set to the number of rows of a.
! n is a positive integer input variable set to the number of columns of a.
! a is an m by n array. On input a contains the matrix for
! which the qr factorization is to be computed. On output
! the strict upper trapezoidal part of a contains the strict
! upper trapezoidal part of r, and the lower trapezoidal
! part of a contains a factored form of q (the non-trivial
! elements of the u vectors described above).
! lda is a positive integer input variable not less than m
! which specifies the leading dimension of the array a.
! pivot is a logical input variable. If pivot is set true,
! then column pivoting is enforced. If pivot is set false,
! then no column pivoting is done.
! ipvt is an integer output array of length lipvt. ipvt
! defines the permutation matrix p such that a*p = q*r.
! Column j of p is column ipvt(j) of the identity matrix.
! If pivot is false, ipvt is not referenced.
! lipvt is a positive integer input variable. If pivot is false,
! then lipvt may be as small as 1. If pivot is true, then
! lipvt must be at least n.
! rdiag is an output array of length n which contains the
! diagonal elements of r.
! acnorm is an output array of length n which contains the norms of the
! corresponding columns of the input matrix a.
! If this information is not needed, then acnorm can coincide with rdiag.
! wa is a work array of length n. If pivot is false, then wa
! can coincide with rdiag.
! subprograms called
! minpack-supplied ... dpmpar,enorm
! fortran-supplied ... MAX,SQRT,MIN
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
! **********
INTEGER :: i, j, jp1, k, kmax, minmn
REAL (dp) :: ajnorm, epsmch, sum, temp, wa(n)
REAL (dp), PARAMETER :: one = 1.0_dp, p05 = 0.05_dp, zero = 0.0_dp
! epsmch is the machine precision.
epsmch = EPSILON(zero)
! compute the initial column norms and initialize several arrays.
DO j = 1, n
acnorm(j) = enorm(m,a(1:,j))
rdiag(j) = acnorm(j)
wa(j) = rdiag(j)
IF (pivot) ipvt(j) = j
END DO
! Reduce a to r with Householder transformations.
minmn = MIN(m,n)
DO j = 1, minmn
IF (.NOT.pivot) GO TO 40
! Bring the column of largest norm into the pivot position.
kmax = j
DO k = j, n
IF (rdiag(k) > rdiag(kmax)) kmax = k
END DO
IF (kmax == j) GO TO 40
DO i = 1, m
temp = a(i,j)
a(i,j) = a(i,kmax)
a(i,kmax) = temp
END DO
rdiag(kmax) = rdiag(j)
wa(kmax) = wa(j)
k = ipvt(j)
ipvt(j) = ipvt(kmax)
ipvt(kmax) = k
! Compute the Householder transformation to reduce the
! j-th column of a to a multiple of the j-th unit vector.
40 ajnorm = enorm(m-j+1, a(j:,j))
IF (ajnorm == zero) CYCLE
IF (a(j,j) < zero) ajnorm = -ajnorm
a(j:m,j) = a(j:m,j)/ajnorm
a(j,j) = a(j,j) + one
! Apply the transformation to the remaining columns and update the norms.
jp1 = j + 1
DO k = jp1, n
sum = DOT_PRODUCT( a(j:m,j), a(j:m,k) )
temp = sum/a(j,j)
a(j:m,k) = a(j:m,k) - temp*a(j:m,j)
IF (.NOT.pivot .OR. rdiag(k) == zero) CYCLE
temp = a(j,k)/rdiag(k)
rdiag(k) = rdiag(k)*SQRT(MAX(zero, one-temp**2))
IF (p05*(rdiag(k)/wa(k))**2 > epsmch) CYCLE
rdiag(k) = enorm(m-j, a(jp1:,k))
wa(k) = rdiag(k)
END DO
rdiag(j) = -ajnorm
END DO
RETURN
! last card of subroutine qrfac.
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 = dp_kind
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: r(:,:)
INTEGER, INTENT(IN) :: ipvt(:)
REAL (dp), INTENT(IN) :: diag(:)
REAL (dp), INTENT(IN) :: qtb(:)
REAL (dp), INTENT(OUT) :: x(:)
REAL (dp), INTENT(OUT) :: sdiag(:)
! **********
! subroutine qrsolv
! Given an m by n matrix a, an n by n diagonal matrix d, and an m-vector b,
! the problem is to determine an x which solves the system
! a*x = b , d*x = 0 ,
! in the least squares sense.
! This subroutine completes the solution of the problem if it is provided
! with the necessary information from the qr factorization, with column
! pivoting, of a. That is, if a*p = q*r, where p is a permutation matrix,
! q has orthogonal columns, and r is an upper triangular matrix with diagonal
! elements of nonincreasing magnitude, then qrsolv expects the full upper
! triangle of r, the permutation matrix p, and the first n components of
! (q transpose)*b. The system a*x = b, d*x = 0, is then equivalent to
! t t
! r*z = q *b , p *d*p*z = 0 ,
! where x = p*z. if this system does not have full rank,
! then a least squares solution is obtained. On output qrsolv
! also provides an upper triangular matrix s such that
! t t t
! p *(a *a + d*d)*p = s *s .
! s is computed within qrsolv and may be of separate interest.
! the subroutine statement is
! subroutine qrsolv(n, r, ldr, ipvt, diag, qtb, x, sdiag, wa)
! N.B. Arguments LDR and WA have been removed in this version.
! where
! n is a positive integer input variable set to the order of r.
! r is an n by n array. On input the full upper triangle must contain
! the full upper triangle of the matrix r.
! On output the full upper triangle is unaltered, and the strict lower
! triangle contains the strict upper triangle (transposed) of the
! upper triangular matrix s.
! ldr is a positive integer input variable not less than n
! which specifies the leading dimension of the array r.
! ipvt is an integer input array of length n which defines the
! permutation matrix p such that a*p = q*r. Column j of p
! is column ipvt(j) of the identity matrix.
! diag is an input array of length n which must contain the
! diagonal elements of the matrix d.
! qtb is an input array of length n which must contain the first
! n elements of the vector (q transpose)*b.
! x is an output array of length n which contains the least
! squares solution of the system a*x = b, d*x = 0.
! sdiag is an output array of length n which contains the
! diagonal elements of the upper triangular matrix s.
! wa is a work array of length n.
! subprograms called
! fortran-supplied ... ABS,SQRT
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
! **********
INTEGER :: i, j, k, kp1, l, nsing
REAL (dp) :: COS, cotan, qtbpj, SIN, sum, TAN, temp, wa(n)
REAL (dp), PARAMETER :: p5 = 0.5_dp, p25 = 0.25_dp, zero = 0.0_dp
! Copy r and (q transpose)*b to preserve input and initialize s.
! In particular, save the diagonal elements of r in x.
DO j = 1, n
r(j:n,j) = r(j,j:n)
x(j) = r(j,j)
wa(j) = qtb(j)
END DO
! Eliminate the diagonal matrix d using a givens rotation.
DO j = 1, n
! Prepare the row of d to be eliminated, locating the
! diagonal element using p from the qr factorization.
l = ipvt(j)
IF (diag(l) == zero) CYCLE
sdiag(j:n) = zero
sdiag(j) = diag(l)
! The transformations to eliminate the row of d modify only a single
! element of (q transpose)*b beyond the first n, which is initially zero.
qtbpj = zero
DO k = j, n
! Determine a givens rotation which eliminates the
! appropriate element in the current row of d.
IF (sdiag(k) == zero) CYCLE
IF (ABS(r(k,k)) < ABS(sdiag(k))) THEN
cotan = r(k,k)/sdiag(k)
SIN = p5/SQRT(p25 + p25*cotan**2)
COS = SIN*cotan
ELSE
TAN = sdiag(k)/r(k,k)
COS = p5/SQRT(p25 + p25*TAN**2)
SIN = COS*TAN
END IF
! Compute the modified diagonal element of r and
! the modified element of ((q transpose)*b,0).
r(k,k) = COS*r(k,k) + SIN*sdiag(k)
temp = COS*wa(k) + SIN*qtbpj
qtbpj = -SIN*wa(k) + COS*qtbpj
wa(k) = temp
! Accumulate the tranformation in the row of s.
kp1 = k + 1
DO i = kp1, n
temp = COS*r(i,k) + SIN*sdiag(i)
sdiag(i) = -SIN*r(i,k) + COS*sdiag(i)
r(i,k) = temp
END DO
END DO
! Store the diagonal element of s and restore
! the corresponding diagonal element of r.
sdiag(j) = r(j,j)
r(j,j) = x(j)
END DO
! Solve the triangular system for z. If the system is singular,
! then obtain a least squares solution.
nsing = n
DO j = 1, n
IF (sdiag(j) == zero .AND. nsing == n) nsing = j - 1
IF (nsing < n) wa(j) = zero
END DO
DO k = 1, nsing
j = nsing - k + 1
sum = DOT_PRODUCT( r(j+1:nsing,j), wa(j+1:nsing) )
wa(j) = (wa(j) - sum)/sdiag(j)
END DO
! Permute the components of z back to components of x.
DO j = 1, n
l = ipvt(j)
x(l) = wa(j)
END DO
RETURN
! last card of subroutine qrsolv.
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 = dp_kind
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: x(:)
REAL (dp) :: fn_val
! **********
! function enorm
! given an n-vector x, this function calculates the euclidean norm of x.
! the euclidean norm is computed by accumulating the sum of squares in
! three different sums. The sums of squares for the small and large
! components are scaled so that no overflows occur. Non-destructive
! underflows are permitted. Underflows and overflows do not occur in the
! computation of the unscaled sum of squares for the intermediate
! components. The definitions of small, intermediate and large components
! depend on two constants, rdwarf and rgiant. The main restrictions on
! these constants are that rdwarf**2 not underflow and rgiant**2 not
! overflow. The constants given here are suitable for every known computer.
! the function statement is
! REAL (dp) function enorm(n,x)
! where
! n is a positive integer input variable.
! x is an input array of length n.
! subprograms called
! fortran-supplied ... ABS,SQRT
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
! **********
INTEGER :: i
REAL (dp) :: agiant, floatn, s1, s2, s3, xabs, x1max, x3max
REAL (dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp, rdwarf = 3.834E-20_dp, &
rgiant = 1.304E+19_dp
s1 = zero
s2 = zero
s3 = zero
x1max = zero
x3max = zero
floatn = n
agiant = rgiant/floatn
DO i = 1, n
xabs = ABS(x(i))
IF (xabs > rdwarf .AND. xabs < agiant) GO TO 70
IF (xabs <= rdwarf) GO TO 30
! sum for large components.
IF (xabs <= x1max) GO TO 10
s1 = one + s1*(x1max/xabs)**2
x1max = xabs
GO TO 20
10 s1 = s1 + (xabs/x1max)**2
20 GO TO 60
! sum for small components.
30 IF (xabs <= x3max) GO TO 40
s3 = one + s3*(x3max/xabs)**2
x3max = xabs
GO TO 60
40 IF (xabs /= zero) s3 = s3 + (xabs/x3max)**2
60 CYCLE
! sum for intermediate components.
70 s2 = s2 + xabs**2
END DO
! calculation of norm.
IF (s1 == zero) GO TO 100
fn_val = x1max*SQRT(s1 + (s2/x1max)/x1max)
GO TO 120
100 IF (s2 == zero) GO TO 110
IF (s2 >= x3max) fn_val = SQRT(s2*(one + (x3max/s2)*(x3max*s3)))
IF (s2 < x3max) fn_val = SQRT(x3max*((s2/x3max) + (x3max*s3)))
GO TO 120
110 fn_val = x3max*SQRT(s3)
120 RETURN
! last card of function enorm.
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 = dp_kind
INTEGER, INTENT(IN) :: m
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN OUT) :: x(n)
REAL (dp), INTENT(IN) :: fvec(m)
REAL (dp), INTENT(OUT) :: fjac(:,:) ! fjac(ldfjac,n)
INTEGER, INTENT(IN OUT) :: iflag
REAL (dp), INTENT(IN) :: epsfcn
INTERFACE
SUBROUTINE fcn(m, n, x, fvec, 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
! **********
! subroutine fdjac2
! this subroutine computes a forward-difference approximation
! to the m by n jacobian matrix associated with a specified
! problem of m functions in n variables.
! the subroutine statement is
! subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
! 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 fdjac2.
! 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 input array of length n.
! fvec is an input array of length m which must contain the
! functions evaluated at x.
! fjac is an output m by n array which contains the
! approximation to the jacobian matrix evaluated at x.
! ldfjac is a positive integer input variable not less than m
! which specifies the leading dimension of the array fjac.
! iflag is an integer variable which can be used to terminate
! the execution of fdjac2. see description of fcn.
! epsfcn is 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.
! wa is a work array of length m.
! subprograms called
! user-supplied ...... fcn
! minpack-supplied ... dpmpar
! fortran-supplied ... ABS,MAX,SQRT
! argonne national laboratory. minpack project. march 1980.
! burton s. garbow, kenneth e. hillstrom, jorge j. more
! **********
INTEGER :: j
REAL (dp) :: eps, epsmch, h, temp, wa(m)
REAL (dp), PARAMETER :: zero = 0.0_dp
! epsmch is the machine precision.
epsmch = EPSILON(zero)
eps = SQRT(MAX(epsfcn, epsmch))
DO j = 1, n
temp = x(j)
h = eps*ABS(temp)
IF (h == zero) h = eps
x(j) = temp + h
CALL fcn(m, n, x, wa, iflag)
IF (iflag < 0) EXIT
x(j) = temp
fjac(1:m,j) = (wa(1:m) - fvec(1:m))/h
END DO
RETURN
! last card of subroutine fdjac2.
END SUBROUTINE fdjac2
end module fvn_misc