SUBROUTINE FDJAC1( &
fcnhyb,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,wa1,wa2)
! This subroutine computes a forward-difference approximation
! to the N by N Jacobian matrix associated with a specified
! problem of N functions in N variables. If the Jacobian has
! a banded form, then function evaluations are saved by only
! approximating the nonzero terms.
!
! The subroutine statement is
!
! subroutine fdjac1(fcnhyb,n,x,fvec,fjac,ldfjac,iflag,ml,mu,epsfcn,
! wa1,wa2)
!
! where
!
! FCNHYB is the name of the user-supplied subroutine which
! calculates the functions. FCNHYB must be declared
! in an external statement in the user calling
! program, and should be written as follows.
!
! subroutine fcnhyb(n,x,fvec,iflag)
! integer n,iflag
! real x(n),fvec(n)
! ----------
! calculate the functions at x and
! return this vector in fvec.
! ----------
! return
! end
!
! The value of IFLAG should not be changed by FCNHYB unless
! the user wants to terminate execution of FDJAC1.
! In this case set IFLAG to a negative integer.
!
! N is a positive integer input variable set to the number
! of functions and variables.
!
! X is an input array of length N.
!
! FVEC is an input array of length N which must contain the
! functions evaluated at X.
!
! FJAC is an output N by N array which contains the
! approximation to the Jacobian matrix evaluated at X.
!
! LDFJAC is a positive integer input variable not less than N
! which specifies the leading dimension of the array FJAC.
!
! IFLAG is an integer variable which can be used to terminate
! the execution of FDJAC1. See description of FCNHYB.
!
! ML is a nonnegative integer input variable which specifies
! the number of subdiagonals within the band of the
! Jacobian matrix. If the Jacobian is not banded, set
! ML to at least N - 1.
!
! EPSFCN is an input variable used in determining a suitable
! step length for the forward-difference approximation. This
! approximation assumes that the relative errors in the
! functions are of the order of EPSFCN. If EPSFCN is less
! than the machine precision, it is assumed that the relative
! errors in the functions are of the order of the machine
! precision.
!
! MU is a nonnegative integer input variable which specifies
! the number of superdiagonals within the band of the
! Jacobian matrix. If the Jacobian is not banded, set
! MU to at least N - 1.
!
! WA1 and WA2 are work arrays of length N. If ML + MU + 1 is at
! least N, then the Jacobian is considered dense, and WA2 is
! not referenced.
!
! Argonne National Laboratory. Minpack project. March 1980.
! Burton S. Garbow, Kenneth E. Hillstrom, Jorge J. More
IMPLICIT NONE
INTEGER n,ldfjac,iflag,ml,mu,i,j,k,msum
real(dp) epsfcn
real(dp) x(n),fvec(n),fjac(ldfjac,n),wa1(n),wa2(n)
real(dp) eps,epsmch,h,temp,zero
EXTERNAL fcnhyb
zero = 0.0D0
! Machine precision
epsmch = spmpar(1)
eps = sqrt(max(epsfcn,epsmch))
msum = ml + mu + 1
if (msum < n) goto 40
! Computation of dense approximate Jacobian.
do j = 1, n
temp = x(j)
h = eps*abs(temp)
if (h == zero) h = eps
x(j) = temp + h
call fcnhyb(n,x,wa1,iflag)
if (iflag < 0) goto 30
x(j) = temp
do i = 1, n
fjac(i,j) = (wa1(i) - fvec(i))/h
end do
end do
30 continue
goto 110
40 continue
! Computation of banded approximate Jacobian.
do k = 1, msum
do j = k, n, msum
wa2(j) = x(j)
h = eps*abs(wa2(j))
if (h == zero) h = eps
x(j) = wa2(j) + h
end do
call fcnhyb(n,x,wa1,iflag)
if (iflag < 0) goto 100
do j = k, n, msum
x(j) = wa2(j)
h = eps*abs(wa2(j))
if (h == zero) h = eps
do i = 1, n
fjac(i,j) = zero
if ((i >= (j - mu)).and.(i <= (j + ml))) &
fjac(i,j) = (wa1(i) - fvec(i))/h
end do
end do
end do
100 continue
110 continue
return
end SUBROUTINE FDJAC1