! Last change: DOS 2 Aug 2000 9:06 pm ! *** copyright 2000 *** ! *** filename chex122a.f95 *** John F. Monahan ** ! ********************** program chex122a ! chex122a -- integration of posterior from Chapter 10's Variance ! Components Example ! ! Importance Sampling based on posterior mode ! ! a -- using normal importance distribution ! implicit none ! pass function name as argument real(kind=8), external :: rlgpst ! integer, parameter :: np = 2500 ! replications integer, parameter :: p = 3 ! parameters real(kind=8), dimension(p) :: t,d1f,typx,eps,tmode real(kind=8), dimension((p*(p+1))/2) :: d2f real(kind=8), dimension(10) :: ss real(kind=8), dimension(np) :: wgt real(kind=8) rll,pstmx,w,zz,det,s1,s2 real(kind=8), parameter :: sr2pi = 2.50662827d0 ! sqrt(2*pi) real ran,gnroul integer i,j,ier,idet,k ! 23 format(f3.0,f7.3,f6.1,f4.0,f3.0,f3.0) 22 format(d16.9) 21 format(6f12.6) 24 format(i4,f12.6,4x,4f12.6/4x,5f12.6) 25 format(2i4,d20.6,f12.6) 26 format(/' Testing Weights, k =',i4,' beta-hat =',f12.6) 27 format(/' Asymptotic cov mx from posterior mode'/ & & 4x,f12.6/4x,2f12.6/4x,3f12.6) 28 format(/'Importance Sampling',i8,' points --Normalization const',& & e14.6/8x,'Posterior mean',10x,'Covariance Matrix') 29 format(' theta(1)',f12.6,6x,f12.6/ & & ' theta(2)',f12.6,6x,2f12.6/ & & ' theta(3)',f12.6,6x,3f12.6) ! ! initial values data tmode/ 1.625d0, 1.357d0, 7.337d0 / ! typical values data typx/ 2*1.d0, 2.d0 / ! epsilon values data eps/ 3*1.d-5 / ! output unit open( unit=6, file='chex122a.out' ) open( unit=8, file='chex122a.wgt' ) ! ! now just use starting value write(6,21) (tmode(j),j=1,3),pstmx pstmx = rlgpst(tmode) write(6,22) pstmx write(6,21) (tmode(j),j=1,3),pstmx ! call del12f(rlgpst,tmode,3,eps,d1f,d2f) write(6,21) (d1f(j),j=1,3) write(6,21) (d2f(j),j=1,6) ! call plum1t(rlgpst,tmode,3,25,typx,pstmx,d1f,d2f,6,ier) ! write(6,24) ier,pstmx write(6,21) tmode ! call del12f(rlgpst,tmode,3,eps,d1f,d2f) write(6,21) (d1f(j),j=1,3) write(6,21) (d2f(j),j=1,6) ! compute Cholesky factorization of the Hessian call chlzky(d2f,3,det,idet) write(6,24) idet,det ! ! have posterior mode and Hessian ! get ready for importance sampling ! ! set the rng zz = ran(5151917) ! initialize sums to zero ss = 0. ! two indicators s1 = 0. s2 = 0. ! ! main replication loop ! do i = 1,np ! generate normal(0,1)'s zz = 0. do j = 1,3 d1f(j) = gnroul(j+i) zz = zz + d1f(j)*d1f(j) end do ! loop on j ! rescale with inverse of Hessian call chlzih(d2f,3,d1f) ! relocate at mode do j = 1,3 t(j) = d1f(j) + tmode(j) end do ! loop on j ! actually getting computing -log(posterior) rll = rlgpst(t) - pstmx - zz/2.d0 ! importance sampling weights w = 0. ! avoid underflow if too small if( rll .lt. 30. ) w = exp(-rll) write(8,22) w wgt(i) = w ! get normalization constant and moments ss(1) = ss(1) + w ss(2) = ss(2) + w*t(1) ss(3) = ss(3) + w*t(2) ss(4) = ss(4) + w*t(3) ss(5) = ss(5) + w*t(1)*t(1) ss(6) = ss(6) + w*t(1)*t(2) ss(7) = ss(7) + w*t(2)*t(2) ss(8) = ss(8) + w*t(1)*t(3) ss(9) = ss(9) + w*t(2)*t(3) ss(10) = ss(10) + w*t(3)*t(3) s1 = dmax1( s1, w ) s2 = s2 + w*w end do ! loop on i (replications) ! posterior moments ss(2) = ss(2) / ss(1) ss(3) = ss(3) / ss(1) ss(4) = ss(4) / ss(1) ss(5) = ss(5) / ss(1) - ss(2)*ss(2) ss(6) = ss(6) / ss(1) - ss(2)*ss(3) ss(7) = ss(7) / ss(1) - ss(3)*ss(3) ss(8) = ss(8) / ss(1) - ss(2)*ss(4) ss(9) = ss(9) / ss(1) - ss(3)*ss(4) ss(10) = ss(10) / ss(1) - ss(4)*ss(4) ! get two indicators and print 'em s1 = s1 / ss(1) s2 = s2 / ( ss(1)*ss(1) ) w = 1. / s2 write(6,21) s1,s2,w ! fix normalization constant do j = 1,3 ss(1) = ss(1) * sr2pi / d2f( (j*(j+1))/2 ) end do ! loop on j ss(1) = ss(1) * exp(-pstmx) / float(np) ! write out results write(6,28) np,ss(1) write(6,29) ss(2),ss(5),ss(3),ss(6),ss(7),ss(4),ss(8),ss(9),ss(10) ! compute inverse of Hessian as covariance matrix call chlzoi(d2f,3) write(6,27) (d2f(j),j=1,6) ! analyze importance sampling weights call hsort(wgt,np) k = int( 4.*real(np)**.33 ) + 1 zz = 0. do i = 1,k zz = zz + log( wgt(np-i+1) ) end do ! loop on i zz = zz/real(k) - log( wgt(np-k) ) write(6,26) k,zz stop end program chex122a real(kind=8) function rlgpst(t) ! - log( posterior ) for exchange variance components problem ! data are 100*(rate-1.6) for USD/GBP rate in May 1997 ! ! t1 error/within variance component ! t2 trader/tmt variance component ! t3 mean parameter ! implicit none real(kind=8), dimension(3), intent(in) :: t integer, parameter :: p = 5 ! number of groups integer, dimension(p) :: ni real(kind=8), dimension(p) :: yb real(kind=8) w,a1,b1,a2,b2,b0,phi0,sf,sg ! integer i,nbig ! data for problem ! p number of groups ! ni number of obs in group i ! nbig sum of ni ! w within/error sum of squares ! a1,b1 prior for error component is inversegamma( a1, b1 ) ! a2,b2 prior for treatment component is inversegamma( a2, b2 ) ! b0,phi0 prior for mean in Normal( b0, phi0 ) ! yb mean for group i data ni/ 2, 1, 7, 3, 2 / data yb/ 6.8d0, 10.3d0, 5.9d0, 6.6d0, 8.5d0/ ! w = 13.10d0 ! within ss nbig = 15 a1 = 1. b1 = 1. a2 = 4. b2 = 4. b0 = 8. phi0 = 16. ! assign outlanding posterior (tiny) if out of range rlgpst = 1.d6 if( t(1) .le. 0.d0 ) return if( t(2) .le. 0.d0 ) return ! first part of posterior rlgpst = (a2 + p/2. + 1)*log(t(2)) + (a1+nbig/2.+1.)*log(t(1)) sg = 0. sf = 0. do i = 1,p sg = sg + log( ni(i)/t(1) + 1./t(2) ) sf = sf + (t(3)-yb(i))*(t(3)-yb(i))/(t(2)+t(1)/ni(i)) end do ! loop on i rlgpst = rlgpst + b2/t(2) + (b1+w/2.)/t(1) & & + (t(3)-b0)*(t(3)-b0)/(2.*phi0) & & + sg/2. + sf/2. ! ! using - log( posterior) *** reminder rlgpst = rlgpst return end function rlgpst 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 SUBROUTINE CHLZOI(A,N) ! ! OVERWRITES L WITH INVERSE(L')*INVERSE(L) WHERE L IS A LOWER ! TRIANGULAR (N BY N) MATRIX STORED IN A ! ! TO COMPUTE THE INVERSE OF A POSITIVE DEFINITE MATRIX A ! FIRST: FACTOR A = L * L' WHERE L IS LOWER TRIANGULAR ! USING CHLSKY(A,N,D,IDET) ! SECOND: COMPUTE THE INVERSE(A) = INVERSE(L') * INVERSE(L) ! USING CHLSOI(A,N) ! ! J F MONAHAN (JULY,1984) DEPT OF STATISTICS, N C S U, RALEIGH, N C 2769 ! 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) S INTEGER I,J,K,IM1,KP1 ! FIRST INVERT L IN PLACE DO K = 1,N ! SOLVE L * X = E(K) A( (K*(K+1))/2 ) = 1./A( (K*(K+1))/2 ) ! BRANCH AROUND IF K EQUALS N IF( K .LT. N ) THEN KP1 = K + 1 DO I = KP1,N S = 0. IM1 = I - 1 DO J = K,IM1 S = S - A( (I*(I-1))/2 + J )*A( (J*(J-1))/2 + K ) ENDDO ! LOOP ON J A( (I*(I-1))/2 + K ) = S/A( (I*(I+1))/2 ) ENDDO ! LOOP ON I ENDIF ENDDO ! LOOP ON K ! DONE WITH INVERSION, NOW COMPUTE INNER PRODUCT DO J = 1,N DO I = J,N ! GET (I,J) ELEMENT OF MATRIX IS THE INNER PRODUCT ! OF COLUMN I AND COLUMN J S = 0. DO K = I,N S = S + A( (K*(K-1))/2 + I )*A( (K*(K-1))/2 + J ) ENDDO ! LOOP ON K ! STORE THE INNER PRODUCT A( (I*(I-1))/2 + J ) = S ENDDO ! LOOP ON I ENDDO ! LOOP ON J RETURN END SUBROUTINE CHLZOI REAL FUNCTION GNROUL(IX) ! RATIO OF UNIFORMS ALGORITHM FOR GENERATING NORMAL(0,1) RV'S ! USING LEVA'S QUADRATIC INNER AND OUTER BOUNDS ! ! A J KINDERMAN AND J F MONAHAN (1977) "COMPUTER GENERATION OF RANDOM ! VARIABLES USING THE RATIO OF UNIFORM DEVIATES," ACM TRANSACTIONS ON ! MATHEMATICAL SOFTWARE, VOLUME 3, PP.257-260 ! ! J L LEVA (1992) "A FAST NORMAL RANDOM NUMBER GENERATOR," ACM ! TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 18, PP. 449-453. ! IMPLICIT NONE INTEGER, INTENT(IN) :: IX ! DUMMY ARGUMENT REAL, PARAMETER :: R = 1.7155277 ! FIRST CONSTANT R = SQRT(2/E) ! CONSTANTS S,T CENTER OF ELLIPSES REAL, PARAMETER :: S = .449871 REAL, PARAMETER :: T = -.386595 ! SHAPE PARAMETERS OF ELLIPSES REAL, PARAMETER :: A = .19600 REAL, PARAMETER :: B = .25472 ! RADII OF ELLIPSES REAL, PARAMETER :: R1 = .27597 REAL, PARAMETER :: R2 = .27846 REAL RAN REAL U,V,X,Y,Q ! START / RESTART HERE DO ! NOTE UNRESTRICTED DO ! GENERATE (U,V) UNIFORMLY IN BOX U = RAN(1) V = R*(RAN(2) - 0.5) X = U - S Y = ABS(V) - T ! COMPUTE QUADRATIC PIECE Q = X*X + Y*(A*Y - B*X) ! COMPUTE DEVIATE BEFORE TESTS GNROUL = V / U ! INNER BOUND -- QUICK ACCEPT CHECK IF( Q .LE. R1 ) RETURN ! OUTER BOUND -- QUICK REJECT CHECK IF( A .GT. R2 ) CYCLE ! FINAL RATIO OF UNIFORMS CHECK IF( GNROUL*GNROUL .LE. -4.*LOG(U) ) RETURN ! REJECT -- START OVER END DO ! END OF UNRESTRICTED DO END FUNCTION GNROUL SUBROUTINE HSORT(K,N) ! HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N ! ! ARGUMENTS ! K REAL VECTOR OF KEYS TO BE SORTED ! N NUMBER OF ITEMS TO BE SORTED ! ! TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT ! ! J F MONAHAN (DEC, 1999) FORTRAN 95 IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(KIND=8), DIMENSION(N), INTENT(IN OUT) :: K REAL(KIND=8) KK INTEGER I,L,NCUR ! ! DO NOTHING IF THERE'S NOTHING TO DO IF( N .LE. 1 ) RETURN ! INITIALIZE TO BUILDHEAP PART (LOOP ON L) L = N/2 + 1 NCUR = N DO I = L,1,-1 CALL HEAPIFY(I) END DO ! LOOP ON I DO I = 2,N ! SWITCH CURRENT LARGEST WITH BOTTOM KK = K(1) K(1) = K(NCUR) K(NCUR) = KK ! REHEAP WITH ONE SHORTER NCUR = NCUR - 1 CALL HEAPIFY(1) END DO ! LOOP ON NCUR RETURN CONTAINS SUBROUTINE HEAPIFY(II) INTEGER, INTENT(IN) :: II INTEGER I,J I = II DO J = 2*I ! IS IT A LEAF OR ARE THERE SONS? IF( J > NCUR ) EXIT ! A LEAF IF( J < NCUR ) THEN ! ANOTHER SON OF I IF( K(J+1) > K(J) ) J = J+1 ! LARGER SON IS K END IF ! A LEAF IF( K(J) > K(I) ) THEN ! EXCHANGE KK = K(J) K(J) = K(I) K(I) = KK I = J ELSE ! EXIT -- HEAP PROPERTY EXIT END IF END DO ! WHILE NOT A LEAF END SUBROUTINE HEAPIFY END SUBROUTINE HSORT REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN ! *** end of filename chex122a.f95 *********************