*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