DDJAC
Approximates the Jacobian of m functions in n unknowns using divided differences.
Required Arguments
FCN — User-supplied subroutine to evaluate functions. The usage is
CALL FCN (INDX, Y, F[,…]) where
Required Arguments
INDX — Index of the variable whose derivative is to be computed. (Input)
DDJAC will set this argument to the index of the variable whose derivative is being computed. In those cases where there is a mix of finite differencing taking place along with additional analytic terms being computed, (see METHOD = 2), DDJAC will make two calls to FCN each time a new function evaluation is needed, once with INDX positive and a second time with INDX negative.
Y — Array containing the point at which the function is to be computed. (Input)
F — Array of length M, where M is the number of functions to be evaluated at point Y, containing the function values of the equations at point Y. (Output)
Normally, the user will return the values of the functions evaluated at point Y in F. However, when the function can be broken into two parts, a part which is known analytically and a part to be differenced, FCN will be called by DDJAC once with INDX positive for the portion to be differenced and again with INDX negative for the portion which is known analytically. In the case where METHOD=2 has been chosen, FCN must be writtten to handle the known analytic portion separate from the part to be differenced. (See Example 4 for an example where METHOD=2 is used.)
Optional Arguments
FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a detailed description of this argument see FCN_DATA below. (Input/Output)
FCN must be declared EXTERNAL in the calling program.
Y — Array of length N containing the point at which the Jacobian is to be evaluated. (Input)
F — Array of length M containing the function values of the equations at point Y. (Output)
FJAC — Two dimensional array of which the first M by N subarray contains the estimated Jacobian. (Input/Output)
On input the user may set entries of columns that are to be accumulated to initial values (See the optional argument METHOD). On final output, FJAC will contain the estimated Jacobian.
Optional Arguments
M — The number of equations. (Input)
Default: M = SIZE (F).
N — The number of variables. (Input)
Default: N = SIZE (Y).
YSCALE — Array of length N containing the diagonal scaling matrix for the variables. (Input)
YSCALE can also be used to provide appropriate signs for the increments.
Default: YSCALE = 1.0.
METHOD — Array of length N containing the methods used to compute the derivatives. (Input)
METHOD(i) is the method to be used for the i-th variable. METHOD(i) can be one of the values in the following table:
|
Value |
Description |
|
0 |
Indicates one-sided differences. |
|
1 |
Indicates central differences. |
|
Indicates the accumulation of the result from whatever type of differences have been specified previously into initial values of the Jacobian |
|
|
3 |
Indicates a variable is to be skipped |
Default: One-sided differences are used for all variables.
FACTOR — Array of length N containing the percentage factor for differencing. (Input)
For each divided difference for variable j the increment used is del. The value of del is computed as follows: First define σ = sign(YSCALE(j)). If the user has set the elements of array YSCALE to non-default values, then define ya = ∣YSCALE(j)∣. Otherwise, ya = ∣Y(j)∣ and σ = 1. Finally, compute del = σya FACTOR(j). By changing the sign of YSCALE(j), the difference del can have any desired orientation, such as staying within bounds on variable j. For central differences, a reduced factor is used for del that normally results in relative errors as small as machine precision to the 2/3 power. The elements of FACTOR must be such that machine precision to the 3/4 power <= FACTOR(j) <= 0.1
Default: All elements of FACTOR are set to sqrt(machine precision).
ISTATUS — Array of length 10 which contains status information that might prove useful to the user wanting to gain better control over the differencing parameters. (Output)
This information can often be ignored. The following table describes the diagnostic information which is returned in each of the entries of ISTATUS:
|
index |
Description |
|
1 |
The number of times a function evaluation was computed. |
|
2 |
The number of columns in which three attempts were made to increase a percentage factor for differencing (i.e. a component in the FACTOR array) but the computed del remained unacceptably small relative to Y[j] or YSCALE[j]. In such cases the percentage factor is set to the square root of machine precision. |
|
3 |
The number of columns in which the computed del was zero to machine precision because Y[j] or YSCALE[j] was zero. In such cases del is set to the square root of machine precision. |
|
4 |
The number of Jacobian columns which had to be recomputed because the largest difference formed in the column was close to zero relative to scale, where scale = max(∣fi(y)∣),∣fi(y + del × ej)∣) and i denotes the row index of the largest difference in the column currently being processed. index = 10 gives the last column where this occurred. |
|
5 |
The number of columns whose largest difference is close to zero relative to scale after the column has been recomputed. |
|
6 |
The number of times scale information was not available for use in the roundoff and truncation error tests. This occurs when min(∣fi(y)∣),∣fi(y + del × ej)∣) = 0 Where i is the index of the largest difference for the column currently being processed. |
|
7 |
The number of times the increment for differencing (del) was computed and had to be increased because (YSCALE[j] + del) ‑ YSCALE[j]) was too small relative to Y[j] or YSCALE[j]. |
|
8 |
The number of times a component of the FACTOR array was reduced because changes in function values were large and excess truncation error was suspected. index = 9 gives the last column in which this occurred. |
|
9 |
The index of the last column where the corresponding component of the FACTOR array had to be reduced because excessive truncation error was suspected. |
|
10 |
The index of the last column where the difference was small and the column had to be recomputed with an adjusted increment (see index = 4). The largest derivative in this column may be inaccurate due to excessive roundoff error. |
FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied subroutine. (Input/Output)
The derived type, s_fcn_data, is defined as:
type s_fcn_data
real(kind(1e0)), pointer, dimension(:) :: rdata
integer, pointer, dimension(:) :: idata
end type
in module mp_types. The double precision counterpart to s_fcn_data is named d_fcn_data. The user must include a use mp_types statement in the calling program to define this derived type.
FORTRAN 90 Interface
Generic: CALL DDJAC (FCN,Y, F, FJAC [, …])
Specific: The specific interface names are S_DDJAC and D_DDJAC.
Description
Computes the Jacobian matrix for a function f(y) with m components in n independent variables. DDJAC uses divided finite differences to compute the Jacobian. This subroutine is designed for use in numerical methods for solving nonlinear problems where a Jacobian is evaluated repeatedly at neighboring arguments. For example this occurs in a Gauss-Newton method for solving non-linear least squares problems or a non-linear optimization method.
DDJAC is suited for applications where the Jacobian is a dense matrix. All cases m < n, m = n, or m > n are allowed. Both one-sided and central divided differences can be used.
The design allows for computation of derivatives in a variety of contexts. Note that a gradient should be considered as the special case with m = 1, n ≥ 1. A derivative of a single function of one variable is the case m = 1, n = 1. Any non-linear solving routine that optionally requests a Jacobian or gradient can use DDJAC. This should be considered if there are special properties or scaling issues associated with f(y). Use the argument METHOD to specify different differencing options for numerical differentiation. These can be combined with some analytic subexpressions or other known relationships.
The divided differences are computed using values of the independent variables at the initial point yj = y, and differenced points ye = y + del × ej . Here the ej, j = 1,…,n, are the unit coordinate vectors.
The value for each difference del depends on the variable j, the differencing method, and the scaling for that variable. This difference is computed internally. See FACTOR for computational details. The evaluation of f(ye) is normally done by the user-provided argument FCN, using the values ye. The index j, values ye, and output F are arguments to FCN.
The computational kernel of DDJAC performs the following steps:
1. Evaluates the equations at the point Y using FCN.
2. Computes the Jacobian.
3. Computes the difference at ye.
There are four examples provided which illustrate various ways to use DDJAC. A discussion of the expected errors for the difference methods is found in A First Course in Numerical Analysis, Anthony Ralston, McGraw-Hill, NY, (1965).
Examples
In this example, the Jacobian matrix of
is estimated by the finite-difference method at the point (1.0, 1.0).
USE DDJAC_INT
USE WRRRN_INT
IMPLICIT NONE
INTEGER, PARAMETER :: N=2, M=2
REAL FJAC(M,N), Y(N), F(M)
EXTERNAL FCN
DATA Y/2*1.0/
! Get Jacobian one-sided difference
! approximation
CALL DDJAC (FCN, Y, F, FJAC)
CALL WRRRN ("The Jacobian is:", FJAC)
END
SUBROUTINE FCN (INDX, Y, F)
INTEGER INDX
REAL Y(*), F(*)
F(1) = Y(1)*Y(2) - 2.0
F(2) = Y(1) - Y(1)*Y(2) + 1.0
RETURN
END
The Jacobian is:
1 2
1 1.000 1.000
2 0.000 -1.000
A simple use of DDJAC is shown. The gradient of the function 
is required for values 
.
The analytic gradient of this function is:
USE DDJAC_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER, PARAMETER :: N=2, M=1
INTEGER J, NOUT
REAL FJAC(M,N), Y(N), F(M), SCALE(N)
EXTERNAL FCN
DATA Y/2.1, 3.2/ SCALE/1.0, 8000.0/
! Get Gradient one-sided difference
! approximation
CALL DDJAC (FCN, Y, F, FJAC, YSCALE=SCALE)
! Print results
CALL UMACH (2, NOUT)
WRITE (NOUT,99999) (FJAC(1,J),J=1,N)
99999 FORMAT (' The Numerical Gradient is (', 2e15.4,' )')
END
SUBROUTINE FCN (INDX, Y, F)
INTEGER INDX
REAL A, B, C, Y(*), F(*)
A = 2500000.
B = 3.4
C = 4.5
F(1) = A * EXP (B * Y(1)) + C * Y(1) * Y(2) * Y(2)
RETURN
END
The Numerical Gradient is ( 0.1073E+11 0.9268E+02 )
This example uses the same data as in Example 2. Here we assume that the second component of the gradient is analytically known. Therefore only the first gradient component needs numerical approximation. The input values of array METHOD specify that numerical differentiation with respect to y2 is skipped.
USE DDJAC_INT
USE UMACH_INT
IMPLICIT NONE
INTEGER, PARAMETER :: N=2, M=1
INTEGER J, NOUT, METHOD(2)
REAL FJAC(M,N), Y(N), F(M), SCALE(N)
EXTERNAL FCN
DATA Y/2.1, 3.2/ SCALE/1.0, 8000.0/
! Initialize second component
! of Jacobian since it is
! known analytically and can be
! skipped
FJAC(1,2) = 2.0 * 4.5 * Y(1) * Y(2)
! Set METHOD to skip the second
! component
METHOD(1) = 0
METHOD(2) = 3
! Get Gradient approximation
CALL DDJAC (FCN, Y, F, FJAC, YSCALE=SCALE, METHOD=METHOD)
! Print results
CALL UMACH (2, NOUT)
WRITE (NOUT,99999) (FJAC(1,J),J=1,N)
99999 FORMAT (' The Numerical Gradient is (', 2e15.4,' )')
END
SUBROUTINE FCN (INDX, Y, F)
INTEGER INDX
REAL A, B, C, Y(*), F(*)
A = 2500000.
B = 3.4
C = 4.5
F(1) = A * EXP (B * Y(1)) + C * Y(1) * Y(2) * Y(2)
RETURN
END
The Numerical Gradient is ( 0.1073E+11 0.6048E+02 )
This example uses the same data as in Example 2. An alternate examination of the function
shows that the first term on the right-hand side need be evaluated just when computing the first partial. The additive term 
occurs when computing the partial with respect to y1. Also the first term does not depend on the second variable. Thus the first term can be left out of the function evaluation when computing the partial with respect to y2, potentially avoiding cancellation errors. The input values of array METHOD allow DDJAC to use these facts and obtain greater accuracy using a minimum number of computations of the exponential function
USE DDJAC_INT
USE UMACH_INT
USE MP_TYPES
IMPLICIT NONE
INTEGER, PARAMETER :: N=2, M=1
INTEGER J, NOUT, METHOD(2)
REAL FJAC(M,N), Y(N), F(M), SCALE(N)
REAL, TARGET :: RDATA(3)
TYPE(S_FCN_DATA) USER_DATA
EXTERNAL FCN
DATA Y/2.1, 3.2/ SCALE/1.0, 8000.0/
! Set up to pass some extra
! information to the function
RDATA(1) = 2500000.0
RDATA(2) = 3.4
RDATA(3) = 4.5
USER_DATA%RDATA => RDATA
! Initialize first component
! of function since it is
! known
FJAC(1,1) = 4.5 * Y(2) * Y(2)
! Set METHOD to accumulate for
! part of the first partial,
! one-sided differences for
! the second
METHOD(1) = 2
METHOD(2) = 0
! Get Gradient approximation
CALL DDJAC (FCN, Y, F, FJAC, YSCALE=SCALE, METHOD=METHOD, &
FCN_DATA=USER_DATA)
! Print results
CALL UMACH (2, NOUT)
WRITE (NOUT,99999) (FJAC(1,J),J=1,N)
99999 FORMAT (' The Numerical Gradient is (', 2e15.4,' )')
END
SUBROUTINE FCN (INDX, Y, F, FCN_DATA)
USE MP_TYPES
IMPLICIT NONE
INTEGER INDX
REAL A, B, C, Y(*), F(*)
TYPE(S_FCN_DATA) FCN_DATA
A = FCN_DATA%RDATA(1)
B = FCN_DATA%RDATA(2)
C = FCN_DATA%RDATA(3)
! Handle both the differenced
! part and the part that is
! known analytically for each
! dependent variable
SELECT CASE(INDX)
CASE (1)
F(1)=A*EXP(B*Y(1))
CASE(-1)
F(1)= C*Y(2)**2
CASE(2)
F(1) = C*Y(1)*Y(2)**2
CASE(-2)
F(1)=0
END SELECT
RETURN
END
The Numerical Gradient is ( 0.1073E+11 0.6046E+02 )
