      SUBROUTINE HESCHK(NR,N,X,FCN,D1FCN,D2FCN,F,G,A,TYPSIZ,SX,RNF,
     +     ANALTL,IAGFLG,UDIAG,WRK1,WRK2,MSG,IPR)
C
C PURPOSE
C -------
C CHECK ANALYTIC HESSIAN AGAINST ESTIMATED HESSIAN
C  (THIS MAY BE DONE ONLY IF THE USER SUPPLIED ANALYTIC HESSIAN
C   D2FCN FILLS ONLY THE LOWER TRIANGULAR PART AND DIAGONAL OF A)
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C X(N)         --> ESTIMATE TO A ROOT OF FCN
C FCN          --> NAME OF SUBROUTINE TO EVALUATE OPTIMIZATION FUNCTION
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C                       FCN:  R(N) --> R(1)
C D1FCN        --> NAME OF SUBROUTINE TO EVALUATE GRADIENT OF FCN.
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C D2FCN        --> NAME OF SUBROUTINE TO EVALUATE HESSIAN OF FCN.
C                  MUST BE DECLARED EXTERNAL IN CALLING ROUTINE
C F            --> FUNCTION VALUE:  FCN(X)
C G(N)        <--  GRADIENT:  G(X)
C A(N,N)      <--  ON EXIT:  HESSIAN IN LOWER TRIANGULAR PART AND DIAG
C TYPSIZ(N)    --> TYPICAL SIZE FOR EACH COMPONENT OF X
C SX(N)        --> DIAGONAL SCALING MATRIX:  SX(I)=1./TYPSIZ(I)
C RNF          --> RELATIVE NOISE IN OPTIMIZATION FUNCTION FCN
C ANALTL       --> TOLERANCE FOR COMPARISON OF ESTIMATED AND
C                  ANALYTICAL GRADIENTS
C IAGFLG       --> =1 IF ANALYTIC GRADIENT SUPPLIED
C UDIAG(N)     --> WORKSPACE
C WRK1(N)      --> WORKSPACE
C WRK2(N)      --> WORKSPACE
C MSG         <--> MESSAGE OR ERROR CODE
C                    ON INPUT : IF =1XX DO NOT COMPARE ANAL + EST HESS
C                    ON OUTPUT: =-22, PROBABLE CODING ERROR OF HESSIAN
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C
      DIMENSION X(N),G(N),A(NR,1)
      DIMENSION TYPSIZ(N),SX(N)
      DIMENSION UDIAG(N),WRK1(N),WRK2(N)
      EXTERNAL FCN,D1FCN
C
C COMPUTE FINITE DIFFERENCE APPROXIMATION A TO THE HESSIAN.
C
      IF(IAGFLG.EQ.1) CALL FSTOFD(NR,N,N,X,D1FCN,G,A,SX,RNF,WRK1,3)
      IF(IAGFLG.NE.1) CALL SNDOFD(NR,N,X,FCN,F,A,SX,RNF,WRK1,WRK2)
C
      KER=0
C
C COPY LOWER TRIANGULAR PART OF "A" TO UPPER TRIANGULAR PART
C AND DIAGONAL OF "A" TO UDIAG
C
      DO 30 J=1,N
        UDIAG(J)=A(J,J)
        IF(J.EQ.N) GO TO 30
        JP1=J+1
        DO 25 I=JP1,N
          A(J,I)=A(I,J)
   25   CONTINUE
   30 CONTINUE
C
C COMPUTE ANALYTIC HESSIAN AND COMPARE TO FINITE DIFFERENCE
C APPROXIMATION.
C
      CALL D2FCN(NR,N,X,A)
      DO 40 J=1,N
C***
        HS=MAX(ABS(G(J)),1.0)/MAX(ABS(X(J)),TYPSIZ(J))
        IF(ABS(A(J,J)-UDIAG(J)).GT.MAX(ABS(UDIAG(J)),HS)*ANALTL)
     +       KER=1
        IF(J.EQ.N) GO TO 40
        JP1=J+1
        DO 35 I=JP1,N
          IF(ABS(A(I,J)-A(J,I)).GT.MAX(ABS(A(I,J)),HS)*ANALTL) KER=1
   35   CONTINUE
   40 CONTINUE
C
      IF(KER.EQ.0) GO TO 90
        WRITE(IPR,901)
        DO 50 I=1,N
          IF(I.EQ.1) GO TO 45
          IM1=I-1
          DO 43 J=1,IM1
            WRITE(IPR,902) I,J,A(I,J),A(J,I)
   43     CONTINUE
   45     WRITE(IPR,902) I,I,A(I,I),UDIAG(I)
   50   CONTINUE
        MSG=-22
C     ENDIF
   90 CONTINUE
      RETURN
  901 FORMAT(47H HESCHK    PROBABLE ERROR IN CODING OF ANALYTIC,
     +       18H HESSIAN FUNCTION./
     +       21H HESCHK      ROW  COL,14X,8HANALYTIC,14X,10H(ESTIMATE))
  902 FORMAT(11H HESCHK    ,2I5,2X,E20.13,2X,1H(,E20.13,1H))
      END
