! Last change: DOS 29 Jul 2000 4:41 pm ! *** copyright 2000 *** ! *** filename ddel12f.f95 *** John F. Monahan ** ! ********************** program ddel12f ! demonstration of del12f -- computes numerical derivatives -- ! gradients and Hessians ! for functions of several variables ! *** double precision version *** ! implicit none integer, parameter :: n = 4 ! dimension of problem real(kind=8), external :: f real(kind=8), dimension(n) :: x,dx,d1a,d1f real(kind=8), dimension( (n*(n+1))/2 ) :: d2a,d2f integer i,j,k ! 21 format(i4,f8.4,f12.6,3x,4f12.6) 25 format(/' first the analytic results'/ & & ' i x(i) delf(i) del2f(i,j)') 26 format(/' now the numerical derivatives, dx''s now',f9.6/ & & ' i x(i) delf(i) del2f(i,j)') ! data x/ 1.d0, 2.d0, 3.d0, 4.d0 / ! deltas for x -- start big data dx/ 4*1.d0 / ! open( unit=6, file='ddel12f.out' ) ! compute analytic derivatives call analyt(x,d1a,d2a) ! write out analytic results write(6,25) do i = 1,4 write(6,21) i,x(i),d1a(i),(d2a( (i*(i-1))/2+j ),j=1,i) end do ! loop on i ! compute numerical derivatives do k = 1,5 ! reduce dx's each time by 4 do i = 1,4 dx(i) = dx(i) / 4. end do ! loop on i ! call del12f call del12f(f,x,4,dx,d1f,d2f) ! write out numerical results write(6,26) dx(1) do i = 1,4 write(6,21) i,x(i),d1f(i),(d2f( (i*(i-1))/2+j ),j=1,i) end do ! loop on i ! end do ! loop on k stop end program ddel12f real(kind=8) function f(x) implicit none real(kind=8), dimension(4), intent(in) :: x real(kind=8) rl ! function we're finding derivatives of rl = log( x(1)*x(2) / ( x(3)*x(4) ) ) f = rl * rl return end function f subroutine analyt(x,d1a,d2a) ! compute the analytic derivatives for log( x1*x2 /(x3*x4) )**2 implicit none real(kind=8), dimension(4), intent(in) :: x real(kind=8), dimension(4), intent(out) :: d1a real(kind=8), dimension(10), intent(out) :: d2a real(kind=8) rl rl = log( x(1)*x(2) / ( x(3)*x(4) ) ) ! start with the first partials d1a(1) = 2.*rl/x(1) d1a(2) = 2.*rl/x(2) d1a(3) = -2.*rl/x(3) d1a(4) = -2.*rl/x(4) ! now the hessian matrix d2a(1) = 2.*(1.-rl)/(x(1)*x(1)) d2a(2) = 2./( x(1)*x(2) ) d2a(3) = 2.*(1.-rl)/(x(2)*x(2)) d2a(4) = -2./( x(1)*x(3) ) d2a(5) = -2./( x(2)*x(3) ) d2a(6) = 2.*(1.+rl)/(x(3)*x(3)) d2a(7) = -2./( x(1)*x(4) ) d2a(8) = -2./( x(2)*x(4) ) d2a(9) = 2./( x(3)*x(4) ) d2a(10) = 2.*(1.+rl)/(x(4)*x(4)) return end subroutine analyt SUBROUTINE DEL12F(F,X,N,DX,D1F,D2F) ! COMPUTES NUMERICAL FIRST AND SECOND DERIVATIVES FOR FUNCTION F ! GRADIENT VECTOR AND HESSIAN MATRIX ! ! F FUNCTION OF INTEREST ! X REAL V POINT AT WHICH GRADIENT AND HESSIAN ARE TO BE COMPUTED ! N INTEGER DIMENSION OF X ! DX REAL V CHANGES IN EACH COMPONENT OF X ! D1F REAL V ON OUTPUT, THE GRADIENT VECTOR ! D2F REAL V ON OUTPUT, THE HESSIAN MATRIX IN SYMMETRIC STORAGE ! IMPLICIT NONE REAL(KIND=8) F INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: X REAL(KIND=8), DIMENSION( N ), INTENT(IN) :: DX REAL(KIND=8), DIMENSION( N ), INTENT(OUT) :: D1F REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(OUT) :: D2F REAL(KIND=8) V0,VP,VN,XH,XG INTEGER I,J,K,KI,KJ ! CENTER POINT V0 = F(X) ! GET CENTERED FIRST DIFFERENCES DO I = 1,N XH = X(I) X(I) = XH + DX(I) D1F(I) = F(X) ! STORE F(X+H*EI) HERE X(I) = XH - DX(I) KI = (I*(I+1))/2 D2F(KI) = F(X) ! STORE F(X-H*EI) HERE X(I) = XH END DO ! LOOP ON I ! ! HAVE AXIS POINTS STORED, NOW DO THE HESSIAN DO I = 1,N XG = X(I) KI = (I*(I+1))/2 ! COMMON FOR EACH I DO J = I,N ! SKIP THE SECOND PARTIALS FOR LAST IF( I .NE. J ) THEN XH = X(J) X(I) = XG + DX(I) X(J) = XH + DX(J) ! EVALUATED AT THIS CORNER ON THE IJ FACE VP = F(X) ! X(I) = XG - DX(I) X(J) = XH - DX(J) ! EVALUATE AT THE NEGATIVE CORNER ON THE IJ FACE VN = F(X) ! X(J) = XH KJ = (J*(J+1))/2 K = KJ - J + I D2F(K) = ( ( (VP-D1F(I)) - (D1F(J)-V0) ) + ( (VN-D2F(KI)) - & & (D2F(KJ)-V0) ) ) / (2.D0*DX(I)*DX(J)) ! ENDIF END DO ! LOOP ON J ! WHEN DONE WITH ALL I CASES, GET GRADIENT AND HESSIAN VP = D1F(I) VN = D2F(KI) D1F(I) = (VP-VN)/(2.D0*DX(I)) D2F(KI) = ( (VP-V0) - (V0-VN) ) / (DX(I)*DX(I)) X(I) = XG ! END DO ! LOOP ON I ! RETURN END SUBROUTINE DEL12F ! *** end of filename ddel12f.f95 *********************