source/maxlik/src-Fmatch/divdif.f

   1       subroutine divdif(n,x,func,second)
   2       implicit none
   3 C
   4 C Compute numerically the second derivatives of a function of n variables.
   5 C The lower triangle of the derivative matrix is stored in the vector second.
   6 C
   7       integer n,nf,i,ii,j,icant,uiparm
   8       double precision x(n),xi,xj,f,fxp,fxm,fxpym,fxmyp,dwaf,
   9      & delta /5.0d-4/, ddelta /2.5d-7/
  10       double precision urparm
  11       double precision second(n*(n+1)/2)
  12       external func,ufparm
  13       call func(n,x,nf,f,uiparm,urparm,ufparm)
  14       dwaf=2*f
  15       do i=1,n
  16         ii=icant(i,i)
  17 c        print*, i,i," ii",ii
  18         xi=x(i)
  19         x(i)=xi+delta
  20         call func(n,x,nf,fxp,uiparm,urparm,ufparm)
  21         x(i)=xi-delta
  22         call func(n,x,nf,fxm,uiparm,urparm,ufparm)
  23         x(i)=xi
  24         second(ii)=(fxp+fxm-dwaf)/ddelta
  25       enddo
  26       do i=1,n
  27         do j=1,i-1
  28           xi=x(i)
  29           xj=x(j)
  30           ii=icant(i,j)
  31 c          print *,i,j," ii",ii
  32           x(i)=xi+delta
  33           x(j)=xj-delta
  34           call func(n,x,nf,fxpym,uiparm,urparm,ufparm)
  35           x(i)=xi-delta
  36           x(j)=xj+delta
  37           call func(n,x,nf,fxmyp,uiparm,urparm,ufparm)
  38           second(ii)=0.5d0*((second(icant(i,i))+second(icant(j,j)))
  39      &      -(fxpym+fxmyp-dwaf)/ddelta)
  40           x(i)=xi
  41           x(j)=xj
  42         enddo
  43       enddo
  44       return
  45       end
  46 c------------------------------------------------------------------------------
  47       subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
  48       implicit none
  49       integer n,nf,uiparm
  50       double precision x(2),urparm,f
  51       external ufparm
  52 c      f = x(1)**2+x(2)**2
  53       f=dcos(x(1)-x(2))
  54       return
  55       end
  56 C-----------------------------------------------------------------------
  57       INTEGER FUNCTION ICANT(I,J)
  58       implicit none
  59       integer i,j
  60       IF (I.GE.J) THEN
  61         ICANT=(I*(I-1))/2+J
  62       ELSE
  63         ICANT=(J*(J-1))/2+I
  64       ENDIF
  65       RETURN
  66       END
  67 c------------------------------------------------------------------------
  68       subroutine ufparm
  69       return
  70       end