##// END OF EJS Templates
schain
RSS
Faraday manual processing improved
Faraday manual processing improved
rflores - - Load All Authors

File last commit:

r1601:3b2826fe0d97
r1774:e7bb496918b2
Show More
lmdif1.f
1567 lines | 48.2 KiB | text/x-fortran | FortranFixedLexer
/ schainf / Ffiles / lmdif1.f
History | Source | Raw |Copy content |Copy permalink
rflores
add Ffiles
 r1601
subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
integer m,n,info,lwa
integer iwa(n)
real tol
real x(n),fvec(m),wa(lwa)
external fcn
c **********
c
c subroutine lmdif1
c
c the purpose of lmdif1 is to minimize the sum of the squares of
c m nonlinear functions in n variables by a modification of the
c levenberg-marquardt algorithm. this is done by using the more
c general least-squares solver lmdif. the user must provide a
c subroutine which calculates the functions. the jacobian is
c then calculated by a forward-difference approximation.
c
c the subroutine statement is
c
c subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa)
c
c where
c
c fcn is the name of the user-supplied subroutine which
c calculates the functions. fcn must be declared
c in an external statement in the user calling
c program, and should be written as follows.
c
c subroutine fcn(m,n,x,fvec,iflag
c integer m,n,iflag
c real x(n),fvec(m)
c ----------
c calculate the functions at x and
c return this vector in fvec.
c ----------
c return
c end
c
c the value of iflag should not be changed by fcn unless
c the user wants to terminate execution of lmdif1.
c in this case set iflag to a negative integer.
c
c m is a positive integer input variable set to the number
c of functions.
c
c n is a positive integer input variable set to the number
c of variables. n must not exceed m.
c
c x is an array of length n. on input x must contain
c an initial estimate of the solution vector. on output x
c contains the final estimate of the solution vector.
c
c fvec is an output array of length m which contains
c the functions evaluated at the output x.
c
c tol is a nonnegative input variable. termination occurs
c when the algorithm estimates either that the relative
c error in the sum of squares is at most tol or that
c the relative error between x and the solution is at
c most tol.
c
c info is an integer output variable. if the user has
c terminated execution, info is set to the (negative)
c value of iflag. see description of fcn. otherwise,
c info is set as follows.
c
c info = 0 improper input parameters.
c
c info = 1 algorithm estimates that the relative error
c in the sum of squares is at most tol.
c
c info = 2 algorithm estimates that the relative error
c between x and the solution is at most tol.
c
c info = 3 conditions for info = 1 and info = 2 both hold.
c
c info = 4 fvec is orthogonal to the columns of the
c jacobian to machine precision.
c
c info = 5 number of calls to fcn has reached or
c exceeded 200*(n+1).
c
c info = 6 tol is too small. no further reduction in
c the sum of squares is possible.
c
c info = 7 tol is too small. no further improvement in
c the approximate solution x is possible.
c
c iwa is an integer work array of length n.
c
c wa is a work array of length lwa.
c
c lwa is a positive integer input variable not less than
c m*n+5*n+m.
c
c subprograms called
c
c user-supplied ...... fcn
c
c minpack-supplied ... lmdif
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer maxfev,mode,mp5n,nfev,nprint
real epsfcn,factor,ftol,gtol,xtol,zero
data factor,zero /1.0e2,0.0e0/
info = 0
c
c check the input parameters for errors.
c
if (n .le. 0 .or. m .lt. n .or. tol .lt. zero
* .or. lwa .lt. m*n + 5*n + m) go to 10
c
c call lmdif.
c
maxfev = 200*(n + 1)
ftol = tol
xtol = tol
gtol = zero
epsfcn = zero
mode = 1
nprint = 0
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))
if (info .eq. 8) info = 4
10 continue
return
c
c last card of subroutine lmdif1.
c
end
real function enorm(n,x)
integer n
real x(n)
c **********
c
c function enorm
c
c given an n-vector x, this function calculates the
c euclidean norm of x.
c
c the euclidean norm is computed by accumulating the sum of
c squares in three different sums. the sums of squares for the
c small and large components are scaled so that no overflows
c occur. non-destructive underflows are permitted. underflows
c and overflows do not occur in the computation of the unscaled
c sum of squares for the intermediate components.
c the definitions of small, intermediate and large components
c depend on two constants, rdwarf and rgiant. the main
c restrictions on these constants are that rdwarf**2 not
c underflow and rgiant**2 not overflow. the constants
c given here are suitable for every known computer.
c
c the function statement is
c
c real function enorm(n,x)
c
c where
c
c n is a positive integer input variable.
c
c x is an input array of length n.
c
c subprograms called
c
c fortran-supplied ... abs,sqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i
real agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,x1max,x3max,
* zero
data one,zero,rdwarf,rgiant /1.0e0,0.0e0,3.834e-20,1.304e19/
s1 = zero
s2 = zero
s3 = zero
x1max = zero
x3max = zero
floatn = n
agiant = rgiant/floatn
do 90 i = 1, n
xabs = abs(x(i))
if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70
if (xabs .le. rdwarf) go to 30
c
c sum for large components.
c
if (xabs .le. x1max) go to 10
s1 = one + s1*(x1max/xabs)**2
x1max = xabs
go to 20
10 continue
s1 = s1 + (xabs/x1max)**2
20 continue
go to 60
30 continue
c
c sum for small components.
c
if (xabs .le. x3max) go to 40
s3 = one + s3*(x3max/xabs)**2
x3max = xabs
go to 50
40 continue
if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2
50 continue
60 continue
go to 80
70 continue
c
c sum for intermediate components.
c
s2 = s2 + xabs**2
80 continue
90 continue
c
c calculation of norm.
c
if (s1 .eq. zero) go to 100
enorm = x1max*sqrt(s1+(s2/x1max)/x1max)
go to 130
100 continue
if (s2 .eq. zero) go to 110
if (s2 .ge. x3max)
* enorm = sqrt(s2*(one+(x3max/s2)*(x3max*s3)))
if (s2 .lt. x3max)
* enorm = sqrt(x3max*((s2/x3max)+(x3max*s3)))
go to 120
110 continue
enorm = x3max*sqrt(s3)
120 continue
130 continue
return
c
c last card of function enorm.
c
end
subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
integer m,n,ldfjac,iflag
real epsfcn
real x(n),fvec(m),fjac(ldfjac,n),wa(m)
c **********
c
c subroutine fdjac2
c
c this subroutine computes a forward-difference approximation
c to the m by n jacobian matrix associated with a specified
c problem of m functions in n variables.
c
c the subroutine statement is
c
c subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
c
c where
c
c fcn is the name of the user-supplied subroutine which
c calculates the functions. fcn must be declared
c in an external statement in the user calling
c program, and should be written as follows.
c
c subroutine fcn(m,n,x,fvec,iflag)
c integer m,n,iflag
c real x(n),fvec(m)
c ----------
c calculate the functions at x and
c return this vector in fvec.
c ----------
c return
c end
c
c the value of iflag should not be changed by fcn unless
c the user wants to terminate execution of fdjac2.
c in this case set iflag to a negative integer.
c
c m is a positive integer input variable set to the number
c of functions.
c
c n is a positive integer input variable set to the number
c of variables. n must not exceed m.
c
c x is an input array of length n.
c
c fvec is an input array of length m which must contain the
c functions evaluated at x.
c
c fjac is an output m by n array which contains the
c approximation to the jacobian matrix evaluated at x.
c
c ldfjac is a positive integer input variable not less than m
c which specifies the leading dimension of the array fjac.
c
c iflag is an integer variable which can be used to terminate
c the execution of fdjac2. see description of fcn.
c
c epsfcn is an input variable used in determining a suitable
c step length for the forward-difference approximation. this
c approximation assumes that the relative errors in the
c functions are of the order of epsfcn. if epsfcn is less
c than the machine precision, it is assumed that the relative
c errors in the functions are of the order of the machine
c precision.
c
c wa is a work array of length m.
c
c subprograms called
c
c user-supplied ...... fcn
c
c minpack-supplied ... spmpar
c
c fortran-supplied ... abs,amax1,sqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j
real eps,epsmch,h,temp,zero
real spmpar
data zero /0.0e0/
c
c epsmch is the machine precision.
c
epsmch = spmpar(1)
write(*,*) "epsfcn,epsmch",epsfcn,epsmch
c call exit
c
c epsmch=1.18999999e-07
eps = sqrt(amax1(epsfcn,epsmch))
do 20 j = 1, n
temp = x(j)
h = eps*abs(temp)
c write(*,*) h,eps,temp
c call exit
if (h .eq. zero) h = eps
x(j) = temp + h
call fcn(m,n,x,wa,iflag)
if (iflag .lt. 0) go to 30
x(j) = temp
do 10 i = 1, m
fjac(i,j) = (wa(i) - fvec(i))/h
c if (j .eq. n) then
c write(*,*) j,i,wa(i),fvec(i),h,fjac(i,j)
c end if
10 continue
c call exit
c write(*,*) fjac(:,1)
c call exit
20 continue
30 continue
c write(*,*) fjac(m-26,n)
c call exit
return
c
c last card of subroutine fdjac2.
c
end
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)
integer m,n,maxfev,mode,nprint,info,nfev,ldfjac
integer ipvt(n)
real ftol,xtol,gtol,epsfcn,factor,err(20)
real x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n),wa1(n),wa2(n),
* wa3(n),wa4(m)
external fcn
c **********
c
c subroutine lmdif
c
c the purpose of lmdif is to minimize the sum of the squares of
c m nonlinear functions in n variables by a modification of
c the levenberg-marquardt algorithm. the user must provide a
c subroutine which calculates the functions. the jacobian is
c then calculated by a forward-difference approximation.
c
c the subroutine statement is
c
c subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
c diag,mode,factor,nprint,info,nfev,fjac,
c ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
c
c where
c
c fcn is the name of the user-supplied subroutine which
c calculates the functions. fcn must be declared
c in an external statement in the user calling
c program, and should be written as follows.
c
c subroutine fcn(m,n,x,fvec,iflag)
c integer m,n,iflag
c real x(n),fvec(m)
c ----------
c calculate the functions at x and
c return this vector in fvec.
c ----------
c return
c end
c
c the value of iflag should not be changed by fcn unless
c the user wants to terminate execution of lmdif.
c in this case set iflag to a negative integer.
c
c m is a positive integer input variable set to the number
c of functions.
c
c n is a positive integer input variable set to the number
c of variables. n must not exceed m.
c
c x is an array of length n. on input x must contain
c an initial estimate of the solution vector. on output x
c contains the final estimate of the solution vector.
c
c fvec is an output array of length m which contains
c the functions evaluated at the output x.
c
c ftol is a nonnegative input variable. termination
c occurs when both the actual and predicted relative
c reductions in the sum of squares are at most ftol.
c therefore, ftol measures the relative error desired
c in the sum of squares.
c
c xtol is a nonnegative input variable. termination
c occurs when the relative error between two consecutive
c iterates is at most xtol. therefore, xtol measures the
c relative error desired in the approximate solution.
c
c gtol is a nonnegative input variable. termination
c occurs when the cosine of the angle between fvec and
c any column of the jacobian is at most gtol in absolute
c value. therefore, gtol measures the orthogonality
c desired between the function vector and the columns
c of the jacobian.
c
c maxfev is a positive integer input variable. termination
c occurs when the number of calls to fcn is at least
c maxfev by the end of an iteration.
c
c epsfcn is an input variable used in determining a suitable
c step length for the forward-difference approximation. this
c approximation assumes that the relative errors in the
c functions are of the order of epsfcn. if epsfcn is less
c than the machine precision, it is assumed that the relative
c errors in the functions are of the order of the machine
c precision.
c
c diag is an array of length n. if mode = 1 (see
c below), diag is internally set. if mode = 2, diag
c must contain positive entries that serve as
c multiplicative scale factors for the variables.
c
c mode is an integer input variable. if mode = 1, the
c variables will be scaled internally. if mode = 2,
c the scaling is specified by the input diag. other
c values of mode are equivalent to mode = 1.
c
c factor is a positive input variable used in determining the
c initial step bound. this bound is set to the product of
c factor and the euclidean norm of diag*x if nonzero, or else
c to factor itself. in most cases factor should lie in the
c interval (.1,100.). 100. is a generally recommended value.
c
c nprint is an integer input variable that enables controlled
c printing of iterates if it is positive. in this case,
c fcn is called with iflag = 0 at the beginning of the first
c iteration and every nprint iterations thereafter and
c immediately prior to return, with x and fvec available
c for printing. if nprint is not positive, no special calls
c of fcn with iflag = 0 are made.
c
c info is an integer output variable. if the user has
c terminated execution, info is set to the (negative)
c value of iflag. see description of fcn. otherwise,
c info is set as follows.
c
c info = 0 improper input parameters.
c
c info = 1 both actual and predicted relative reductions
c in the sum of squares are at most ftol.
c
c info = 2 relative error between two consecutive iterates
c is at most xtol.
c
c info = 3 conditions for info = 1 and info = 2 both hold.
c
c info = 4 the cosine of the angle between fvec and any
c column of the jacobian is at most gtol in
c absolute value.
c
c info = 5 number of calls to fcn has reached or
c exceeded maxfev.
c
c info = 6 ftol is too small. no further reduction in
c the sum of squares is possible.
c
c info = 7 xtol is too small. no further improvement in
c the approximate solution x is possible.
c
c info = 8 gtol is too small. fvec is orthogonal to the
c columns of the jacobian to machine precision.
c
c nfev is an integer output variable set to the number of
c calls to fcn.
c
c fjac is an output m by n array. the upper n by n submatrix
c of fjac contains an upper triangular matrix r with
c diagonal elements of nonincreasing magnitude such that
c
c t t t
c p *(jac *jac)*p = r *r,
c
c where p is a permutation matrix and jac is the final
c calculated jacobian. column j of p is column ipvt(j)
c (see below) of the identity matrix. the lower trapezoidal
c part of fjac contains information generated during
c the computation of r.
c
c ldfjac is a positive integer input variable not less than m
c which specifies the leading dimension of the array fjac.
c
c ipvt is an integer output array of length n. ipvt
c defines a permutation matrix p such that jac*p = q*r,
c where jac is the final calculated jacobian, q is
c orthogonal (not stored), and r is upper triangular
c with diagonal elements of nonincreasing magnitude.
c column j of p is column ipvt(j) of the identity matrix.
c
c qtf is an output array of length n which contains
c the first n elements of the vector (q transpose)*fvec.
c
c wa1, wa2, and wa3 are work arrays of length n.
c
c wa4 is a work array of length m.
c
c subprograms called
c
c user-supplied ...... fcn
c
c minpack-supplied ... spmpar,enorm,fdjac2,lmpar,qrfac
c
c fortran-supplied ... abs,amax1,amin1,sqrt,mod
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,iflag,iter,j,l
real actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm,one,par,
* pnorm,prered,p1,p5,p25,p75,p0001,ratio,sum,temp,temp1,
* temp2,xnorm,zero
real spmpar,enorm
C common /errors/err
data one,p1,p5,p25,p75,p0001,zero
* /1.0e0,1.0e-1,5.0e-1,2.5e-1,7.5e-1,1.0e-4,0.0e0/
c
c epsmch is the machine precision.
c
epsmch = spmpar(1)
c
info = 0
iflag = 0
nfev = 0
c
c check the input parameters for errors.
c
if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m
* .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero
* .or. maxfev .le. 0 .or. factor .le. zero) go to 300
if (mode .ne. 2) go to 20
do 10 j = 1, n
if (diag(j) .le. zero) go to 300
10 continue
20 continue
c
c evaluate the function at the starting point
c and calculate its norm.
c
iflag = 1
c write(*,*) "before calling fcn"
c write(*,*) "fvec_before: ",fvec(222)
c write(*,*) "x_before: ",x
c write(*,*) fjac(1,1)
c call exit
call fcn(m,n,x,fvec,iflag)
c write(*,*) "x_after: ",x
c write(*,*) "fvec(2): ",fvec(222)
c write(*,*) fjac(1,1)
c call exit
c call exit
c write(*,*) "fvec_after: ",fvec(1)
c write(*,*) "after calling fcn"
c call exit
nfev = 1
if (iflag .lt. 0) go to 300
fnorm = enorm(m,fvec)
c
c initialize levenberg-marquardt parameter and iteration counter.
c
par = zero
iter = 1
c
c beginning of the outer loop.
c
30 continue
c
c calculate the jacobian matrix.
c
c write(*,*) "x",x
c call exit
iflag = 2
write(*,*) fjac(1,1)
call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4)
nfev = nfev + n
write(*,*) "2",fjac(1,1)
c write(*,*) wa4
c call exit
c write(*,*) "fvec(2): ",fvec(2)
c call exit
if (iflag .lt. 0) go to 300
c
c if requested, call fcn to enable printing of iterates.
c
if (nprint .le. 0) go to 40
iflag = 0
if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag)
if (iflag .lt. 0) go to 300
40 continue
c
c compute the qr factorization of the jacobian.
c
c write(*,*) fjac(1,1)
c call exit
call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3)
c write(*,*) fjac(12,13)
c write(*,*) "after qrfac"
c call exit
c
c on the first iteration and if mode is 1, scale according
c to the norms of the columns of the initial jacobian.
c
if (iter .ne. 1) go to 80
if (mode .eq. 2) go to 60
do 50 j = 1, n
diag(j) = wa2(j)
if (wa2(j) .eq. zero) diag(j) = one
50 continue
60 continue
c
c on the first iteration, calculate the norm of the scaled x
c and initialize the step bound delta.
c
do 70 j = 1, n
wa3(j) = diag(j)*x(j)
70 continue
xnorm = enorm(n,wa3)
delta = factor*xnorm
if (delta .eq. zero) delta = factor
80 continue
c
c form (q transpose)*fvec and store the first n components in
c qtf.
c
c write(*,*) fjac(12,13)
c call exit
do 90 i = 1, m
wa4(i) = fvec(i)
90 continue
c write(*,*) wa4(222)
c call exit
do 130 j = 1, n
if (fjac(j,j) .eq. zero) go to 120
sum = zero
do 100 i = j, m
c write(*,*) fjac(i,j),wa4(i)
c call exit
sum = sum + fjac(i,j)*wa4(i)
100 continue
temp = -sum/fjac(j,j)
c write(*,*) sum
c call exit
do 110 i = j, m
wa4(i) = wa4(i) + fjac(i,j)*temp
110 continue
120 continue
fjac(j,j) = wa1(j)
qtf(j) = wa4(j)
130 continue
c write(*,*) wa4(222)
c call exit
c
c compute the norm of the scaled gradient.
c
gnorm = zero
if (fnorm .eq. zero) go to 170
do 160 j = 1, n
l = ipvt(j)
if (wa2(l) .eq. zero) go to 150
sum = zero
do 140 i = 1, j
sum = sum + fjac(i,j)*(qtf(i)/fnorm)
140 continue
gnorm = amax1(gnorm,abs(sum/wa2(l)))
150 continue
160 continue
170 continue
c
c test for convergence of the gradient norm.
c
if (gnorm .le. gtol) info = 4
if (info .ne. 0) go to 300
c
c rescale if necessary.
c
if (mode .eq. 2) go to 190
do 180 j = 1, n
diag(j) = amax1(diag(j),wa2(j))
180 continue
190 continue
c
c beginning of the inner loop.
c
200 continue
c
c determine the levenberg-marquardt parameter.
c
write(*,*) "Starting lmpar"
write(*,*) "x(2)",x(2),"wa4(222)",wa4(222)
call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2,
* wa3,wa4)
c
c store the direction p and x + p. calculate the norm of p.
c
do 210 j = 1, n
wa1(j) = -wa1(j)
wa2(j) = x(j) + wa1(j)
wa3(j) = diag(j)*wa1(j)
210 continue
pnorm = enorm(n,wa3)
c
c on the first iteration, adjust the initial step bound.
c
if (iter .eq. 1) delta = amin1(delta,pnorm)
c
c evaluate the function at x + p and calculate its norm.
c
iflag = 1
write(*,*) "starting last-1 fcn"
write(*,*) "wa4_before",wa4(222)
call fcn(m,n,wa2,wa4,iflag)
nfev = nfev + 1
c call exit
if (iflag .lt. 0) go to 300
fnorm1 = enorm(m,wa4)
write(*,*) "wa4",wa4(222)
c call exit
c
c compute the scaled actual reduction.
c
write(*,*) "fnorm",fnorm,"fnorm1",fnorm1
actred = -one
if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
c
c compute the scaled predicted reduction and
c the scaled directional derivative.
c
do 230 j = 1, n
wa3(j) = zero
l = ipvt(j)
temp = wa1(l)
do 220 i = 1, j
wa3(i) = wa3(i) + fjac(i,j)*temp
220 continue
230 continue
temp1 = enorm(n,wa3)/fnorm
temp2 = (sqrt(par)*pnorm)/fnorm
prered = temp1**2 + temp2**2/p5
dirder = -(temp1**2 + temp2**2)
c
c compute the ratio of the actual to the predicted
c reduction.
c
write(*,*) "ratio",ratio
ratio = zero
write(*,*) "ratio2",ratio
c call exit
if (prered .ne. zero) ratio = actred/prered
write(*,*) "ratio3",ratio,actred,prered
c
c update the step bound.
c
if (ratio .gt. p25) go to 240
if (actred .ge. zero) temp = p5
if (actred .lt. zero)
* temp = p5*dirder/(dirder + p5*actred)
if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1
delta = temp*amin1(delta,pnorm/p1)
par = par/temp
go to 260
240 continue
if (par .ne. zero .and. ratio .lt. p75) go to 250
delta = pnorm/p5
par = p5*par
250 continue
260 continue
c
c test for successful iteration.
c
if (ratio .lt. p0001) go to 290
c
c successful iteration. update x, fvec, and their norms.
c
do 270 j = 1, n
C err(j)=abs(wa2(j)-x(j))
x(j) = wa2(j)
wa2(j) = diag(j)*x(j)
270 continue
do 280 i = 1, m
fvec(i) = wa4(i)
280 continue
xnorm = enorm(n,wa2)
fnorm = fnorm1
iter = iter + 1
290 continue
c
c tests for convergence.
c
if (abs(actred) .le. ftol .and. prered .le. ftol
* .and. p5*ratio .le. one) info = 1
if (delta .le. xtol*xnorm) info = 2
if (abs(actred) .le. ftol .and. prered .le. ftol
* .and. p5*ratio .le. one .and. info .eq. 2) info = 3
write(*,*) "info",info
if (info .ne. 0) go to 300
c
c tests for termination and stringent tolerances.
c
write(*,*) delta,epsmch,xnorm
if (nfev .ge. maxfev) info = 5
if (abs(actred) .le. epsmch .and. prered .le. epsmch
* .and. p5*ratio .le. one) info = 6
if (delta .le. epsmch*xnorm) info = 7
if (gnorm .le. epsmch) info = 8
write(*,*) "info2",info
if (info .ne. 0) go to 300
c
c end of the inner loop. repeat if iteration unsuccessful.
c
if (ratio .lt. p0001) go to 200
c
c end of the outer loop.
c
go to 30
300 continue
c
c termination, either normal or user imposed.
c
write(*,*) "starting last fcn"
if (iflag .lt. 0) info = iflag
iflag = 0
if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag)
c call exit
return
c
c last card of subroutine lmdif.
c
end
subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1,
* wa2)
integer n,ldr
integer ipvt(n)
real delta,par
real r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa1(n),wa2(n)
c **********
c
c subroutine lmpar
c
c given an m by n matrix a, an n by n nonsingular diagonal
c matrix d, an m-vector b, and a positive number delta,
c the problem is to determine a value for the parameter
c par such that if x solves the system
c
c a*x = b , sqrt(par)*d*x = 0 ,
c
c in the least squares sense, and dxnorm is the euclidean
c norm of d*x, then either par is zero and
c
c (dxnorm-delta) .le. 0.1*delta ,
c
c or par is positive and
c
c abs(dxnorm-delta) .le. 0.1*delta .
c
c this subroutine completes the solution of the problem
c if it is provided with the necessary information from the
c qr factorization, with column pivoting, of a. that is, if
c a*p = q*r, where p is a permutation matrix, q has orthogonal
c columns, and r is an upper triangular matrix with diagonal
c elements of nonincreasing magnitude, then lmpar expects
c the full upper triangle of r, the permutation matrix p,
c and the first n components of (q transpose)*b. on output
c lmpar also provides an upper triangular matrix s such that
c
c t t t
c p *(a *a + par*d*d)*p = s *s .
c
c s is employed within lmpar and may be of separate interest.
c
c only a few iterations are generally needed for convergence
c of the algorithm. if, however, the limit of 10 iterations
c is reached, then the output par will contain the best
c value obtained so far.
c
c the subroutine statement is
c
c subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,
c wa1,wa2)
c
c where
c
c n is a positive integer input variable set to the order of r.
c
c r is an n by n array. on input the full upper triangle
c must contain the full upper triangle of the matrix r.
c on output the full upper triangle is unaltered, and the
c strict lower triangle contains the strict upper triangle
c (transposed) of the upper triangular matrix s.
c
c ldr is a positive integer input variable not less than n
c which specifies the leading dimension of the array r.
c
c ipvt is an integer input array of length n which defines the
c permutation matrix p such that a*p = q*r. column j of p
c is column ipvt(j) of the identity matrix.
c
c diag is an input array of length n which must contain the
c diagonal elements of the matrix d.
c
c qtb is an input array of length n which must contain the first
c n elements of the vector (q transpose)*b.
c
c delta is a positive input variable which specifies an upper
c bound on the euclidean norm of d*x.
c
c par is a nonnegative variable. on input par contains an
c initial estimate of the levenberg-marquardt parameter.
c on output par contains the final estimate.
c
c x is an output array of length n which contains the least
c squares solution of the system a*x = b, sqrt(par)*d*x = 0,
c for the output par.
c
c sdiag is an output array of length n which contains the
c diagonal elements of the upper triangular matrix s.
c
c wa1 and wa2 are work arrays of length n.
c
c subprograms called
c
c minpack-supplied ... spmpar,enorm,qrsolv
c
c fortran-supplied ... abs,amax1,amin1,sqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,iter,j,jm1,jp1,k,l,nsing
real dxnorm,dwarf,fp,gnorm,parc,parl,paru,p1,p001,sum,temp,zero
real spmpar,enorm
data p1,p001,zero /1.0e-1,1.0e-3,0.0e0/
c
c dwarf is the smallest positive magnitude.
c
dwarf = spmpar(2)
c
c compute and store in x the gauss-newton direction. if the
c jacobian is rank-deficient, obtain a least squares solution.
c
nsing = n
do 10 j = 1, n
wa1(j) = qtb(j)
if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1
if (nsing .lt. n) wa1(j) = zero
10 continue
if (nsing .lt. 1) go to 50
do 40 k = 1, nsing
j = nsing - k + 1
wa1(j) = wa1(j)/r(j,j)
temp = wa1(j)
jm1 = j - 1
if (jm1 .lt. 1) go to 30
do 20 i = 1, jm1
wa1(i) = wa1(i) - r(i,j)*temp
20 continue
30 continue
40 continue
50 continue
do 60 j = 1, n
l = ipvt(j)
x(l) = wa1(j)
60 continue
c
c initialize the iteration counter.
c evaluate the function at the origin, and test
c for acceptance of the gauss-newton direction.
c
iter = 0
do 70 j = 1, n
wa2(j) = diag(j)*x(j)
70 continue
dxnorm = enorm(n,wa2)
fp = dxnorm - delta
if (fp .le. p1*delta) go to 220
c
c if the jacobian is not rank deficient, the newton
c step provides a lower bound, parl, for the zero of
c the function. otherwise set this bound to zero.
c
parl = zero
if (nsing .lt. n) go to 120
do 80 j = 1, n
l = ipvt(j)
wa1(j) = diag(l)*(wa2(l)/dxnorm)
80 continue
do 110 j = 1, n
sum = zero
jm1 = j - 1
if (jm1 .lt. 1) go to 100
do 90 i = 1, jm1
sum = sum + r(i,j)*wa1(i)
90 continue
100 continue
wa1(j) = (wa1(j) - sum)/r(j,j)
110 continue
temp = enorm(n,wa1)
parl = ((fp/delta)/temp)/temp
120 continue
c
c calculate an upper bound, paru, for the zero of the function.
c
do 140 j = 1, n
sum = zero
do 130 i = 1, j
sum = sum + r(i,j)*qtb(i)
130 continue
l = ipvt(j)
wa1(j) = sum/diag(l)
140 continue
gnorm = enorm(n,wa1)
paru = gnorm/delta
if (paru .eq. zero) paru = dwarf/amin1(delta,p1)
c
c if the input par lies outside of the interval (parl,paru),
c set par to the closer endpoint.
c
par = amax1(par,parl)
par = amin1(par,paru)
if (par .eq. zero) par = gnorm/dxnorm
c
c beginning of an iteration.
c
150 continue
iter = iter + 1
c
c evaluate the function at the current value of par.
c
if (par .eq. zero) par = amax1(dwarf,p001*paru)
temp = sqrt(par)
do 160 j = 1, n
wa1(j) = temp*diag(j)
160 continue
call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2)
do 170 j = 1, n
wa2(j) = diag(j)*x(j)
170 continue
dxnorm = enorm(n,wa2)
temp = fp
fp = dxnorm - delta
c
c if the function is small enough, accept the current value
c of par. also test for the exceptional cases where parl
c is zero or the number of iterations has reached 10.
c
if (abs(fp) .le. p1*delta
* .or. parl .eq. zero .and. fp .le. temp
* .and. temp .lt. zero .or. iter .eq. 10) go to 220
c
c compute the newton correction.
c
do 180 j = 1, n
l = ipvt(j)
wa1(j) = diag(l)*(wa2(l)/dxnorm)
180 continue
do 210 j = 1, n
wa1(j) = wa1(j)/sdiag(j)
temp = wa1(j)
jp1 = j + 1
if (n .lt. jp1) go to 200
do 190 i = jp1, n
wa1(i) = wa1(i) - r(i,j)*temp
190 continue
200 continue
210 continue
temp = enorm(n,wa1)
parc = ((fp/delta)/temp)/temp
c
c depending on the sign of the function, update parl or paru.
c
if (fp .gt. zero) parl = amax1(parl,par)
if (fp .lt. zero) paru = amin1(paru,par)
c
c compute an improved estimate for par.
c
par = amax1(parl,par+parc)
c
c end of an iteration.
c
go to 150
220 continue
c
c termination.
c
if (iter .eq. 0) par = zero
return
c
c last card of subroutine lmpar.
c
end
subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
integer m,n,lda,lipvt
integer ipvt(lipvt)
logical pivot
real a(lda,n),rdiag(n),acnorm(n),wa(n)
c **********
c
c subroutine qrfac
c
c this subroutine uses householder transformations with column
c pivoting (optional) to compute a qr factorization of the
c m by n matrix a. that is, qrfac determines an orthogonal
c matrix q, a permutation matrix p, and an upper trapezoidal
c matrix r with diagonal elements of nonincreasing magnitude,
c such that a*p = q*r. the householder transformation for
c column k, k = 1,2,...,min(m,n), is of the form
c
c t
c i - (1/u(k))*u*u
c
c where u has zeros in the first k-1 positions. the form of
c this transformation and the method of pivoting first
c appeared in the corresponding linpack subroutine.
c
c the subroutine statement is
c
c subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
c
c where
c
c m is a positive integer input variable set to the number
c of rows of a.
c
c n is a positive integer input variable set to the number
c of columns of a.
c
c a is an m by n array. on input a contains the matrix for
c which the qr factorization is to be computed. on output
c the strict upper trapezoidal part of a contains the strict
c upper trapezoidal part of r, and the lower trapezoidal
c part of a contains a factored form of q (the non-trivial
c elements of the u vectors described above).
c
c lda is a positive integer input variable not less than m
c which specifies the leading dimension of the array a.
c
c pivot is a logical input variable. if pivot is set true,
c then column pivoting is enforced. if pivot is set false,
c then no column pivoting is done.
c
c ipvt is an integer output array of length lipvt. ipvt
c defines the permutation matrix p such that a*p = q*r.
c column j of p is column ipvt(j) of the identity matrix.
c if pivot is false, ipvt is not referenced.
c
c lipvt is a positive integer input variable. if pivot is false,
c then lipvt may be as small as 1. if pivot is true, then
c lipvt must be at least n.
c
c rdiag is an output array of length n which contains the
c diagonal elements of r.
c
c acnorm is an output array of length n which contains the
c norms of the corresponding columns of the input matrix a.
c if this information is not needed, then acnorm can coincide
c with rdiag.
c
c wa is a work array of length n. if pivot is false, then wa
c can coincide with rdiag.
c
c subprograms called
c
c minpack-supplied ... spmpar,enorm
c
c fortran-supplied ... amax1,sqrt,min0
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,jp1,k,kmax,minmn
real ajnorm,epsmch,one,p05,sum,temp,zero
real spmpar,enorm
data one,p05,zero /1.0e0,5.0e-2,0.0e0/
c
c epsmch is the machine precision.
c
epsmch = spmpar(1)
c epsmch=1.18999999e-07
c
c compute the initial column norms and initialize several arrays.
c
c write(*,*) a(1,1)
c a(1,1)=12.0003815
c write(*,*) a(1,1)
c write(*,*) enorm(m,a(1,1)),a(1,1)
c call exit
do 10 j = 1, n
c write(*,*) m,a(1,j),enorm(m,a(1,j))
acnorm(j) = enorm(m,a(1,j))
c write(*,*) acnorm(j)
rdiag(j) = acnorm(j)
c write(*,*) "rdiag",rdiag(j)
wa(j) = rdiag(j)
if (pivot) ipvt(j) = j
10 continue
c
c reduce a to r with householder transformations.
c
c call exit
minmn = min0(m,n)
do 110 j = 1, minmn
if (.not.pivot) go to 40
c
c bring the column of largest norm into the pivot position.
c
c if (j .eq. 9) then
c write(*,*) "a(j,j)",a(j,j)
c call exit
c end if
kmax = j
do 20 k = j, n
c write(*,*) k,kmax,rdiag(k),rdiag(kmax)
if (rdiag(k) .gt. rdiag(kmax)) kmax = k
20 continue
c write(*,*) "kmax",kmax
c if (j .eq. 12) then
c kmax=88
c call exit
c end if
if (kmax .eq. j) go to 40
do 30 i = 1, m
c write(*,*) a(i,j),a(i,kmax)
temp = a(i,j)
a(i,j) = a(i,kmax)
a(i,kmax) = temp
30 continue
c write(*,*) "a(j,j)_2",a(j,j)
c if (j .eq. 13) then
c write(*,*) "a(j,j)_2",a(j,j)
c call exit
c end if
rdiag(kmax) = rdiag(j)
wa(kmax) = wa(j)
k = ipvt(j)
ipvt(j) = ipvt(kmax)
ipvt(kmax) = k
40 continue
c
c compute the householder transformation to reduce the
c j-th column of a to a multiple of the j-th unit vector.
c
ajnorm = enorm(m-j+1,a(j,j))
c write(*,*) "ajnorm",ajnorm,a(j,j)
if (ajnorm .eq. zero) go to 100
if (a(j,j) .lt. zero) ajnorm = -ajnorm
c write(*,*) "ajnorm",ajnorm
do 50 i = j, m
a(i,j) = a(i,j)/ajnorm
c if (j .eq. 13) then
c write(*,*) "a(12,13)",a(12,13)
c call exit
c end if
50 continue
a(j,j) = a(j,j) + one
c call exit
c
c apply the transformation to the remaining columns
c and update the norms.
c
jp1 = j + 1
if (n .lt. jp1) go to 100
do 90 k = jp1, n
sum = zero
do 60 i = j, m
sum = sum + a(i,j)*a(i,k)
60 continue
temp = sum/a(j,j)
do 70 i = j, m
a(i,k) = a(i,k) - temp*a(i,j)
70 continue
if (.not.pivot .or. rdiag(k) .eq. zero) go to 80
temp = a(j,k)/rdiag(k)
rdiag(k) = rdiag(k)*sqrt(amax1(zero,one-temp**2))
if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80
rdiag(k) = enorm(m-j,a(jp1,k))
wa(k) = rdiag(k)
80 continue
90 continue
100 continue
rdiag(j) = -ajnorm
110 continue
c write(*,*) "a_3",a(12,13)
c call exit
return
c
c last card of subroutine qrfac.
c
end
subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
integer n,ldr
integer ipvt(n)
real r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n)
c **********
c
c subroutine qrsolv
c
c given an m by n matrix a, an n by n diagonal matrix d,
c and an m-vector b, the problem is to determine an x which
c solves the system
c
c a*x = b , d*x = 0 ,
c
c in the least squares sense.
c
c this subroutine completes the solution of the problem
c if it is provided with the necessary information from the
c qr factorization, with column pivoting, of a. that is, if
c a*p = q*r, where p is a permutation matrix, q has orthogonal
c columns, and r is an upper triangular matrix with diagonal
c elements of nonincreasing magnitude, then qrsolv expects
c the full upper triangle of r, the permutation matrix p,
c and the first n components of (q transpose)*b. the system
c a*x = b, d*x = 0, is then equivalent to
c
c t t
c r*z = q *b , p *d*p*z = 0 ,
c
c where x = p*z. if this system does not have full rank,
c then a least squares solution is obtained. on output qrsolv
c also provides an upper triangular matrix s such that
c
c t t t
c p *(a *a + d*d)*p = s *s .
c
c s is computed within qrsolv and may be of separate interest.
c
c the subroutine statement is
c
c subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
c
c where
c
c n is a positive integer input variable set to the order of r.
c
c r is an n by n array. on input the full upper triangle
c must contain the full upper triangle of the matrix r.
c on output the full upper triangle is unaltered, and the
c strict lower triangle contains the strict upper triangle
c (transposed) of the upper triangular matrix s.
c
c ldr is a positive integer input variable not less than n
c which specifies the leading dimension of the array r.
c
c ipvt is an integer input array of length n which defines the
c permutation matrix p such that a*p = q*r. column j of p
c is column ipvt(j) of the identity matrix.
c
c diag is an input array of length n which must contain the
c diagonal elements of the matrix d.
c
c qtb is an input array of length n which must contain the first
c n elements of the vector (q transpose)*b.
c
c x is an output array of length n which contains the least
c squares solution of the system a*x = b, d*x = 0.
c
c sdiag is an output array of length n which contains the
c diagonal elements of the upper triangular matrix s.
c
c wa is a work array of length n.
c
c subprograms called
c
c fortran-supplied ... abs,sqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,jp1,k,kp1,l,nsing
real cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero
data p5,p25,zero /5.0e-1,2.5e-1,0.0e0/
c
c copy r and (q transpose)*b to preserve input and initialize s.
c in particular, save the diagonal elements of r in x.
c
do 20 j = 1, n
do 10 i = j, n
r(i,j) = r(j,i)
10 continue
x(j) = r(j,j)
wa(j) = qtb(j)
20 continue
c
c eliminate the diagonal matrix d using a givens rotation.
c
do 100 j = 1, n
c
c prepare the row of d to be eliminated, locating the
c diagonal element using p from the qr factorization.
c
l = ipvt(j)
if (diag(l) .eq. zero) go to 90
do 30 k = j, n
sdiag(k) = zero
30 continue
sdiag(j) = diag(l)
c
c the transformations to eliminate the row of d
c modify only a single element of (q transpose)*b
c beyond the first n, which is initially zero.
c
qtbpj = zero
do 80 k = j, n
c
c determine a givens rotation which eliminates the
c appropriate element in the current row of d.
c
if (sdiag(k) .eq. zero) go to 70
if (abs(r(k,k)) .ge. abs(sdiag(k))) go to 40
cotan = r(k,k)/sdiag(k)
sin = p5/sqrt(p25+p25*cotan**2)
cos = sin*cotan
go to 50
40 continue
tan = sdiag(k)/r(k,k)
cos = p5/sqrt(p25+p25*tan**2)
sin = cos*tan
50 continue
c
c compute the modified diagonal element of r and
c the modified element of ((q transpose)*b,0).
c
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
c
c accumulate the tranformation in the row of s.
c
kp1 = k + 1
if (n .lt. kp1) go to 70
do 60 i = kp1, n
temp = cos*r(i,k) + sin*sdiag(i)
sdiag(i) = -sin*r(i,k) + cos*sdiag(i)
r(i,k) = temp
60 continue
70 continue
80 continue
90 continue
c
c store the diagonal element of s and restore
c the corresponding diagonal element of r.
c
sdiag(j) = r(j,j)
r(j,j) = x(j)
100 continue
c
c solve the triangular system for z. if the system is
c singular, then obtain a least squares solution.
c
nsing = n
do 110 j = 1, n
if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1
if (nsing .lt. n) wa(j) = zero
110 continue
if (nsing .lt. 1) go to 150
do 140 k = 1, nsing
j = nsing - k + 1
sum = zero
jp1 = j + 1
if (nsing .lt. jp1) go to 130
do 120 i = jp1, nsing
sum = sum + r(i,j)*wa(i)
120 continue
130 continue
wa(j) = (wa(j) - sum)/sdiag(j)
140 continue
150 continue
c
c permute the components of z back to components of x.
c
do 160 j = 1, n
l = ipvt(j)
x(l) = wa(j)
160 continue
return
c
c last card of subroutine qrsolv.
c
end
real function spmpar(i)
integer i,j
if(i.eq.1) then
j=4
else if(i.eq.2) then
j=1
else if(i.eq.3) then
j=2
else
stop
end if
spmpar=r1mach(j)
c write(*,*) spmpar,j
c call exit
return
end