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