*DECK DFDJC1
SUBROUTINE DFDJC1 (FCN, N, X, FVEC, FJAC, LDFJAC, IFLAG, ML, MU,
+ EPSFCN, WA1, WA2)
C***BEGIN PROLOGUE DFDJC1
C***SUBSIDIARY
C***PURPOSE Subsidiary to DNSQ and DNSQE
C***LIBRARY SLATEC
C***TYPE DOUBLE PRECISION (FDJAC1-S, DFDJC1-D)
C***AUTHOR (UNKNOWN)
C***DESCRIPTION
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 DFDJC1(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
C The value of IFLAG should not be changed by FCN unless
C the user wants to terminate execution of DFDJC1.
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 DFDJC1. 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***SEE ALSO DNSQ, DNSQE
C***ROUTINES CALLED D1MACH
C***REVISION HISTORY (YYMMDD)
C 800301 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890831 Modified array declarations. (WRB)
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 900326 Removed duplicate information from DESCRIPTION section.
C (WRB)
C 900328 Added TYPE section. (WRB)
C***END PROLOGUE DFDJC1
DOUBLE PRECISION D1MACH
INTEGER I, IFLAG, J, K, LDFJAC, ML, MSUM, MU, N
DOUBLE PRECISION EPS, EPSFCN, EPSMCH, FJAC(LDFJAC,*),
1 FVEC(*), H, TEMP, WA1(*), WA2(*), X(*), ZERO
SAVE ZERO
DATA ZERO /0.0D0/
C
C EPSMCH IS THE MACHINE PRECISION.
C
C***FIRST EXECUTABLE STATEMENT DFDJC1
EPSMCH = D1MACH(4)
C
EPS = SQRT(MAX(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*ABS(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*ABS(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*ABS(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)
1 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 DFDJC1.
C
END