CHGRD

Checks a user-supplied gradient of a function.

Required Arguments

FCN — User-supplied subroutine to evaluate the function of which the gradient will be checked. The usage is

CALL FCN (N, X, F), where

CALL FCN (N, X, F), where

N – Length of X. (Input)

X – The point at which the function is evaluated. (Input)

X should not be changed by FCN.

X should not be changed by FCN.

F – The computed function value at the point X. (Output)

FCN must be declared EXTERNAL in the calling program.

GRAD — Vector of length N containing the estimated gradient at X. (Input)

X — Vector of length N containing the point at which the gradient is to be checked. (Input)

INFO — Integer vector of length N. (Output)

INFO(I) = 0 means the user-supplied gradient is a poor estimate of the numerical gradient at the point X(I).

INFO(I) = 1 means the user-supplied gradient is a good estimate of the numerical gradient at the point X(I).

INFO(I) = 2 means the user-supplied gradient disagrees with the numerical gradient at the point X(I), but it might be impossible to calculate the numerical gradient.

INFO(I) = 3 means the user-supplied gradient and the numerical gradient are both zero at X(I), and, therefore, the gradient should be rechecked at a different point.

Optional Arguments

N — Dimension of the problem. (Input)

Default: N = SIZE (X,1).

Default: N = SIZE (X,1).

FORTRAN 90 Interface

Generic: CALL CHGRD (FCN, GRAD, X, INFO [, …])

Specific: The specific interface names are S_CHGRD and D_CHGRD.

FORTRAN 77 Interface

Single: CALL CHGRD (FCN, GRAD, N, X, INFO)

Double: The double precision name is DCHGRD.

Description

The routine CHGRD uses the following finite-difference formula to estimate the gradient of a function of n variables at x:

where hi = ɛ1/2 max{∣xi∣, 1/si} sign(xi), ɛ is the machine epsilon, ei is the i-th unit vector, and si is the scaling factor of the i-th variable.

The routine CHGRD checks the user-supplied gradient ∇f(x) by comparing it with the finite-difference gradient g(x). If

where = ɛ1/4, then (∇f(x))i, which is the i-th element of ∇f(x), is declared correct; otherwise, CHGRD computes the bounds of calculation error and approximation error. When both bounds are too small to account for the difference, (∇f(x))i is reported as incorrect. In the case of a large error bound, CHGRD uses a nearly optimal stepsize to recompute gi(x) and reports that (∇f(x))i is correct if

Otherwise, (∇f(x))i is considered incorrect unless the error bound for the optimal step is greater than ∣(∇f(x))i∣. In this case, the numeric gradient may be impossible to compute correctly. For more details, see Schnabel (1985).

Comments

1. Workspace may be explicitly provided, if desired, by use of C2GRD/DC2GRD. The reference is:

CALL C2GRD (FCN, GRAD, N, X, INFO, FX, XSCALE, EPSFCN, XNEW)

The additional arguments are as follows:

FX — The functional value at X.

XSCALE — Real vector of length N containing the diagonal scaling matrix.

EPSFCN — The relative “noise” of the function FCN.

XNEW — Real work vector of length N.

2. Informational errors

Type | Code | Description |
---|---|---|

4 | 1 | The user-supplied gradient is a poor estimate of the numerical gradient. |

Example

The user-supplied gradient of

at (625, 1, 3.125, 0.25) is checked where t = 2.125.

USE CHGRD_INT

USE WRIRN_INT

IMPLICIT NONE

! Declare variables

INTEGER N

PARAMETER (N=4)

!

INTEGER INFO(N)

REAL GRAD(N), X(N)

EXTERNAL DRIV, FCN

!

! Input values for point X

! X = (625.0, 1.0, 3.125, .25)

!

DATA X/625.0E0, 1.0E0, 3.125E0, 0.25E0/

!

CALL DRIV (N, X, GRAD)

!

CALL CHGRD (FCN, GRAD, X, INFO)

CALL WRIRN (’The information vector’, INFO, 1, N, 1)

!

END

!

SUBROUTINE FCN (N, X, FX)

INTEGER N

REAL X(N), FX

!

REAL EXP

INTRINSIC EXP

!

FX = X(1) + X(2)*EXP(-1.0E0*(2.125E0-X(3))**2/X(4))

RETURN

END

!

SUBROUTINE DRIV (N, X, GRAD)

INTEGER N

REAL X(N), GRAD(N)

!

REAL EXP

INTRINSIC EXP

!

GRAD(1) = 1.0E0

GRAD(2) = EXP(-1.0E0*(2.125E0-X(3))**2/X(4))

GRAD(3) = X(2)*EXP(-1.0E0*(2.125E0-X(3))**2/X(4))*2.0E0/X(4)* &

(2.125-X(3))

GRAD(4) = X(2)*EXP(-1.0E0*(2.125E0-X(3))**2/X(4))* &

(2.125E0-X(3))**2/(X(4)*X(4))

RETURN

END

Output

The information vector

1 2 3 4

1 1 1 1