Program g05xcfe ! G05XCF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: g05xaf, g05xbf, g05xcf, g05xdf, g05xef, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nout = 6 ! .. Local Scalars .. Real (Kind=nag_wp) :: t0, tend Integer :: a, bgord, d, ifail, ldb, ldc, & ldz, nmove, npaths, ntimes, rcord ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: bb(:,:), bd(:,:), c(:,:), & diff(:), intime(:), rcommb(:), & rcommd(:), start(:), term(:), & times(:), zb(:,:), zd(:,:) Integer, Allocatable :: move(:) ! .. Intrinsic Procedures .. Intrinsic :: size ! .. Executable Statements .. ! Get information required to set up the bridge Call get_bridge_init_data(bgord,t0,tend,ntimes,intime,nmove,move) ! Make the bridge construction bgord Allocate (times(ntimes)) ifail = 0 Call g05xef(bgord,t0,tend,ntimes,intime,nmove,move,times,ifail) ! Initialize the Brownian bridge generator Allocate (rcommb(12*(ntimes+1)),rcommd(12*(ntimes+1))) ifail = 0 Call g05xaf(t0,tend,times,ntimes,rcommb,ifail) ifail = 0 Call g05xcf(t0,tend,times,ntimes,rcommd,ifail) ! Get additional information required by the bridge generator Call get_bridge_gen_data(npaths,rcord,d,start,a,term,c) Allocate (diff(d)) diff(:) = term(:) - start(:) Call allocate_arrays(rcord,npaths,d,a,ntimes,zb,bb,zd,bd) ! Generate the Z values Call get_z(npaths,d,a,ntimes,zb) zd(:,:) = zb(:,:) ! Leading dimensions for the various input arrays ldz = size(zb,1) ldc = size(c,1) ldb = size(bb,1) ! Call the Brownian bridge generator routine ifail = 0 Call g05xbf(npaths,rcord,d,start,a,term,zb,ldz,c,ldc,bb,ldb,rcommb, & ifail) ifail = 0 Call g05xdf(npaths,rcord,d,a,diff,zd,ldz,c,ldc,bd,ldb,rcommd,ifail) ! Display the results Call display_results(rcord,t0,tend,ntimes,intime,d,start,bb,bd) Contains Subroutine get_bridge_init_data(bgord,t0,tend,ntimes,intime,nmove,move) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (Out) :: t0, tend Integer, Intent (Out) :: bgord, nmove, ntimes ! .. Array Arguments .. Real (Kind=nag_wp), Allocatable, Intent (Out) :: intime(:) Integer, Allocatable, Intent (Out) :: move(:) ! .. Local Scalars .. Integer :: i ! .. Intrinsic Procedures .. Intrinsic :: real ! .. Executable Statements .. ! Set the basic parameters for a Wiener process ntimes = 10 t0 = 0.0_nag_wp Allocate (intime(ntimes)) ! We want to generate the Wiener process at these time points Do i = 1, ntimes intime(i) = t0 + real(i,kind=nag_wp) End Do tend = t0 + real(ntimes+1,kind=nag_wp) nmove = 0 Allocate (move(nmove)) bgord = 3 End Subroutine get_bridge_init_data Subroutine get_bridge_gen_data(npaths,rcord,d,start,a,term,c) ! .. Use Statements .. Use nag_library, Only: dpotrf ! .. Scalar Arguments .. Integer, Intent (Out) :: a, d, npaths, rcord ! .. Array Arguments .. Real (Kind=nag_wp), Allocatable, Intent (Out) :: c(:,:), start(:), & term(:) ! .. Local Scalars .. Integer :: info ! .. Executable Statements .. ! Set the basic parameters for a free Wiener process npaths = 2 rcord = 2 d = 2 a = 0 Allocate (start(d),term(d),c(d,d)) start(1:d) = (/0.0_nag_wp,2.0_nag_wp/) term(1:d) = (/1.0_nag_wp,0.0_nag_wp/) ! We want the following covariance matrix c(:,1) = (/6.0_nag_wp,-1.0_nag_wp/) c(:,2) = (/-1.0_nag_wp,5.0_nag_wp/) ! G05XBF works with the Cholesky factorization of the covariance matrix C ! so perform the decomposition Call dpotrf('Lower',d,c,d,info) If (info/=0) Then Write (nout,*) & 'Specified covariance matrix is not positive definite: info=', & info Stop End If End Subroutine get_bridge_gen_data Subroutine allocate_arrays(rcord,npaths,d,a,ntimes,zb,bb,zd,bd) ! .. Scalar Arguments .. Integer, Intent (In) :: a, d, npaths, ntimes, rcord ! .. Array Arguments .. Real (Kind=nag_wp), Allocatable, Intent (Out) :: bb(:,:), bd(:,:), & zb(:,:), zd(:,:) ! .. Local Scalars .. Integer :: idim ! .. Executable Statements .. idim = d*(ntimes+1-a) If (rcord==1) Then Allocate (zb(idim,npaths),zd(idim,npaths)) Allocate (bb(d*(ntimes+1),npaths),bd(d*(ntimes+1),npaths)) Else Allocate (zb(npaths,idim),zd(npaths,idim)) Allocate (bb(npaths,d*(ntimes+1)),bd(npaths,d*(ntimes+1))) End If End Subroutine allocate_arrays Subroutine get_z(npaths,d,a,ntimes,z) ! .. Use Statements .. Use nag_library, Only: g05skf ! .. Scalar Arguments .. Integer, Intent (In) :: a, d, npaths, ntimes ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: z(:,:) ! .. Local Scalars .. Integer :: idim, ifail ! .. Local Arrays .. Integer :: seed(1) Integer, Allocatable :: state(:) ! .. Executable Statements .. idim = d*(ntimes+1-a) ! We now need to generate the input pseudorandom points ! First initialize the base pseudorandom number generator seed(1) = 1023401 Call initialize_prng(6,0,seed,state) ! Generate the pseudorandom points from N(0,1) ifail = 0 Call g05skf(idim*npaths,0.0_nag_wp,1.0_nag_wp,state,z,ifail) End Subroutine get_z Subroutine initialize_prng(genid,subid,seed,state) ! .. Use Statements .. Use nag_library, Only: g05kff ! .. Scalar Arguments .. Integer, Intent (In) :: genid, subid ! .. Array Arguments .. Integer, Intent (In) :: seed(:) Integer, Allocatable, Intent (Out) :: state(:) ! .. Local Scalars .. Integer :: ifail, lseed, lstate ! .. Executable Statements .. lseed = size(seed,1) ! Initial call to initializer to get size of STATE array lstate = 0 Allocate (state(lstate)) ifail = 0 Call g05kff(genid,subid,seed,lseed,state,lstate,ifail) ! Reallocate STATE Deallocate (state) Allocate (state(lstate)) ! Initialize the generator to a repeatable sequence ifail = 0 Call g05kff(genid,subid,seed,lseed,state,lstate,ifail) End Subroutine initialize_prng Subroutine display_results(rcord,t0,tend,ntimes,intime,d,start,bb,bd) ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: t0, tend Integer, Intent (In) :: d, ntimes, rcord ! .. Array Arguments .. Real (Kind=nag_wp), Intent (In) :: bb(:,:), bd(:,:), intime(:), & start(:) ! .. Local Scalars .. Integer :: i, n ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: cum(:), unscaled(:) ! .. Executable Statements .. Allocate (cum(d),unscaled(d)) Write (nout,*) 'G05XCF Example Program Results' Write (nout,*) Do n = 1, npaths Write (nout,99999) 'Weiner Path ', n, ', ', ntimes + 1, & ' time steps, ', d, ' dimensions' Write (nout,'(A)') & ' Output of G05XBF Output of G05XDF Sum of G05XDF' cum(:) = start(:) If (rcord==1) Then unscaled(:) = bd(1:d,n)*(intime(1)-t0) cum(:) = cum(:) + unscaled(:) Write (nout,99998) 1, bb(1:d,n), bd(1:d,n), cum(1:d) Else unscaled(:) = bd(n,1:d)*(intime(1)-t0) cum(:) = cum(:) + unscaled(:) Write (nout,99998) 1, bb(n,1:d), bd(n,1:d), cum(1:d) End If Do i = 2, ntimes If (rcord==1) Then unscaled(:) = bd(1+(i-1)*d:i*d,n)*(intime(i)-intime(i-1)) cum(:) = cum(:) + unscaled(:) Write (nout,99998) i, bb(1+(i-1)*d:i*d,n), bd(1+(i-1)*d:i*d,n), & cum(1:d) Else unscaled(:) = bd(n,1+(i-1)*d:i*d)*(intime(i)-intime(i-1)) cum(:) = cum(:) + unscaled(:) Write (nout,99998) i, bb(n,1+(i-1)*d:i*d), bd(n,1+(i-1)*d:i*d), & cum(1:d) End If End Do i = ntimes + 1 If (rcord==1) Then unscaled(:) = bd(1+(i-1)*d:i*d,n)*(tend-intime(i-1)) cum(:) = cum(:) + unscaled(:) Write (nout,99998) i, bb(1+(i-1)*d:i*d,n), bd(1+(i-1)*d:i*d,n), & cum(1:d) Else unscaled(:) = bd(n,1+(i-1)*d:i*d)*(tend-intime(i-1)) cum(:) = cum(:) + unscaled(:) Write (nout,99998) i, bb(n,1+(i-1)*d:i*d), bd(n,1+(i-1)*d:i*d), & cum(1:d) End If Write (nout,*) End Do 99999 Format (1X,A,I0,A,I0,A,I0,A) 99998 Format (1X,I2,1X,2(1X,F8.4),2X,2(1X,F8.4),2X,2(1X,F8.4)) End Subroutine display_results End Program g05xcfe