Replacing calls to E04DGF and E04DKF

E04DGF performs an unconstrained minimisation using function values and first derivatives, E04DKF is just an auxiliary routine. To replace these, PDA_UNCMND would be used, which does not make use of derivatives.

The existing NAG code might look as shown below. Two modules are involved, one controls the NAG fit routine, the other serves it by providing the value and gradient of the merit function (the function to be minimised). Information is passed to the merit function in two ways. Scalars and constant-size arrays are in a common block, while one variable-length array is passed as the USER argument. Its length is passed in IUSER. The controlling function may also call the merit function directly.

      SUBROUTINE DOAFIT( ... )
      INTEGER N
      INTEGER IUSER
      INTEGER ITER, IFAIL
      INTEGER IWORK(N+1)
      REAL USER(IUSER)
      DOUBLE PRECISION X(N), FVAL, FGRAD(N)
      DOUBLE PRECISION WORK(13*N)
      EXTERNAL MERIT
  +-  COMMON / ABLOCK / constant-size arrays
  |   CALL MERIT( 0, N, X, FVAL, FGRAD, 0, IUSER, USER )
  |   IFAIL = 1
  |   CALL E04DGF( N, MERIT, ITER, FVAL, FGRAD, X, IWORK, WORK,
  |  :   IUSER, USER, IFAIL )
  |   IF ( IFAIL .NE. 0 ) THEN
  |      An error has occurred
  |   END IF
  |   END
  |
  |   SUBROUTINE MERIT( MODE, N, X, FVAL, FGRAD, NSTATE, IUSER, USER )
  |   INTEGER MODE, N, NSTATE
  |   INTEGER IUSER
  |   REAL USER(IUSER)
  |   DOUBLE PRECISION X(N), FVAL, FGRAD(N)
  +-  COMMON / ABLOCK / constant-size arrays
      FVAL = ...
      DO 1 I = 1, N
         FGRAD(I) = ...
    1 CONTINUE
      END

PDA_UNCMND has separate arguments for the guess and the fit result. With PDA_UNCMND the merit function is simpler in that it need not calculate the gradient. It also has a simpler interface and cannot be passed a variable-length array as above. Thus a pointer to such an array must be passed in the common block. In order to ease de-referencing this pointer, the merit function is split into two modules. MERIT1 receives the passed arguments from its caller and it receives the common block from the master routine DOAFIT. That apart MERIT1 does nothing but call MERIT2. In this call the array pointer is de-referenced.

      SUBROUTINE DOAFIT( ... )
      INTEGER N
      INTEGER IUSER, POINTR
      INTEGER ITER, IFAIL
      INTEGER IWORK(N+1)
      REAL USER(IUSER)
      DOUBLE PRECISION GUESS(N), FIT(N), FVAL
      DOUBLE PRECISION WORK( N*(N+10) )
      EXTERNAL MERIT1
  +-  COMMON / BBLOCK / constant-size arrays,
  |  :                  IUSER, POINTR
  |   POINTR = %LOC(USER)
  |   CALL MERIT1( N, GUESS, FVAL )
  |   IFAIL = 0
  |   CALL PDA_UNCMND( N, GUESS, MERIT1, FIT, FVAL, IFAIL, WORK, N*(N+10) )
  |   IF ( IFAIL .LT. 0 .OR. IFAIL .GT. 3 ) THEN
  |      An error has occurred
  |   END IF
  |   END
  |
  |   SUBROUTINE MERIT1( N, X, FVAL )
  |   INTEGER N
  |   INTEGER IUSER, POINTR
  |   DOUBLE PRECISION X(N), FVAL
  +-  COMMON / BBLOCK / constant-size arrays,
     :                  IUSER, POINTR
      CALL MERIT2( N, X, FVAL, IUSER, %VAL(POINTR) )
      END

      SUBROUTINE MERIT2( N, X, FVAL, IUSER, USER )
      INTEGER N
      INTEGER IUSER
      REAL USER(IUSER)
      DOUBLE PRECISION X(N), FVAL
      FVAL = ...
      END

An alternative to a common block is to have a subroutine with SAVE variables as a reservoir. One routine can call the reservoir routine to set values and another can call it to retrieve values.

%LOC and %VAL are not standard Fortran 77. The only way around it would be to use a different programming language such as Fortran 90 or C.