subroutine hybrd1(fcn,n,x,fvec,tol,info,wa,lwa)
integer n,info,lwa
double precision tol
double precision x(n),fvec(n),wa(lwa)
external fcn
c **********
c
c subroutine hybrd1
c
c the purpose of hybrd1 is to find a zero of a system of
c n nonlinear functions in n variables by a modification
c of the powell hybrid method. this is done by using the
c more general nonlinear equation solver hybrd. the user
c must provide a subroutine which calculates the functions.
c the jacobian is then calculated by a forward-difference
c approximation.
c
c the subroutine statement is
c
c subroutine hybrd1(fcn,n,x,fvec,tol,info,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(n,x,fvec,iflag)
c integer n,iflag
c double precision x(n),fvec(n)
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 hybrd1.
c in this case set iflag to a negative integer.
c
c n is a positive integer input variable set to the number
c of functions and variables.
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 n 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 that the relative error
c between x and the solution is at 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 between x and the solution is at most tol.
c
c info = 2 number of calls to fcn has reached or exceeded
c 200*(n+1).
c
c info = 3 tol is too small. no further improvement in
c the approximate solution x is possible.
c
c info = 4 iteration is not making good progress.
c
c wa is a work array of length lwa.
c
c lwa is a positive integer input variable not less than
c (n*(3*n+13))/2.
c
c subprograms called
c
c user-supplied ...... fcn
c
c minpack-supplied ... hybrd
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer index,j,lr,maxfev,ml,mode,mu,nfev,nprint
double precision epsfcn,factor,one,xtol,zero
data factor,one,zero /1.0d2,1.0d0,0.0d0/
info = 0
c
c check the input parameters for errors.
c
if (n .le. 0 .or. tol .lt. zero .or. lwa .lt. (n*(3*n + 13))/2)
* go to 20
c
c call hybrd.
c
maxfev = 200*(n + 1)
xtol = tol
ml = n - 1
mu = n - 1
epsfcn = zero
mode = 2
do 10 j = 1, n
wa(j) = one
10 continue
nprint = 0
lr = (n*(n + 1))/2
index = 6*n + lr
call hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,wa(1),mode,
* factor,nprint,info,nfev,wa(index+1),n,wa(6*n+1),lr,
* wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1))
if (info .eq. 5) info = 4
20 continue
return
c
c last card of subroutine hybrd1.
c
end
subroutine dogleg(n,r,lr,diag,qtb,delta,x,wa1,wa2)
integer n,lr
double precision delta
double precision r(lr),diag(n),qtb(n),x(n),wa1(n),wa2(n)
c **********
c
c subroutine dogleg
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, the
c problem is to determine the convex combination x of the
c gauss-newton and scaled gradient directions that minimizes
c (a*x - b) in the least squares sense, subject to the
c restriction that the euclidean norm of d*x be at most 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 of a. that is, if a = q*r, where q has
c orthogonal columns and r is an upper triangular matrix,
c then dogleg expects the full upper triangle of r and
c the first n components of (q transpose)*b.
c
c the subroutine statement is
c
c subroutine dogleg(n,r,lr,diag,qtb,delta,x,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 input array of length lr which must contain the upper
c triangular matrix r stored by rows.
c
c lr is a positive integer input variable not less than
c (n*(n+1))/2.
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 x is an output array of length n which contains the desired
c convex combination of the gauss-newton direction and the
c scaled gradient direction.
c
c wa1 and wa2 are work arrays of length n.
c
c subprograms called
c
c minpack-supplied ... dpmpar,enorm
c
c fortran-supplied ... dabs,dmax1,dmin1,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,jj,jp1,k,l
double precision alpha,bnorm,epsmch,gnorm,one,qnorm,sgnorm,sum,
* temp,zero
double precision dpmpar,enorm
data one,zero /1.0d0,0.0d0/
c
c epsmch is the machine precision.
c
epsmch = dpmpar(1)
c
c first, calculate the gauss-newton direction.
c
jj = (n*(n + 1))/2 + 1
do 50 k = 1, n
j = n - k + 1
jp1 = j + 1
jj = jj - k
l = jj + 1
sum = zero
if (n .lt. jp1) go to 20
do 10 i = jp1, n
sum = sum + r(l)*x(i)
l = l + 1
10 continue
20 continue
temp = r(jj)
if (temp .ne. zero) go to 40
l = j
do 30 i = 1, j
temp = dmax1(temp,dabs(r(l)))
l = l + n - i
30 continue
temp = epsmch*temp
if (temp .eq. zero) temp = epsmch
40 continue
x(j) = (qtb(j) - sum)/temp
50 continue
c
c test whether the gauss-newton direction is acceptable.
c
do 60 j = 1, n
wa1(j) = zero
wa2(j) = diag(j)*x(j)
60 continue
qnorm = enorm(n,wa2)
if (qnorm .le. delta) go to 140
c
c the gauss-newton direction is not acceptable.
c next, calculate the scaled gradient direction.
c
l = 1
do 80 j = 1, n
temp = qtb(j)
do 70 i = j, n
wa1(i) = wa1(i) + r(l)*temp
l = l + 1
70 continue
wa1(j) = wa1(j)/diag(j)
80 continue
c
c calculate the norm of the scaled gradient and test for
c the special case in which the scaled gradient is zero.
c
gnorm = enorm(n,wa1)
sgnorm = zero
alpha = delta/qnorm
if (gnorm .eq. zero) go to 120
c
c calculate the point along the scaled gradient
c at which the quadratic is minimized.
c
do 90 j = 1, n
wa1(j) = (wa1(j)/gnorm)/diag(j)
90 continue
l = 1
do 110 j = 1, n
sum = zero
do 100 i = j, n
sum = sum + r(l)*wa1(i)
l = l + 1
100 continue
wa2(j) = sum
110 continue
temp = enorm(n,wa2)
sgnorm = (gnorm/temp)/temp
c
c test whether the scaled gradient direction is acceptable.
c
alpha = zero
if (sgnorm .ge. delta) go to 120
c
c the scaled gradient direction is not acceptable.
c finally, calculate the point along the dogleg
c at which the quadratic is minimized.
c
bnorm = enorm(n,qtb)
temp = (bnorm/gnorm)*(bnorm/qnorm)*(sgnorm/delta)
temp = temp - (delta/qnorm)*(sgnorm/delta)**2
* + dsqrt((temp-(delta/qnorm))**2
* +(one-(delta/qnorm)**2)*(one-(sgnorm/delta)**2))
alpha = ((delta/qnorm)*(one - (sgnorm/delta)**2))/temp
120 continue
c
c form appropriate convex combination of the gauss-newton
c direction and the scaled gradient direction.
c
temp = (one - alpha)*dmin1(sgnorm,delta)
do 130 j = 1, n
x(j) = temp*wa1(j) + alpha*x(j)
130 continue
140 continue
return
c
c last card of subroutine dogleg.
c
end
double precision function dpmpar(i)
integer i
c **********
c
c Function dpmpar
c
c This function provides double precision machine parameters
c when the appropriate set of data statements is activated (by
c removing the c from column 1) and all other data statements are
c rendered inactive. Most of the parameter values were obtained
c from the corresponding Bell Laboratories Port Library function.
c
c The function statement is
c
c double precision function dpmpar(i)
c
c where
c
c i is an integer input variable set to 1, 2, or 3 which
c selects the desired machine parameter. If the machine has
c t base b digits and its smallest and largest exponents are
c emin and emax, respectively, then these parameters are
c
c dpmpar(1) = b**(1 - t), the machine precision,
c
c dpmpar(2) = b**(emin - 1), the smallest magnitude,
c
c dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
c
c Argonne National Laboratory. MINPACK Project. November 1996.
c Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More'
c
c **********
integer mcheps(4)
integer minmag(4)
integer maxmag(4)
double precision dmach(3)
equivalence (dmach(1),mcheps(1))
equivalence (dmach(2),minmag(1))
equivalence (dmach(3),maxmag(1))
c
c Machine constants for the IBM 360/370 series,
c the Amdahl 470/V6, the ICL 2900, the Itel AS/6,
c the Xerox Sigma 5/7/9 and the Sel systems 85/86.
c
c data mcheps(1),mcheps(2) / z34100000, z00000000 /
c data minmag(1),minmag(2) / z00100000, z00000000 /
c data maxmag(1),maxmag(2) / z7fffffff, zffffffff /
c
c Machine constants for the Honeywell 600/6000 series.
c
c data mcheps(1),mcheps(2) / o606400000000, o000000000000 /
c data minmag(1),minmag(2) / o402400000000, o000000000000 /
c data maxmag(1),maxmag(2) / o376777777777, o777777777777 /
c
c Machine constants for the CDC 6000/7000 series.
c
c data mcheps(1) / 15614000000000000000b /
c data mcheps(2) / 15010000000000000000b /
c
c data minmag(1) / 00604000000000000000b /
c data minmag(2) / 00000000000000000000b /
c
c data maxmag(1) / 37767777777777777777b /
c data maxmag(2) / 37167777777777777777b /
c
c Machine constants for the PDP-10 (KA processor).
c
c data mcheps(1),mcheps(2) / "114400000000, "000000000000 /
c data minmag(1),minmag(2) / "033400000000, "000000000000 /
c data maxmag(1),maxmag(2) / "377777777777, "344777777777 /
c
c Machine constants for the PDP-10 (KI processor).
c
c data mcheps(1),mcheps(2) / "104400000000, "000000000000 /
c data minmag(1),minmag(2) / "000400000000, "000000000000 /
c data maxmag(1),maxmag(2) / "377777777777, "377777777777 /
c
c Machine constants for the PDP-11.
c
c data mcheps(1),mcheps(2) / 9472, 0 /
c data mcheps(3),mcheps(4) / 0, 0 /
c
c data minmag(1),minmag(2) / 128, 0 /
c data minmag(3),minmag(4) / 0, 0 /
c
c data maxmag(1),maxmag(2) / 32767, -1 /
c data maxmag(3),maxmag(4) / -1, -1 /
c
c Machine constants for the Burroughs 6700/7700 systems.
c
c data mcheps(1) / o1451000000000000 /
c data mcheps(2) / o0000000000000000 /
c
c data minmag(1) / o1771000000000000 /
c data minmag(2) / o7770000000000000 /
c
c data maxmag(1) / o0777777777777777 /
c data maxmag(2) / o7777777777777777 /
c
c Machine constants for the Burroughs 5700 system.
c
c data mcheps(1) / o1451000000000000 /
c data mcheps(2) / o0000000000000000 /
c
c data minmag(1) / o1771000000000000 /
c data minmag(2) / o0000000000000000 /
c
c data maxmag(1) / o0777777777777777 /
c data maxmag(2) / o0007777777777777 /
c
c Machine constants for the Burroughs 1700 system.
c
c data mcheps(1) / zcc6800000 /
c data mcheps(2) / z000000000 /
c
c data minmag(1) / zc00800000 /
c data minmag(2) / z000000000 /
c
c data maxmag(1) / zdffffffff /
c data maxmag(2) / zfffffffff /
c
c Machine constants for the Univac 1100 series.
c
c data mcheps(1),mcheps(2) / o170640000000, o000000000000 /
c data minmag(1),minmag(2) / o000040000000, o000000000000 /
c data maxmag(1),maxmag(2) / o377777777777, o777777777777 /
c
c Machine constants for the Data General Eclipse S/200.
c
c Note - it may be appropriate to include the following card -
c static dmach(3)
c
c data minmag/20k,3*0/,maxmag/77777k,3*177777k/
c data mcheps/32020k,3*0/
c
c Machine constants for the Harris 220.
c
c data mcheps(1),mcheps(2) / '20000000, '00000334 /
c data minmag(1),minmag(2) / '20000000, '00000201 /
c data maxmag(1),maxmag(2) / '37777777, '37777577 /
c
c Machine constants for the Cray-1.
c
c data mcheps(1) / 0376424000000000000000b /
c data mcheps(2) / 0000000000000000000000b /
c
c data minmag(1) / 0200034000000000000000b /
c data minmag(2) / 0000000000000000000000b /
c
c data maxmag(1) / 0577777777777777777777b /
c data maxmag(2) / 0000007777777777777776b /
c
c Machine constants for the Prime 400.
c
c data mcheps(1),mcheps(2) / :10000000000, :00000000123 /
c data minmag(1),minmag(2) / :10000000000, :00000100000 /
c data maxmag(1),maxmag(2) / :17777777777, :37777677776 /
c
c Machine constants for the VAX-11.
c
c data mcheps(1),mcheps(2) / 9472, 0 /
c data minmag(1),minmag(2) / 128, 0 /
c data maxmag(1),maxmag(2) / -32769, -1 /
c
c Machine constants for IEEE machines.
c
data dmach(1) /2.22044604926d-16/
data dmach(2) /2.22507385852d-308/
data dmach(3) /1.79769313485d+308/
c
dpmpar = dmach(i)
return
c
c Last card of function dpmpar.
c
end
double precision function enorm(n,x)
integer n
double precision 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 double precision 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 ... dabs,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i
double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs,
* x1max,x3max,zero
data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/
s1 = zero
s2 = zero
s3 = zero
x1max = zero
x3max = zero
floatn = n
agiant = rgiant/floatn
do 90 i = 1, n
xabs = dabs(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*dsqrt(s1+(s2/x1max)/x1max)
go to 130
100 continue
if (s2 .eq. zero) go to 110
if (s2 .ge. x3max)
* enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3)))
if (s2 .lt. x3max)
* enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3)))
go to 120
110 continue
enorm = x3max*dsqrt(s3)
120 continue
130 continue
return
c
c last card of function enorm.
c
end
subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
* wa1,wa2)
integer n,ldfjac,iflag,ml,mu
double precision epsfcn
double precision x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
c **********
c
c subroutine fdjac1
c
c this subroutine computes a forward-difference approximation
c to the n by n jacobian matrix associated with a specified
c problem of n functions in n variables. if the jacobian has
c a banded form, then function evaluations are saved by only
c approximating the nonzero terms.
c
c the subroutine statement is
c
c subroutine fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
c wa1,wa2)
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(n,x,fvec,iflag)
c integer n,iflag
c double precision x(n),fvec(n)
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 fdjac1.
c in this case set iflag to a negative integer.
c
c n is a positive integer input variable set to the number
c of functions and variables.
c
c x is an input array of length n.
c
c fvec is an input array of length n which must contain the
c functions evaluated at x.
c
c fjac is an output n 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 n
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 fdjac1. see description of fcn.
c
c ml is a nonnegative integer input variable which specifies
c the number of subdiagonals within the band of the
c jacobian matrix. if the jacobian is not banded, set
c ml to at least n - 1.
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 mu is a nonnegative integer input variable which specifies
c the number of superdiagonals within the band of the
c jacobian matrix. if the jacobian is not banded, set
c mu to at least n - 1.
c
c wa1 and wa2 are work arrays of length n. if ml + mu + 1 is at
c least n, then the jacobian is considered dense, and wa2 is
c not referenced.
c
c subprograms called
c
c minpack-supplied ... dpmpar
c
c fortran-supplied ... dabs,dmax1,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,k,msum
double precision eps,epsmch,h,temp,zero
double precision dpmpar
data zero /0.0d0/
c
c epsmch is the machine precision.
c
epsmch = dpmpar(1)
c
eps = dsqrt(dmax1(epsfcn,epsmch))
msum = ml + mu + 1
if (msum .lt. n) go to 40
c
c computation of dense approximate jacobian.
c
do 20 j = 1, n
temp = x(j)
h = eps*dabs(temp)
if (h .eq. zero) h = eps
x(j) = temp + h
call fcn(n,x,wa1,iflag)
if (iflag .lt. 0) go to 30
x(j) = temp
do 10 i = 1, n
fjac(i,j) = (wa1(i) - fvec(i))/h
10 continue
20 continue
30 continue
go to 110
40 continue
c
c computation of banded approximate jacobian.
c
do 90 k = 1, msum
do 60 j = k, n, msum
wa2(j) = x(j)
h = eps*dabs(wa2(j))
if (h .eq. zero) h = eps
x(j) = wa2(j) + h
60 continue
call fcn(n,x,wa1,iflag)
if (iflag .lt. 0) go to 100
do 80 j = k, n, msum
x(j) = wa2(j)
h = eps*dabs(wa2(j))
if (h .eq. zero) h = eps
do 70 i = 1, n
fjac(i,j) = zero
if (i .ge. j - mu .and. i .le. j + ml)
* fjac(i,j) = (wa1(i) - fvec(i))/h
70 continue
80 continue
90 continue
100 continue
110 continue
return
c
c last card of subroutine fdjac1.
c
end
subroutine hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,diag,
* mode,factor,nprint,info,nfev,fjac,ldfjac,r,lr,
* qtf,wa1,wa2,wa3,wa4)
integer n,maxfev,ml,mu,mode,nprint,info,nfev,ldfjac,lr
double precision xtol,epsfcn,factor
double precision x(n),fvec(n),diag(n),fjac(ldfjac,n),r(lr),
* qtf(n),wa1(n),wa2(n),wa3(n),wa4(n)
external fcn
c **********
c
c subroutine hybrd
c
c the purpose of hybrd is to find a zero of a system of
c n nonlinear functions in n variables by a modification
c of the powell hybrid method. 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 hybrd(fcn,n,x,fvec,xtol,maxfev,ml,mu,epsfcn,
c diag,mode,factor,nprint,info,nfev,fjac,
c ldfjac,r,lr,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(n,x,fvec,iflag)
c integer n,iflag
c double precision x(n),fvec(n)
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 hybrd.
c in this case set iflag to a negative integer.
c
c n is a positive integer input variable set to the number
c of functions and variables.
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 n which contains
c the functions evaluated at the output x.
c
c xtol is a nonnegative input variable. termination
c occurs when the relative error between two consecutive
c iterates is at most xtol.
c
c maxfev is a positive integer input variable. termination
c occurs when the number of calls to fcn is at least maxfev
c by the end of an iteration.
c
c ml is a nonnegative integer input variable which specifies
c the number of subdiagonals within the band of the
c jacobian matrix. if the jacobian is not banded, set
c ml to at least n - 1.
c
c mu is a nonnegative integer input variable which specifies
c the number of superdiagonals within the band of the
c jacobian matrix. if the jacobian is not banded, set
c mu to at least n - 1.
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 relative error between two consecutive iterates
c is at most xtol.
c
c info = 2 number of calls to fcn has reached or exceeded
c maxfev.
c
c info = 3 xtol is too small. no further improvement in
c the approximate solution x is possible.
c
c info = 4 iteration is not making good progress, as
c measured by the improvement from the last
c five jacobian evaluations.
c
c info = 5 iteration is not making good progress, as
c measured by the improvement from the last
c ten iterations.
c
c nfev is an integer output variable set to the number of
c calls to fcn.
c
c fjac is an output n by n array which contains the
c orthogonal matrix q produced by the qr factorization
c of the final approximate jacobian.
c
c ldfjac is a positive integer input variable not less than n
c which specifies the leading dimension of the array fjac.
c
c r is an output array of length lr which contains the
c upper triangular matrix produced by the qr factorization
c of the final approximate jacobian, stored rowwise.
c
c lr is a positive integer input variable not less than
c (n*(n+1))/2.
c
c qtf is an output array of length n which contains
c the vector (q transpose)*fvec.
c
c wa1, wa2, wa3, and wa4 are work arrays of length n.
c
c subprograms called
c
c user-supplied ...... fcn
c
c minpack-supplied ... dogleg,dpmpar,enorm,fdjac1,
c qform,qrfac,r1mpyq,r1updt
c
c fortran-supplied ... dabs,dmax1,dmin1,min0,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,jm1,l,msum,ncfail,ncsuc,nslow1,nslow2
integer iwa(1)
logical jeval,sing
double precision actred,delta,epsmch,fnorm,fnorm1,one,pnorm,
* prered,p1,p5,p001,p0001,ratio,sum,temp,xnorm,
* zero
double precision dpmpar,enorm
data one,p1,p5,p001,p0001,zero
* /1.0d0,1.0d-1,5.0d-1,1.0d-3,1.0d-4,0.0d0/
c
c epsmch is the machine precision.
c
epsmch = dpmpar(1)
c
info = 0
iflag = 0
nfev = 0
c
c check the input parameters for errors.
c
if (n .le. 0 .or. xtol .lt. zero .or. maxfev .le. 0
* .or. ml .lt. 0 .or. mu .lt. 0 .or. factor .le. zero
* .or. ldfjac .lt. n .or. lr .lt. (n*(n + 1))/2) 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
call fcn(n,x,fvec,iflag)
nfev = 1
if (iflag .lt. 0) go to 300
fnorm = enorm(n,fvec)
c
c determine the number of calls to fcn needed to compute
c the jacobian matrix.
c
msum = min0(ml+mu+1,n)
c
c initialize iteration counter and monitors.
c
iter = 1
ncsuc = 0
ncfail = 0
nslow1 = 0
nslow2 = 0
c
c beginning of the outer loop.
c
30 continue
jeval = .true.
c
c calculate the jacobian matrix.
c
iflag = 2
call fdjac1(fcn,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,
* wa2)
nfev = nfev + msum
if (iflag .lt. 0) go to 300
c
c compute the qr factorization of the jacobian.
c
call qrfac(n,n,fjac,ldfjac,.false.,iwa,1,wa1,wa2,wa3)
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 70
if (mode .eq. 2) go to 50
do 40 j = 1, n
diag(j) = wa2(j)
if (wa2(j) .eq. zero) diag(j) = one
40 continue
50 continue
c
c on the first iteration, calculate the norm of the scaled x
c and initialize the step bound delta.
c
do 60 j = 1, n
wa3(j) = diag(j)*x(j)
60 continue
xnorm = enorm(n,wa3)
delta = factor*xnorm
if (delta .eq. zero) delta = factor
70 continue
c
c form (q transpose)*fvec and store in qtf.
c
do 80 i = 1, n
qtf(i) = fvec(i)
80 continue
do 120 j = 1, n
if (fjac(j,j) .eq. zero) go to 110
sum = zero
do 90 i = j, n
sum = sum + fjac(i,j)*qtf(i)
90 continue
temp = -sum/fjac(j,j)
do 100 i = j, n
qtf(i) = qtf(i) + fjac(i,j)*temp
100 continue
110 continue
120 continue
c
c copy the triangular factor of the qr factorization into r.
c
sing = .false.
do 150 j = 1, n
l = j
jm1 = j - 1
if (jm1 .lt. 1) go to 140
do 130 i = 1, jm1
r(l) = fjac(i,j)
l = l + n - i
130 continue
140 continue
r(l) = wa1(j)
if (wa1(j) .eq. zero) sing = .true.
150 continue
c
c accumulate the orthogonal factor in fjac.
c
call qform(n,n,fjac,ldfjac,wa1)
c
c rescale if necessary.
c
if (mode .eq. 2) go to 170
do 160 j = 1, n
diag(j) = dmax1(diag(j),wa2(j))
160 continue
170 continue
c
c beginning of the inner loop.
c
180 continue
c
c if requested, call fcn to enable printing of iterates.
c
if (nprint .le. 0) go to 190
iflag = 0
if (mod(iter-1,nprint) .eq. 0) call fcn(n,x,fvec,iflag)
if (iflag .lt. 0) go to 300
190 continue
c
c determine the direction p.
c
call dogleg(n,r,lr,diag,qtf,delta,wa1,wa2,wa3)
c
c store the direction p and x + p. calculate the norm of p.
c
do 200 j = 1, n
wa1(j) = -wa1(j)
wa2(j) = x(j) + wa1(j)
wa3(j) = diag(j)*wa1(j)
200 continue
pnorm = enorm(n,wa3)
c
c on the first iteration, adjust the initial step bound.
c
if (iter .eq. 1) delta = dmin1(delta,pnorm)
c
c evaluate the function at x + p and calculate its norm.
c
iflag = 1
call fcn(n,wa2,wa4,iflag)
nfev = nfev + 1
if (iflag .lt. 0) go to 300
fnorm1 = enorm(n,wa4)
c
c compute the scaled actual reduction.
c
actred = -one
if (fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2
c
c compute the scaled predicted reduction.
c
l = 1
do 220 i = 1, n
sum = zero
do 210 j = i, n
sum = sum + r(l)*wa1(j)
l = l + 1
210 continue
wa3(i) = qtf(i) + sum
220 continue
temp = enorm(n,wa3)
prered = zero
if (temp .lt. fnorm) prered = one - (temp/fnorm)**2
c
c compute the ratio of the actual to the predicted
c reduction.
c
ratio = zero
if (prered .gt. zero) ratio = actred/prered
c
c update the step bound.
c
if (ratio .ge. p1) go to 230
ncsuc = 0
ncfail = ncfail + 1
delta = p5*delta
go to 240
230 continue
ncfail = 0
ncsuc = ncsuc + 1
if (ratio .ge. p5 .or. ncsuc .gt. 1)
* delta = dmax1(delta,pnorm/p5)
if (dabs(ratio-one) .le. p1) delta = pnorm/p5
240 continue
c
c test for successful iteration.
c
if (ratio .lt. p0001) go to 260
c
c successful iteration. update x, fvec, and their norms.
c
do 250 j = 1, n
x(j) = wa2(j)
wa2(j) = diag(j)*x(j)
fvec(j) = wa4(j)
250 continue
xnorm = enorm(n,wa2)
fnorm = fnorm1
iter = iter + 1
260 continue
c
c determine the progress of the iteration.
c
nslow1 = nslow1 + 1
if (actred .ge. p001) nslow1 = 0
if (jeval) nslow2 = nslow2 + 1
if (actred .ge. p1) nslow2 = 0
c
c test for convergence.
c
if (delta .le. xtol*xnorm .or. fnorm .eq. zero) info = 1
if (info .ne. 0) go to 300
c
c tests for termination and stringent tolerances.
c
if (nfev .ge. maxfev) info = 2
if (p1*dmax1(p1*delta,pnorm) .le. epsmch*xnorm) info = 3
if (nslow2 .eq. 5) info = 4
if (nslow1 .eq. 10) info = 5
if (info .ne. 0) go to 300
c
c criterion for recalculating jacobian approximation
c by forward differences.
c
if (ncfail .eq. 2) go to 290
c
c calculate the rank one modification to the jacobian
c and update qtf if necessary.
c
do 280 j = 1, n
sum = zero
do 270 i = 1, n
sum = sum + fjac(i,j)*wa4(i)
270 continue
wa2(j) = (sum - wa3(j))/pnorm
wa1(j) = diag(j)*((diag(j)*wa1(j))/pnorm)
if (ratio .ge. p0001) qtf(j) = sum
280 continue
c
c compute the qr factorization of the updated jacobian.
c
call r1updt(n,n,r,lr,wa1,wa2,wa3,sing)
call r1mpyq(n,n,fjac,ldfjac,wa2,wa3)
call r1mpyq(1,n,qtf,1,wa2,wa3)
c
c end of the inner loop.
c
jeval = .false.
go to 180
290 continue
c
c end of the outer loop.
c
go to 30
300 continue
c
c termination, either normal or user imposed.
c
if (iflag .lt. 0) info = iflag
iflag = 0
if (nprint .gt. 0) call fcn(n,x,fvec,iflag)
return
c
c last card of subroutine hybrd.
c
end
subroutine qform(m,n,q,ldq,wa)
integer m,n,ldq
double precision q(ldq,m),wa(m)
c **********
c
c subroutine qform
c
c this subroutine proceeds from the computed qr factorization of
c an m by n matrix a to accumulate the m by m orthogonal matrix
c q from its factored form.
c
c the subroutine statement is
c
c subroutine qform(m,n,q,ldq,wa)
c
c where
c
c m is a positive integer input variable set to the number
c of rows of a and the order of q.
c
c n is a positive integer input variable set to the number
c of columns of a.
c
c q is an m by m array. on input the full lower trapezoid in
c the first min(m,n) columns of q contains the factored form.
c on output q has been accumulated into a square matrix.
c
c ldq is a positive integer input variable not less than m
c which specifies the leading dimension of the array q.
c
c wa is a work array of length m.
c
c subprograms called
c
c fortran-supplied ... 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,jm1,k,l,minmn,np1
double precision one,sum,temp,zero
data one,zero /1.0d0,0.0d0/
c
c zero out upper triangle of q in the first min(m,n) columns.
c
minmn = min0(m,n)
if (minmn .lt. 2) go to 30
do 20 j = 2, minmn
jm1 = j - 1
do 10 i = 1, jm1
q(i,j) = zero
10 continue
20 continue
30 continue
c
c initialize remaining columns to those of the identity matrix.
c
np1 = n + 1
if (m .lt. np1) go to 60
do 50 j = np1, m
do 40 i = 1, m
q(i,j) = zero
40 continue
q(j,j) = one
50 continue
60 continue
c
c accumulate q from its factored form.
c
do 120 l = 1, minmn
k = minmn - l + 1
do 70 i = k, m
wa(i) = q(i,k)
q(i,k) = zero
70 continue
q(k,k) = one
if (wa(k) .eq. zero) go to 110
do 100 j = k, m
sum = zero
do 80 i = k, m
sum = sum + q(i,j)*wa(i)
80 continue
temp = sum/wa(k)
do 90 i = k, m
q(i,j) = q(i,j) - temp*wa(i)
90 continue
100 continue
110 continue
120 continue
return
c
c last card of subroutine qform.
c
end
subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
integer m,n,lda,lipvt
integer ipvt(lipvt)
logical pivot
double precision 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 ... dpmpar,enorm
c
c fortran-supplied ... dmax1,dsqrt,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
double precision ajnorm,epsmch,one,p05,sum,temp,zero
double precision dpmpar,enorm
data one,p05,zero /1.0d0,5.0d-2,0.0d0/
c
c epsmch is the machine precision.
c
epsmch = dpmpar(1)
c
c compute the initial column norms and initialize several arrays.
c
do 10 j = 1, n
acnorm(j) = enorm(m,a(1,j))
rdiag(j) = acnorm(j)
wa(j) = rdiag(j)
if (pivot) ipvt(j) = j
10 continue
c
c reduce a to r with householder transformations.
c
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
kmax = j
do 20 k = j, n
if (rdiag(k) .gt. rdiag(kmax)) kmax = k
20 continue
if (kmax .eq. j) go to 40
do 30 i = 1, m
temp = a(i,j)
a(i,j) = a(i,kmax)
a(i,kmax) = temp
30 continue
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))
if (ajnorm .eq. zero) go to 100
if (a(j,j) .lt. zero) ajnorm = -ajnorm
do 50 i = j, m
a(i,j) = a(i,j)/ajnorm
50 continue
a(j,j) = a(j,j) + one
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)*dsqrt(dmax1(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
return
c
c last card of subroutine qrfac.
c
end
subroutine r1mpyq(m,n,a,lda,v,w)
integer m,n,lda
double precision a(lda,n),v(n),w(n)
c **********
c
c subroutine r1mpyq
c
c given an m by n matrix a, this subroutine computes a*q where
c q is the product of 2*(n - 1) transformations
c
c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
c
c and gv(i), gw(i) are givens rotations in the (i,n) plane which
c eliminate elements in the i-th and n-th planes, respectively.
c q itself is not given, rather the information to recover the
c gv, gw rotations is supplied.
c
c the subroutine statement is
c
c subroutine r1mpyq(m,n,a,lda,v,w)
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 must contain the matrix
c to be postmultiplied by the orthogonal matrix q
c described above. on output a*q has replaced a.
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 v is an input array of length n. v(i) must contain the
c information necessary to recover the givens rotation gv(i)
c described above.
c
c w is an input array of length n. w(i) must contain the
c information necessary to recover the givens rotation gw(i)
c described above.
c
c subroutines called
c
c fortran-supplied ... dabs,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c **********
integer i,j,nmj,nm1
double precision cos,one,sin,temp
data one /1.0d0/
c
c apply the first set of givens rotations to a.
c
nm1 = n - 1
if (nm1 .lt. 1) go to 50
do 20 nmj = 1, nm1
j = n - nmj
if (dabs(v(j)) .gt. one) cos = one/v(j)
if (dabs(v(j)) .gt. one) sin = dsqrt(one-cos**2)
if (dabs(v(j)) .le. one) sin = v(j)
if (dabs(v(j)) .le. one) cos = dsqrt(one-sin**2)
do 10 i = 1, m
temp = cos*a(i,j) - sin*a(i,n)
a(i,n) = sin*a(i,j) + cos*a(i,n)
a(i,j) = temp
10 continue
20 continue
c
c apply the second set of givens rotations to a.
c
do 40 j = 1, nm1
if (dabs(w(j)) .gt. one) cos = one/w(j)
if (dabs(w(j)) .gt. one) sin = dsqrt(one-cos**2)
if (dabs(w(j)) .le. one) sin = w(j)
if (dabs(w(j)) .le. one) cos = dsqrt(one-sin**2)
do 30 i = 1, m
temp = cos*a(i,j) + sin*a(i,n)
a(i,n) = -sin*a(i,j) + cos*a(i,n)
a(i,j) = temp
30 continue
40 continue
50 continue
return
c
c last card of subroutine r1mpyq.
c
end
subroutine r1updt(m,n,s,ls,u,v,w,sing)
integer m,n,ls
logical sing
double precision s(ls),u(m),v(n),w(m)
c **********
c
c subroutine r1updt
c
c given an m by n lower trapezoidal matrix s, an m-vector u,
c and an n-vector v, the problem is to determine an
c orthogonal matrix q such that
c
c t
c (s + u*v )*q
c
c is again lower trapezoidal.
c
c this subroutine determines q as the product of 2*(n - 1)
c transformations
c
c gv(n-1)*...*gv(1)*gw(1)*...*gw(n-1)
c
c where gv(i), gw(i) are givens rotations in the (i,n) plane
c which eliminate elements in the i-th and n-th planes,
c respectively. q itself is not accumulated, rather the
c information to recover the gv, gw rotations is returned.
c
c the subroutine statement is
c
c subroutine r1updt(m,n,s,ls,u,v,w,sing)
c
c where
c
c m is a positive integer input variable set to the number
c of rows of s.
c
c n is a positive integer input variable set to the number
c of columns of s. n must not exceed m.
c
c s is an array of length ls. on input s must contain the lower
c trapezoidal matrix s stored by columns. on output s contains
c the lower trapezoidal matrix produced as described above.
c
c ls is a positive integer input variable not less than
c (n*(2*m-n+1))/2.
c
c u is an input array of length m which must contain the
c vector u.
c
c v is an array of length n. on input v must contain the vector
c v. on output v(i) contains the information necessary to
c recover the givens rotation gv(i) described above.
c
c w is an output array of length m. w(i) contains information
c necessary to recover the givens rotation gw(i) described
c above.
c
c sing is a logical output variable. sing is set true if any
c of the diagonal elements of the output s are zero. otherwise
c sing is set false.
c
c subprograms called
c
c minpack-supplied ... dpmpar
c
c fortran-supplied ... dabs,dsqrt
c
c argonne national laboratory. minpack project. march 1980.
c burton s. garbow, kenneth e. hillstrom, jorge j. more,
c john l. nazareth
c
c **********
integer i,j,jj,l,nmj,nm1
double precision cos,cotan,giant,one,p5,p25,sin,tan,tau,temp,
* zero
double precision dpmpar
data one,p5,p25,zero /1.0d0,5.0d-1,2.5d-1,0.0d0/
c
c giant is the largest magnitude.
c
giant = dpmpar(3)
c
c initialize the diagonal element pointer.
c
jj = (n*(2*m - n + 1))/2 - (m - n)
c
c move the nontrivial part of the last column of s into w.
c
l = jj
do 10 i = n, m
w(i) = s(l)
l = l + 1
10 continue
c
c rotate the vector v into a multiple of the n-th unit vector
c in such a way that a spike is introduced into w.
c
nm1 = n - 1
if (nm1 .lt. 1) go to 70
do 60 nmj = 1, nm1
j = n - nmj
jj = jj - (m - j + 1)
w(j) = zero
if (v(j) .eq. zero) go to 50
c
c determine a givens rotation which eliminates the
c j-th element of v.
c
if (dabs(v(n)) .ge. dabs(v(j))) go to 20
cotan = v(n)/v(j)
sin = p5/dsqrt(p25+p25*cotan**2)
cos = sin*cotan
tau = one
if (dabs(cos)*giant .gt. one) tau = one/cos
go to 30
20 continue
tan = v(j)/v(n)
cos = p5/dsqrt(p25+p25*tan**2)
sin = cos*tan
tau = sin
30 continue
c
c apply the transformation to v and store the information
c necessary to recover the givens rotation.
c
v(n) = sin*v(j) + cos*v(n)
v(j) = tau
c
c apply the transformation to s and extend the spike in w.
c
l = jj
do 40 i = j, m
temp = cos*s(l) - sin*w(i)
w(i) = sin*s(l) + cos*w(i)
s(l) = temp
l = l + 1
40 continue
50 continue
60 continue
70 continue
c
c add the spike from the rank 1 update to w.
c
do 80 i = 1, m
w(i) = w(i) + v(n)*u(i)
80 continue
c
c eliminate the spike.
c
sing = .false.
if (nm1 .lt. 1) go to 140
do 130 j = 1, nm1
if (w(j) .eq. zero) go to 120
c
c determine a givens rotation which eliminates the
c j-th element of the spike.
c
if (dabs(s(jj)) .ge. dabs(w(j))) go to 90
cotan = s(jj)/w(j)
sin = p5/dsqrt(p25+p25*cotan**2)
cos = sin*cotan
tau = one
if (dabs(cos)*giant .gt. one) tau = one/cos
go to 100
90 continue
tan = w(j)/s(jj)
cos = p5/dsqrt(p25+p25*tan**2)
sin = cos*tan
tau = sin
100 continue
c
c apply the transformation to s and reduce the spike in w.
c
l = jj
do 110 i = j, m
temp = cos*s(l) + sin*w(i)
w(i) = -sin*s(l) + cos*w(i)
s(l) = temp
l = l + 1
110 continue
c
c store the information necessary to recover the
c givens rotation.
c
w(j) = tau
120 continue
c
c test for zero diagonal elements in the output s.
c
if (s(jj) .eq. zero) sing = .true.
jj = jj + (m - j + 1)
130 continue
140 continue
c
c move w back into the last column of the output s.
c
l = jj
do 150 i = n, m
s(l) = w(i)
l = l + 1
150 continue
if (s(jj) .eq. zero) sing = .true.
return
c
c last card of subroutine r1updt.
c
end
(
geocities.com/wilfred_tang)