Program g05xcfe
! G05XCF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
! .. 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(rcord,npaths,d*(ntimes+1-a),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(rcord,npaths,idim,z)
! .. Use Statements ..
Use nag_library, Only: g05skf
! .. Scalar Arguments ..
Integer, Intent (In) :: idim, npaths, rcord
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: z(npaths*idim)
! .. Local Scalars ..
Integer :: ifail
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: ztmp(:,:), ztmp2(:,:)
Integer :: seed(1)
Integer, Allocatable :: state(:)
! .. Intrinsic Procedures ..
Intrinsic :: reshape, transpose
! .. Executable Statements ..
! 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,size(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)
If (rcord==1) Then
Allocate (ztmp(npaths,idim),ztmp2(idim,npaths))
ztmp(1:npaths,1:idim) = reshape(z,(/npaths,idim/))
ztmp2(1:idim,1:npaths) = transpose(ztmp)
z(1:npaths*idim) = reshape(ztmp2,(/npaths*idim/))
End If
End Subroutine get_z
Subroutine initialize_prng(genid,subid,seed,lseed,state)
! .. Use Statements ..
Use nag_library, Only: g05kff
! .. Scalar Arguments ..
Integer, Intent (In) :: genid, lseed, subid
! .. Array Arguments ..
Integer, Intent (In) :: seed(lseed)
Integer, Allocatable, Intent (Out) :: state(:)
! .. Local Scalars ..
Integer :: ifail, lstate
! .. Executable Statements ..
! 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