subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,mode,
* factor,nprint,info,nfev,njev,r,lr,qtf,wa1,wa2,
* wa3,wa4)
integer n,ldfjac,maxfev,mode,nprint,info,nfev,njev,lr
real xtol,factor
real x(n),fvec(n),fjac(ldfjac,n),diag(n),r(lr),qtf(n),wa1(n),
* wa2(n),wa3(n),wa4(n)
c **********
c
c subroutine hybrj
c
c the purpose of hybrj 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 and the jacobian.
c
c the subroutine statement is
c
c subroutine hybrj(fcn,n,x,fvec,fjac,ldfjac,xtol,maxfev,diag,
c mode,factor,nprint,info,nfev,njev,r,lr,qtf,
c wa1,wa2,wa3,wa4)
c
c where
c
c fcn is the name of the user-supplied subroutine which
c calculates the functions and the jacobian. fcn must
c be declared in an external statement in the user
c calling program, and should be written as follows.
c
c subroutine fcn(n,x,fvec,fjac,ldfjac,iflag)
c integer n,ldfjac,iflag
c real x(n),fvec(n),fjac(ldfjac,n)
c ----------
c if iflag = 1 calculate the functions at x and
c return this vector in fvec. do not alter fjac.
c if iflag = 2 calculate the jacobian at x and
c return this matrix in fjac. do not alter 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 hybrj.
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 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 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 with iflag = 1
c has reached maxfev.
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. fvec and fjac should not be altered.
c if nprint is not positive, no special calls of fcn
c 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 with iflag = 1 has
c reached 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 with iflag = 1.
c
c njev is an integer output variable set to the number of
c calls to fcn with iflag = 2.
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,spmpar,enorm,
c qform,qrfac,r1mpyq,r1updt
c
c fortran-supplied ... abs,amax1,amin1,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,ncfail,ncsuc,nslow1,nslow2
integer iwa(1)
logical jeval,sing
real actred,delta,epsmch,fnorm,fnorm1,one,pnorm,prered,p1,p5,
* p001,p0001,ratio,sum,temp,xnorm,zero
real spmpar,enorm
data one,p1,p5,p001,p0001,zero
* /1.0e0,1.0e-1,5.0e-1,1.0e-3,1.0e-4,0.0e0/
c
c epsmch is the machine precision.
c
epsmch = spmpar(1)
c
info = 0
iflag = 0
nfev = 0
njev = 0
c
c check the input parameters for errors.
c
if (n .le. 0 .or. ldfjac .lt. n .or. xtol .lt. zero
* .or. maxfev .le. 0 .or. factor .le. zero
* .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,fjac,ldfjac,iflag)
nfev = 1
if (iflag .lt. 0) go to 300
fnorm = enorm(n,fvec)
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 fcn(n,x,fvec,fjac,ldfjac,iflag)
njev = njev + 1
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) = amax1(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,fjac,ldfjac,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 = amin1(delta,pnorm)
c
c evaluate the function at x + p and calculate its norm.
c
iflag = 1
call fcn(n,wa2,wa4,fjac,ldfjac,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 = amax1(delta,pnorm/p5)
if (abs(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*amax1(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.
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,fjac,ldfjac,iflag)
return
c
c last card of subroutine hybrj.
c
end