! Last change: DOS 30 Jul 2000 1:52 pm ! *** copyright 2000 *** ! *** filename augrain2.f95 *** John F. Monahan ** ! ********************** ! ! SOME COMPILERS REQUIRE MODULE TO PREVIOUSLY EXIST OR BE ! AT THE BEGINNING MODULE passinfo real(KIND=8) sx,slx integer, parameter :: n = 45 ! number of years END MODULE passinfo ! program augrain2 ! augrain2 -- mle via regular optimizer ! for fitting two parameter gamma ! data are August rainfall data for Raleigh, NC ! declare functions passed as arguments ! pass data to functions using module use passinfo implicit none real(KIND=8), external :: ull real(KIND=8), dimension(2) :: v,typx,del real(KIND=8), dimension(3) :: h real(KIND=8) x,sxx,a0,b0,al,ar,f integer i,icode ! formats 20 format(' Fit two param gamma to rainfall data'/ & & ' get mle via plum1t') 21 format(2x,f5.2) 22 format(' sum, variance, sum of logs',3f12.6) 23 format(' initial estimates from method of moments, '/ & & ' alpha=',f12.6,' beta=',f12.6,' initial ll=',f12.6) 24 format(/' final estimates from plum1t, '/ & & ' alpha=',f12.6,' beta=',f12.6,' final ll=',f12.6) 25 format(' gradient vector ',2f12.6/ & & ' hessian matrix (aa, ab, bb)',3f12.6) 26 format(/' covariance matrix of estimates '/' alpha ',2e16.4/ & & ' beta ',2e16.4) 27 format(' std err(alpha)= ',e16.6,9x,'std err(beta)= ',e16.6) ! open( unit=6, file='augrain2.out') ! read August rainfall data -- 45 years open( unit=8, file='augrain.dat') ! write(6,20) ! initialize sums sx = 0.d0 sxx = 0.d0 slx = 0.d0 ! read data and compute sums do i = 1,n read(8,21) x sx = sx + x if( i .gt. 1 ) & & sxx = sxx + (real(i,8)*x-sx)**2 / real(i*(i-1),8) !KIND slx = slx + log(x) end do ! loop on i ! get sample variance sxx = sxx / real(n-1,8) !KIND write(6,22) sx,sxx,slx ! get starting values from method of moments b0 = real(n,8)*sxx / sx !KIND a0 = sx*sx / (real(n*n,8)*sxx) !KIND ! how good is prelim est v(1) = a0 v(2) = b0 f = ull(v) write(6,23) a0,b0,f ! max likelihood via plum1t typx(1) = v(1)*1.d-4 typx(2) = v(2)*1.d-4 call plum1t(ull,v,2,20,typx,f,del,h,6,icode) f = ull(v) write(6,24) v(1),v(2),f ! now get gradient and hessian typx(1) = v(1)*1.d-4 typx(2) = v(2)*1.d-4 call del12f(ull,v,2,typx,del,h) write(6,25) del,h ! get inverse x = h(1)*h(3) - h(2)*h(2) h(1) = h(1)/x h(3) = h(3)/x h(2) = - h(2)/x write(6,26) h(3),h(2),h(2),h(1) ! standard errors al = sqrt(h(1)) ar = sqrt(h(3)) write(6,27) ar,al stop end program augrain2 real(KIND=8) function ull(v) ! likelihood function for two parameter gamma use passinfo implicit none real(KIND=8), dimension(2), intent(in) :: v real(KIND=8) a,b,dlgama a = v(1) b = v(2) ull = - n*a*log(b) - n*dlgama(a) + (a-1.d0)*slx - sx/b ! unlikelihood to be minimized ull = - ull return end function ull REAL(KIND=8) FUNCTION DLGAMA(X) ! COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: X REAL(KIND=8), DIMENSION(7) :: P ! COEFS FOR RATIONAL FUNCTION REAL(KIND=8), DIMENSION(6) :: Q REAL(KIND=8), DIMENSION(3) :: PT ! COEFS FOR STIRLING CORRECTION REAL(KIND=8), DIMENSION(2) :: QT REAL(KIND=8) XN,A,V REAL(KIND=8), PARAMETER :: LSQR2PI = .918938533204672742D0 DATA P/ .378601050348257245d4, .207745979389418732d4, & & .893581804523749814d3, .222112396168011795d3, & & .489543462279099381d2, .61260674503360843d1, & & .778079585613300576d0 / DATA Q/ .378601050348257187d4, .476793860503687915d3, & & -.867230987531102994d3, .835500586679197696d2, & & .507884753288954097d2, -.134004147857813483d2 / DATA PT/ .288119283935546015d0, .498030766924499634d0, & & .691561607375687d-1 / DATA QT/ .345743140722674507d1, .609161691641660296d1 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0.D0 ) THEN DLGAMA = -1.D0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8.D0 ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = (PT(1)+(PT(2)+PT(3)/XN)/XN) / (QT(1)+(QT(2)+1.D0/XN)/XN) DLGAMA = (X-0.5D0)*LOG(X) - X + V/X + LSQR2PI ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.D0 XN = X DO WHILE( XN .LT. 2.D0 ) V = V/XN XN = XN + 1.D0 END DO ! WHILE( XN .LT. 2.D0 ) DO WHILE( XN .GT. 3.D0 ) XN = XN - 1.D0 V = V * XN END DO ! WHILE( XN .GT. 3.D0 ) A = XN-2.D0 ! RATIONAL FUNCTION APPROXIMATION V =(P(1)+A*(P(2)+A*(P(3)+A*(P(4)+A*(P(5)+A*(P(6)+A*P(7))))))) & & * V / (Q(1)+A*(Q(2)+A*(Q(3)+A*(Q(4)+A*(Q(5)+A*(Q(6)+A)))))) DLGAMA = LOG(V) ENDIF RETURN END FUNCTION DLGAMA subroutine plum1t(f,x,n,mit,typx,fnew,grad,hess, & & outcod,termcd) ! subroutine to minimize f ! requiring only function evaluations f(x) ! computes gradient and Hessian using numerical derivatives ! designed following recommendations of Dennis and Schnabel (1983) ! Numerical Method for Unconstrained Optimization and Nonlinear ! Equations, Prentice-Hall ! ! f function to be minimized (function subprogram) ! x real vector in: starting point for search ! out: minimum ! n integer dimension of x ! mit integer maximum number of iterations ! typx real vector typical values for |x|, used for scaling ! fnew real out: value of function at minimum ! grad real vector out: gradient of f at minimum ! hess real vector out: Hessian matrix at minimum in symmetric ! storage mode ** declare n*(n+1)/2 length ** ! outcod integer output code for writing messages ! termcd integer termination code ! termcd=1 means stopped on small gradient (success) ! termcd=2 means stop for small step size ! termcd=3 means stop for failing backtracks ! termcd=4 means stop -- too many iterations ! termcd=5 means too many big steps, maybe unbounded ! termcd=6 means Hessian not positive definite ! termcd=7 means descent direction goes up (not posdf) ! 1 is success, 2 or 3 maybe done, perhaps limited by accuracy of f ! 4, 5, 6, 7 are failures ! ! required external routines: ! del12f computes numerical gradient and Hessian ! chlzky computes Cholesky decomposition of pos def matrix ! chlzhi solves system of equations using Cholesky ! chlzih factor as matrix ! adjust called by chlzky to control size of determinant ! ! j f monahan (july, 1993) ! recoded October 1999, April 2000 for Fortran 95 ! IMPLICIT NONE REAL(KIND=8) f INTEGER, INTENT(IN) :: n,mit REAL(KIND=8), DIMENSION( n ), INTENT(IN OUT) :: x,typx,grad REAL(KIND=8), DIMENSION( (n*(n+1))/2 ), INTENT(OUT) :: hess REAL(KIND=8), INTENT(IN OUT) :: fnew INTEGER, INTENT(IN) :: outcod INTEGER, INTENT(OUT) :: termcd REAL(KIND=8), DIMENSION( n ) :: step,bstep REAL(KIND=8), DIMENSION( (n*(n+1))/2 ) :: thess REAL(KIND=8) stepmx,stsize,fold,flast,gradmx,det,c,slope,rlleng,& & lam2sm,sl,sltry,slast,v1,v2,a,b,disc ! integer i,nn,kit,idet,obig,mbig ! ! *** these constants should change with computer *** ! designed for ieee binary double precision ! machine epsilon REAL(KIND=8), PARAMETER :: eta = 4.44d-16 ! if the evaluations of function f are only accurate to d digits, ! then a better value would be eta = 10.**(-d) ! ! other constants are functions of this machine epsilon ! ! gradtl used for stopping on small gradient REAL(KIND=8), PARAMETER :: gradtl = 7.63d-6 !**d gradtl = cbrt(eta) ! steptl used for stopping on small steps REAL(KIND=8), PARAMETER :: steptl = 5.82d-11 !**d steptl = cbrt(eta)*cbrt(eta) ! ddiff used for finite differences REAL(KIND=8), PARAMETER :: ddiff = 7.63d-6 !**d ddiff = gradtl ! alphmn is required improvement in step REAL(KIND=8), PARAMETER :: alphmn = 1.d-3 ! *** done with constants *** ! 21 format(' begin iteration ',i4,' out of max ',i4,' current x'& & /2x,6f12.6/4x,6f12.6/4x,6f12.6/6x,6f12.6) 22 format(' stepsize',f10.4,' stepmx ',f10.4,' lambda ',f5.3, & & ' fnew ',f10.4) 23 format(' quadratic backtrack') 24 format(' cubic backtrack ') 43 format(' stopped -- too many iterations') 44 format(' stopped -- backtracking failed, but might be done') 45 format(' stopped -- too many big steps, probably unbounded') 46 format(' stopped -- taking steps too small') 47 format(' stopped -- Hessian no longer positive definite') 48 format(' stopped -- gradient small -- success') 49 format(' stopped -- descent direction going up ') ! ! need this later nn = (n*(n+1))/2 ! initialize iteration count kit = 0 ! find maximum stepsize stepmx = 0. stsize = 0. do i = 1,n typx(i) = abs(typx(i)) stepmx = stepmx + ( x(i)/typx(i) ) * ( x(i)/typx(i) ) stsize = stsize + 1./( typx(i) * typx(i) ) end do ! loop on i stepmx = 400.*sqrt( max( stepmx, stsize ) ) ! stepmx limits size of any step ! and keeps from running far away ! ! first evaluate function here fold = f(x) ! get gradients (and Hessian) ! but first get differences do i = 1,n step(i) = ddiff*max( abs(x(i)), typx(i) ) end do ! loop on i ! now call del12f call del12f(f,x,n,step,grad,hess) do i = 1,nn thess(i) = hess(i) ! save Hessian end do ! loop on i ! have gradient and Hessian ! test if gradient too small to continue gradmx = 0. do i = 1,n ! also copy gradient for next step step(i) = - grad(i) c = max( abs(x(i)), typx(i) )/abs( fold ) gradmx = max( gradmx, abs(grad(i))*c ) end do ! loop on i ! are we done already? if( gradmx .lt. gradtl/1000. ) then ! successful completion ! stop on gradient too small write(outcod,48) termcd = 1 return end if ! (gradmx.lt.gradtl/1000) ! divide by 1000 -- stricter at start than finish ! ! gradient not zero -- where to move? ! compute tenative step using Newton ! ! first Cholesky factorization of hess call chlzky(thess,n,det,idet) ! if not pos def then die quickly if( idet .lt. -2000000000 ) then ! Hessian not positive definite! write(outcod,47) termcd = 6 return end if ! ( Cholesky failed ) ! compute step call chlzhi(thess,n,step) call chlzih(thess,n,step) ! ! now the big section of line search code ! ! start main loop do kit = 1,mit ! write(outcod,21) kit,mit,(x(i),i=1,n) ! first step -- see how well Newton step does ! ! how big is this step stsize = 0. do i = 1,n stsize = stsize + ( step(i)/typx(i) ) * ( step(i)/typx(i) ) end do ! loop on i stsize = sqrt( stsize ) ! if step is too big, shorten it if( stsize .gt. stepmx ) then c = stepmx/stsize do i = 1,n step(i) = c * step(i) end do ! loop on i stsize = stepmx ! end if ! ( stsize .gt. stepmx ) ! compute slope of f in step direction slope = 0. rlleng = 0. do i = 1,n slope = slope + step(i)*grad(i) rlleng = max( rlleng, abs(step(i))/max(abs(x(i)),typx(i)) ) end do ! loop on i if( slope .gt. 0. ) then ! steepest descent seems to go up ! quadratic form in inverse of Hessian < 0 ! really step*grad > 0 write(outcod,49) termcd = 7 return end if ! ( slope .gt. 0. ) ! temporary branch if something is goofy ! lam2sm is smallest relative stepsize lam2sm = steptl / rlleng ! ! go to new point and test the water ! ! sl = 1 means full step sl = 1. ! make new point do ! *** backtrack repeat do *** do i = 1,n bstep(i) = x(i) + sl*step(i) end do ! loop on i ! evaluate function fnew = f(bstep) write(outcod,22) stsize,stepmx,sl,fnew ! how much improvement?? if( fnew .lt. fold + alphmn*slope*sl ) exit ! from backtracking ! if not an improvement, then backtrack ! ! have backtracked too much already? if( sl .lt. lam2sm ) then ! stop because of failing backtracks write(outcod,44) termcd = 3 return end if ! ( sl .lt. lam2sm ) ! backtracking didn't work -- quit ! ! backtrack code ! if( sl .eq. 1. ) then ! first fit quadratic write(outcod,23) sltry = - slope/( 2.*(fnew - fold - slope) ) flast = fnew slast = sl sl = max( sl/10., sltry ) ! go try this one else write(outcod,24) ! next try with cubic if failed before ! this code is messy v1 = fnew - fold - sl*slope v2 = flast - fold - slast*slope a = ( v1/(sl*sl) - v2/(slast*slast) )/(sl-slast) b = ( v2*sl/(slast*slast) - v1*slast/(sl*sl) )/(sl-slast) disc = b*b - 3.*a*slope if ( a .eq. 0. ) then ! cubic is a quadratic sltry = - slope/(2.*b) else ! solve cubic sltry = (sqrt(disc)-b)/(3.*a) end if ! ( a .eq. 0 ) ! shorten at least by half sltry = min( sltry, sl/2. ) flast = fnew slast = sl sl = max( sl/10., sltry ) ! go and try this one end if ! ( sl .eq. 1. ) end do ! *** backtrack loop **** ! ! last step was fruitful, where next? ! was the last step as big as could? obig = 0 if( (sl .eq. 1.) .and. (stsize .gt. .99*stepmx ) ) obig = 1 ! several things at once ! copy new point, new differences, stepsize stsize = 0. do i = 1,n step(i) = x(i) - bstep(i) x(i) = bstep(i) c = max( abs(x(i)), typx(i) ) stsize = max( stsize, abs(step(i))/c ) step(i) = ddiff*max( abs(x(i)), typx(i) ) end do ! loop on i ! ! get gradient and Hessian at new point call del12f(f,x,n,step,grad,hess) do i = 1,nn thess(i) = hess(i) ! save Hessian end do ! loop on i ! do I stop now? ! first see if gradient too small to continue gradmx = 0. do i = 1,n ! also copy gradient for next step step(i) = - grad(i) c = max( abs(x(i)), typx(i) )/abs(fold) gradmx = max( gradmx, abs(grad(i)*c) ) end do ! loop on i fold = fnew if ( gradmx .lt. gradtl ) then ! successful completion ! stop on gradient too small write(outcod,48) termcd = 1 return end if ! ( gradmx .lt. gradtl ) ! gradient not too small, what about step? if ( stsize .lt. steptl ) then ! stop because step size is too small write(outcod,46) termcd = 2 return end if ! ( stsize .lt. steptl ) ! too many maxsteps? are we running away? if ( obig .gt. 0 ) then ! things are ok -- work on next step mbig = mbig + 1 ! we've done mbig maxsteps in a row if ( mbig .eq. 5 ) then ! stop because of too many maxsteps ! probably unbounded minimum write(outcod,45) termcd = 5 return end if ! ( mbig .eq. 5 ) ! work on next step end if ! ( obig .gt. 0 ) if( obig .eq. 0 ) mbig = 0 ! reset count if last one was ok ! ! compute tentative step using Newton call chlzky(thess,n,det,idet) ! if not pos def die quickly if( idet .lt. -2000000000 ) then ! Hessian not positive definite! write(outcod,47) termcd = 6 return end if ! ( Cholesky failed ) ! compute step call chlzhi(thess,n,step) call chlzih(thess,n,step) ! now go back to the line search code ! too many iterations? end do ! loop on kit ! ! stop because of too many iterations write(outcod,43) termcd = 4 return end subroutine plum1t SUBROUTINE ADJUST(D,I) ! NORMALIZES D WHILE KEEPING CONSTANT VALUE OF D * (2**I) INTEGER, PARAMETER :: IBIG = -2147483644 REAL(KIND=8), INTENT (IN OUT) :: D INTEGER, INTENT(IN OUT) :: I IF ( D .GT. 0.0 ) THEN DO WHILE ( D .LT. 1.0 ) D = D*16. I = I-4 ENDDO DO WHILE ( D .GT. 16.0 ) D = D/16. I = I+4 ENDDO ELSE I = IBIG ! IF D < 0 THEN I = - BIG ENDIF RETURN END SUBROUTINE ADJUST SUBROUTINE CHLZKY(A,N,DET,IDET) ! CHLSKY COMPUTES THE CHOLESKY (SQUARE-ROOT) FACTORIZATION ! A = L * ( L - TRANSPOSE ) WHERE L IS LOWER TRIANGULAR ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! A IS OVERWRITTEN BY L IN ITS LOWER TRIANGULAR PART ! SUBROUTINE ADJUST KEEPS DETERMINANT FROM EXPLODING USING ! 1 LE DET LE 16 AND DETERMINANT OF A IS DET*2**IDET ! ! J F MONAHAN (DEC,1983) DEPT OF STAT, N C S U, RALEIGH, N C 27650 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL(KIND=8), INTENT(OUT) :: DET INTEGER, INTENT(OUT) :: IDET REAL(KIND=8) T,S INTEGER I,J,K,IM1,JM1 ! INITIALIZE FOR DETERMINANT DET = 1. IDET = 0 ! DO I-TH ROW DO I = 1,N T = A( (I*(I+1))/2 ) ! FIRST ROW IS TRIVIAL IF(I .GT. 1) THEN IM1 = I-1 DO J=1,IM1 ! WORK ON (I,J)-TH ELEMENT S = A( (I*(I-1))/2 + J ) IF(J .GT. 1) THEN JM1 = J-1 DO K = 1,JM1 S = S - A( (J*(J-1))/2 + K ) * A( (I*(I-1))/2 + K ) ENDDO ! LOOP ON K ENDIF A( (I*(I-1))/2 + J )= S / A( (J*(J+1))/2 ) ! WORK ON DIAGONAL ELEMENT T = T - A( (I*(I-1))/2 + J ) * A( (I*(I-1))/2 + J ) ENDDO ! LOOP ON J ! FINISHED WITH LOOP ON J ! UPDATE DETERMINANT WITH DIAGONAL ENDIF DET = DET * T ! KEEP DET FROM EXPLODING CALL ADJUST(DET,IDET) ! ABORT IF DET IS NEGATIVE OR ZERO IF(IDET.LT.-2000000000) RETURN ! FINISH A POSITIVE DIAGONAL A( (I*(I+1))/2 ) = SQRT(T) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZKY SUBROUTINE CHLZHI(A,N,Y) ! ! OVERWRITES Y WITH INVERSE( L ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE QUADRATIC FORM IN THE INVERSE OF A POSITIVE DEFINITE ! MATRIX A, THAT IS Y' INVERSE( A ) Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE THE SUM OF SQUARES OF THE ELEMENTS OF Y ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,J,IM1 ! ! SINCE THE MATRIX IS LOWER TRIANGULAR, SOLVE FROM THE TOP DOWN ! DO I = 1,N ! BRANCH AROUND ON THE FIRST ONE IF( I .GT. 1 ) THEN IM1 = I - 1 DO J = 1,IM1 Y(I) = Y(I) - A( (I*(I-1))/2 + J )*Y(J) ENDDO ENDIF Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON I RETURN END SUBROUTINE CHLZHI SUBROUTINE CHLZIH(A,N,Y) ! *** NOTICE TRANSPOSE *** ! OVERWRITES Y WITH INVERSE( L' ) * Y ! WHERE L IS A LOWER TRIANGULAR MATRIX (N BY N) STORED IN A ! *** USES SYMMETRIC STORAGE FORMAT *** ! ! TO COMPUTE THE PRODUCT OF THE INVERSE OF A POSITIVE DEFINITE MATRIX ! AND A VECTOR Y, THAT IS INVERSE( A ) * Y THEN ! FIRST: FACTOR A BY CHLSKY(A,N,DET,IDET) ! THIS COMPUTES LOWER TRIANGULAR L SO THAT A = L * L' ! SECOND: COMPUTE INVERSE(L) * Y BY CHLSHI(A,N,Y) ! THIRD: COMPUTE INVERSE(L') * ( INVERSE(L)*Y ) BY CHLSIH(A,N,Y) ! ! J F MONAHAN (JULY,1984) DEPT OF STAT, N C S U, RALEIGH, N C 27695 USA ! REVISED JULY 1999, APRIL 2000 FOR FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION( (N*(N+1))/2 ), INTENT(IN) :: A REAL(KIND=8), DIMENSION( N ), INTENT(IN OUT) :: Y INTEGER I,II,IP1,J ! ! SINCE THE MATRIX IS UPPER TRIANGULAR, SOLVE FROM THE BOTTOM UP ! ! DO THE LAST ONE SEPARATELY Y(N) = Y(N) / A( (N*(N+1))/2 ) ! IF N = 1 THEN QUIT ELSE START LOOPING IF( N .EQ. 1 ) RETURN DO II = 2,N ! I BELOW NOW GOES FROM N-1 DOWN TO 1 I = N + 1 - II IP1 = I + 1 DO J = IP1,N Y(I) = Y(I) - A( (J*(J-1))/2 + I )*Y(J) ENDDO ! LOOP ON J Y(I) = Y(I) / A( (I*(I+1))/2 ) ENDDO ! LOOP ON II,I RETURN END SUBROUTINE CHLZIH 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 augrain2.f95 *********************