      SUBROUTINE  DAVIDN( FUNCT,X,N,NDIF )
C
C  ...  6/20/83, 12/19/92
C     MODIFIED  2/15/93
C
      IMPLICIT  REAL*8( A-H,O-Z )
      DIMENSION  X(40), DX(40), G(40), G0(40), Y(40)
      DIMENSION  H(40,40), WRK(40), S(40)
      COMMON     / CCC /  ISW, IPR, ISMT, IDIF
      COMMON     / DDD /  XM , AIC , SD
      EXTERNAL  FUNCT
      DATA        TAU2  /          1.0D-6  /
      DATA  EPS1, EPS2  / 1.0D-6 , 1.0D-6  /
      RAMDA  = 0.5D0
      CONST1 = 1.0D-30
      IPR = 0
C
C          INITIAL ESTIMATE OF INVERSE OF HESSIAN
C
      ICOUNT = 0
 1000 DO 20  I=1,N
      DO 10  J=1,N
   10 H(I,J) = 0.0D0
      S(I)   = 0.0D0
      DX(I)  = 0.0D0
   20 H(I,I) = 1.0D0
      ISW = 0
C
      IF( NDIF.EQ.0 )  CALL  FUNCT( N,X,XM,G,IG )
      IF( NDIF.GE.1 )  CALL  FUNCND( FUNCT,N,X,XM,G,IG )
C
      IF( IPR .GE. 2 )   WRITE( 6,640 )   XM, SD, AIC
      WRITE(6,650) (X(I),I=1,N), XM, (G(I),I=1,N)
C
      ICC = 0
C      ITERATION
 2000 CONTINUE
      ICC = ICC + 1
      IF( ICC .EQ. 1 )  GO TO 120
C
      DO 40  I=1,N
   40 Y(I) = G(I) - G0(I)
      DO 60  I=1,N
      SUM = 0.0D0
      DO 50  J=1,N
   50 SUM = SUM + Y(J)*H(I,J)
   60 WRK(I) = SUM
      S1 = 0.0D0
      S2 = 0.0D0
      DO 70  I=1,N
      S1 = S1 + WRK(I)*Y(I)
   70 S2 = S2 + DX(I) *Y(I)
      IF( S1.LE.CONST1 .OR. S2.LE.CONST1 )  GO TO 900
C     IF( S1 .LE. S2 )   GO TO 100
C
C          UPDATE THE INVERSE OF HESSIAN MATRIX
C
C               ---  BROYDEN-FLETCHER-GOLDFARB-SHANNO TYPE CORRECTION  -
C
  100 CONTINUE
      STEM = S1 / S2 + 1.0D0
      DO 110  I=1,N
      DO 110  J=I,N
      H(I,J) = H(I,J)-(DX(I)*WRK(J)+WRK(I)*DX(J)-DX(I)*DX(J)*STEM)/S2
  110 H(J,I) = H(I,J)
C
  120 CONTINUE
      SS = 0.0D0
      DO 140  I=1,N
      SUM = 0.0D0
      DO 130  J=1,N
  130 SUM = SUM + H(I,J)*G(J)
      SS  = SS + SUM**2
  140 S(I) = -SUM
C
      S1 = 0.0D0
      S2 = 0.0D0
      DO 150  I=1,N
      S1 = S1 + S(I)*G(I)
  150 S2 = S2 + G(I)*G(I)
C     DS2 = DSQRT(S2)
C     GTEM = DABS(S1)/DS2
C     IF( GTEM.LE.TAU1 .AND. DS2.LE.TAU2 )  GO TO 900
      IF( S1.LT.0.0D0 )  GO TO 200
      DO 170  I=1,N
      DO 160  J=1,N
  160 H(I,J) = 0.0D0
      H(I,I) = 1.0D0
  170 S(I) = -S(I)
  200 CONTINUE
C
      ED = XM
C
C          LINEAR  SEARCH
C
      CALL  LINEAR( FUNCT,X,S,RAMDA,ED,N,IG )
C
      IF( IPR .GE. 2 )  WRITE( 6,630 )  RAMDA, ED, SD, AIC
C
      S1 = 0.0D0
      DO 210  I=1,N
      DX(I) = S(I)*RAMDA
      S1 = S1 + DX(I)*DX(I)
      G0(I) = G(I)
  210 X(I) = X(I) + DX(I)
      XMB = XM
      ISW = 0
C
      IF( NDIF.EQ.0 )  CALL  FUNCT( N,X,XM,G,IG )
      IF( NDIF.GE.1 )  CALL  FUNCND( FUNCT,N,X,XM,G,IG )
      WRITE(6,650) (X(I),I=1,N), XM, (G(I),I=1,N)
C
      S2 = 0.D0
      DO 220  I=1,N
  220 S2 = S2 + G(I)*G(I)
      IF( DSQRT(S2) .LT. TAU2 )  GO TO 900
      IF( XMB/XM-1.D0.LT.EPS1 .AND. DSQRT(S1).LT.EPS2 )  GO TO 900
C     IF( ICC .GE. 5 )  GO TO 900
      GO TO 2000
  900 CONTINUE
      IF( IPR .LE. 0 )  RETURN
      WRITE( 6,600 )
      WRITE( 6,610 )  (X(I),I=1,N)
      WRITE( 6,620 )
      WRITE( 6,610 )  (G(I),I=1,N)
C
      ICOUNT  = ICOUNT + 1
      S2 = 0.D0
      DO 910  I=1,N
  910 S2 = S2 + G(I)*G(I)
      IF( S2.GT.1.0D0.AND.ICOUNT.LE.1 )   GO TO 1000
      RETURN
  600 FORMAT( 1H0,'-----  X  -----' )
  610 FORMAT( 1H ,10D13.5 )
  620 FORMAT( 1H0,'***  GRADIENT  ***' )
  630 FORMAT( 1H ,'LAMBDA =',D15.7,3X,'(-1)LOG LIKELIHOOD =',D23.15,
     *        3X,'SD =',D22.15,5X,'AIC =',D23.15 )
  640 FORMAT( 1H ,26X,'(-1)LOG-LIKELIHOOD =',D23.15,3X,'SD =',D22.15,
     *        5X,'AIC =',D23.15 )
  650 FORMAT( 10X,5F12.7 )
      E N D
      SUBROUTINE  FUNCND( FUNCT,M,A,F,G,IFG )
C
C  ...  FUNCTION EVALUATION AND NUMERICAL DIFFERENCING  ...
C     MODIFIED  2/15/93
C
      IMPLICIT   REAL*8( A-H,O-Z )
      DIMENSION  A(M) , G(M) , B(20),GD(5)
      COMMON  / CCC /  ISW , IPR, ISMT, IDIF
      COMMON  /CMFUNC/ DJACOB,FC,SIG2,AIC,FI,SIG2I,AICI,GI(20),GC(20)
C     DATA       ICNT /0/
      CONST = 0.00001D0
C
      CALL  FUNCT( M,A,F,GD,IFG )
      FB = F
      IF( ISW .GE. 1 )   RETURN
C
C     WRITE( 6,600 )   (A(I),I=1,M)
      DO 10  I=1,M
   10 B(I) = A(I)
C
      DO 30  II=1,M
      B(II) = A(II) + CONST
      CALL  FUNCT( M,B,FF,GD,IFG )
      IF( IDIF .EQ. 1 )  GO TO 20
      B(II) = A(II) - CONST
      CALL  FUNCT( M,B,FB,GD,IFG )
   20 G(II) = (FF-FB)/(CONST*IDIF)
      IF( G(II) .GT. 1.0D20 )  G(II) = (F-FB)/CONST
      IF( G(II) .LT.-1.0D20 )  G(II) = (FF-F)/CONST
      IF( FB.GT.F .AND. FF.GT.F )  G(II) = 0.0D0
   30 B(II) = A(II)
C
      RETURN
C
  600 FORMAT( 3X,'---  PARAMETER  ---',(/,3X,5D13.5) )
  610 FORMAT( 3X,'---  GRADIENT  ---',(/,3X,5D13.5) )
      E N D
      SUBROUTINE  LINEAR( FUNCT,X,H,RAM,EE,K,IG )
C
C  ...  LINEAR SEARCH  ...
C     MODIFIED  2/15/93
C
      IMPLICIT  REAL*8( A-H,O-Z )
      INTEGER  RETURN,SUB
      DIMENSION  X(1), H(1), X1(40)
      DIMENSION  G(40)
      COMMON     / CCC /  ISW , IPR, ISMT, IDIF
C     IPR = 7
C     C1 = 1.0D-7
C     C2 = 1.0D-30
C
      ISW = 1
      IG  = 0
      RAM = 0.1D0
      CONST2 = 1.0D-60
   14 DO 15 I=1,K
   15 IF( DABS(H(I)) .GT. 1.0D10 )  GO TO 16
      GO TO 18
   16 DO 17 I=1,K
   17 H(I) = H(I)*1.0D-10
      GO TO 14
   18 CONTINUE
      HNORM = 0.D0
      DO 10  I=1,K
   10 HNORM = HNORM + H(I)**2
      HNORM = DSQRT( HNORM )
C
      RAM2 = RAM
      E1 =EE
      RAM1 = 0.D0
C
      DO 20  I=1,K
   20 X1(I) = X(I) + RAM2*H(I)
      CALL  FUNCT( K,X1,E2,G,IG )
      IF(IPR.GE.7)  WRITE(6,2)  RAM2,E2
C
      IF( IG .EQ. 1  )  GO TO 50
      IF( E2 .GT. E1 )  GO TO 50
   30 RAM3 = RAM2*4.D0
      DO 40  I=1,K
   40 X1(I) = X(I) + RAM3*H(I)
      CALL  FUNCT( K,X1,E3,G,IG )
      IF( IG.EQ.1 )  GO TO  500
      IF( IPR.GE.7 )  WRITE(6,3)  RAM3,E3
      IF( E3 .GT. E2 )  GO TO 70
      IF(RAM3.GT.1.0D10 .AND. E3.LT.E1)  GO TO 45
      IF(RAM3.GT.1.0D10 .AND. E3.GE.E1)  GO TO 46
      RAM1 = RAM2
      RAM2 = RAM3
      E1 = E2
      E2 = E3
      GO TO 30
C
   45 RAM = RAM3
      EE = E3
      RETURN
C
   46 RAM = 0.0D0
      RETURN
C
   50 RAM3 = RAM2
      E3 = E2
      RAM2 = RAM3*0.1D0
      IF( RAM2*HNORM .LT. CONST2 )  GO TO  400
      DO 60  I=1,K
   60 X1(I) = X(I) + RAM2*H(I)
      CALL  FUNCT( K,X1,E2,G,IG )
      IF(IPR.GE.7)  WRITE(6,4)  RAM2,E2
      IF( E2.GT.E1 )  GO TO 50
C
   70 ASSIGN 80 TO RETURN
      GO TO 200
C
   80 DO 90  I=1,K
   90 X1(I) = X(I) + RAM*H(I)
      CALL  FUNCT( K,X1,EE,G,IG )
      IF(IPR.GE.7)  WRITE(6,5)  RAM,EE
C
      IFG = 0
      ASSIGN  300 TO  SUB
      ASSIGN 200 TO SUB
   95 ASSIGN 130 TO RETURN
      IF( RAM .GT. RAM2 )  GO TO 110
      IF( EE .GE. E2 )  GO TO 100
      RAM3 = RAM2
      RAM2 = RAM
      E3 =E2
      E2 =EE
      GO TO  SUB,( 200,300 )
C
  100 RAM1 = RAM
      E1 = EE
      GO TO  SUB,( 200,300 )
C
  110 IF( EE .LE. E2 )  GO TO 120
      RAM3 = RAM
      E3 = EE
      GO TO  SUB,( 200,300 )
C
  120 RAM1 = RAM2
      RAM2 = RAM
      E1 = E2
      E2 = EE
      GO TO  SUB,( 200,300 )
C
  130 DO 140  I=1,K
  140 X1(I) = X(I) + RAM*H(I)
      CALL  FUNCT( K,X1,EE,G,IG )
      IF( IPR.GE.7 )  WRITE(6,6)  RAM,EE
      ASSIGN 200 TO SUB
      IFG = IFG+1
  600 FORMAT( 1H ,6D20.13 )
C     IF( HNORM*(RAM3-RAM1) .GT. 1.0D-12 )  GO TO 95
      IF( RAM3-RAM1 .GT. RAM2*1.0D-1 )  GO TO 95
C     IF( E1-EE .GT. 1.0D-15 )  GO TO 95
C     IF( E3-EE .GT. 1.0D-15 )  GO TO 95
C     IF(DX.LE.C1*XN+C2 .AND. DF.LE.C1*DABS(EE)+C2)  RETURN
C     IF( IFG .LE. 2 )  GO TO 95
C
      IF( E2 .LT. EE )  RAM = RAM2
      RETURN
C
C      -------  INTERNAL SUBROUTINE SUB1  -------
  200 IF( RAM3-RAM2 .GT. 5.0D0*(RAM2-RAM1) )  GO TO 202
      IF( RAM2-RAM1 .GT. 5.0D0*(RAM3-RAM2) )  GO TO 204
      A1 = (RAM3-RAM2)*E1
      A2 = (RAM1-RAM3)*E2
      A3 = (RAM2-RAM1)*E3
      B2 = (A1+A2+A3)*2.D0
      B1 = A1*(RAM3+RAM2) + A2*(RAM1+RAM3) + A3*(RAM2+RAM1)
      IF( B2 .EQ. 0.D0 )  GO TO 210
      RAM = B1 /B2
C
      IF( RAM .LE. RAM1 )  RAM = (RAM1 + RAM2)/2.0D0
      IF( RAM .GE. RAM3 )  RAM = (RAM2 + RAM3)/2.0D0
      IF( DABS(RAM-RAM2) .LE. 1.0D-15 )  RAM = (RAM2*4.D0 + RAM3)/5.0D0
      GO TO RETURN ,( 80,130 )
  202 RAM = (4.0D0*RAM2 + RAM3)/5.0D0
      GO TO RETURN, (80,130)
  204 RAM = (RAM1 + 4.0D0*RAM2)/5.0D0
      GO TO RETURN, (80,130)
C
  210 IG = 1
      RAM = RAM2
      RETURN
C
C      -------  INTERNAL SUBROUTINE SUB2  -------
C
  300 IF( RAM3-RAM2 .GT. RAM2-RAM1 )  GO TO 310
      RAM = (RAM1+RAM2)*0.5D0
      GO TO RETURN ,( 80,130 )
C
  310 RAM = (RAM2+RAM3)*0.5D0
      GO TO RETURN ,( 80,130 )
C
  400 RAM = 0.D0
      RETURN
C
  500 RAM = (RAM2+RAM3)*0.5D0
  510 DO 520  I=1,K
  520 X1(I) = X(I) + RAM*H(I)
      CALL  FUNCT( K,X1,E3,G,IG )
      IF( IPR.GE.7 )  WRITE(6,7)  RAM,E3
      IF( IG.EQ.1 )  GO TO 540
      IF( E3.GT.E2 )  GO TO 530
      RAM1 = RAM2
      RAM2 = RAM
      E1 = E2
      E2 = E3
      GO TO 500
C
  530 RAM3 = RAM
      GO TO 70
C
  540 RAM = (RAM2+RAM)*0.5D0
      GO TO 510
C
C ------------------------------------------------------------
    1 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E1 =',D25.17 )
    2 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E2 =',D25.17 )
    3 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E3 =',D25.17 )
    4 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E4 =',D25.17 )
    5 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E5 =',D25.17 )
    6 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E6 =',D25.17 )
    7 FORMAT( 1H ,'LAMBDA =',D18.10, 10X,'E7 =',D25.17 )
      E N D
