Program g05xefe
! G05XEF Example Program Text
! Mark 27.2 Release. NAG Copyright 2021.
! .. Use Statements ..
Use nag_library, Only: g05xaf, g05xbf, 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 :: b(:,:), c(:,:), intime(:), rcomm(:), &
start(:), term(:), times(:), z(:,:)
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 (rcomm(12*(ntimes+1)))
ifail = 0
Call g05xaf(t0,tend,times,ntimes,rcomm,ifail)
! Get additional information required by the bridge generator
Call get_bridge_gen_data(npaths,rcord,d,start,a,term,c)
! Generate the Z values
Call get_z(rcord,npaths,d,a,ntimes,z,b)
! Leading dimensions for the various input arrays
ldz = size(z,1)
ldc = size(c,1)
ldb = size(b,1)
! Call the Brownian bridge generator routine
ifail = 0
Call g05xbf(npaths,rcord,d,start,a,term,z,ldz,c,ldc,b,ldb,rcomm,ifail)
! Display the results
Call display_results(rcord,ntimes,d,b)
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 + 1.71_nag_wp*real(i,kind=nag_wp)
End Do
tend = t0 + 1.71_nag_wp*real(ntimes+1,kind=nag_wp)
! We suppose the following 3 times are very important and should be
! constructed first. Note: these are indices into INTIME
nmove = 3
Allocate (move(nmove))
move(1:nmove) = (/3,5,4/)
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 = 3
a = 0
Allocate (start(d),term(d),c(d,d))
start(1:d) = 0.0_nag_wp
! As A = 0, TERM need not be initialized
! We want the following covariance matrix
c(:,1) = (/6.0_nag_wp,1.0_nag_wp,-0.2_nag_wp/)
c(:,2) = (/1.0_nag_wp,5.0_nag_wp,0.3_nag_wp/)
c(:,3) = (/-0.2_nag_wp,0.3_nag_wp,4.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 get_z(rcord,npaths,d,a,ntimes,z,b)
! .. Use Statements ..
Use nag_library, Only: g05yjf
! .. Scalar Arguments ..
Integer, Intent (In) :: a, d, npaths, ntimes, rcord
! .. Array Arguments ..
Real (Kind=nag_wp), Allocatable, Intent (Out) :: b(:,:), z(:,:)
! .. Local Scalars ..
Integer :: idim, ifail
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: std(:), tz(:,:), xmean(:)
Integer, Allocatable :: iref(:), state(:)
Integer :: seed(1)
! .. Intrinsic Procedures ..
Intrinsic :: transpose
! .. Executable Statements ..
idim = d*(ntimes+1-a)
! Allocate Z
If (rcord==1) Then
Allocate (z(idim,npaths))
Allocate (b(d*(ntimes+1),npaths))
Else
Allocate (z(npaths,idim))
Allocate (b(npaths,d*(ntimes+1)))
End If
! We now need to generate the input quasi-random points
! First initialize the base pseudorandom number generator
seed(1) = 1023401
Call initialize_prng(6,0,seed,size(seed),state)
! Scrambled quasi-random sequences preserve the good discrepancy
! properties of quasi-random sequences while counteracting the bias
! some applications experience when using quasi-random sequences.
! Initialize the scrambled quasi-random generator.
Call initialize_scrambled_qrng(1,2,idim,state,iref)
! Generate the quasi-random points from N(0,1)
Allocate (xmean(idim),std(idim))
xmean(1:idim) = 0.0_nag_wp
std(1:idim) = 1.0_nag_wp
If (rcord==1) Then
Allocate (tz(npaths,idim))
ifail = 0
Call g05yjf(xmean,std,npaths,tz,iref,ifail)
z(:,:) = transpose(tz)
Else
ifail = 0
Call g05yjf(xmean,std,npaths,z,iref,ifail)
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 initialize_scrambled_qrng(genid,stype,idim,state,iref)
! .. Use Statements ..
Use nag_library, Only: g05ynf
! .. Scalar Arguments ..
Integer, Intent (In) :: genid, idim, stype
! .. Array Arguments ..
Integer, Allocatable, Intent (Out) :: iref(:)
Integer, Intent (Inout) :: state(*)
! .. Local Scalars ..
Integer :: ifail, iskip, liref, nsdigits
! .. Executable Statements ..
liref = 32*idim + 7
iskip = 0
nsdigits = 32
Allocate (iref(liref))
ifail = 0
Call g05ynf(genid,stype,idim,iref,liref,iskip,nsdigits,state,ifail)
End Subroutine initialize_scrambled_qrng
Subroutine display_results(rcord,ntimes,d,b)
! .. Scalar Arguments ..
Integer, Intent (In) :: d, ntimes, rcord
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: b(:,:)
! .. Local Scalars ..
Integer :: i, j, k
! .. Executable Statements ..
Write (nout,*) 'G05XEF Example Program Results'
Write (nout,*)
Do i = 1, npaths
Write (nout,99999) 'Weiner Path ', i, ', ', ntimes + 1, &
' time steps, ', d, ' dimensions'
Write (nout,99997)(j,j=1,d)
k = 1
Do j = 1, ntimes + 1
If (rcord==1) Then
Write (nout,99998) j, b(k:k+d-1,i)
Else
Write (nout,99998) j, b(i,k:k+d-1)
End If
k = k + d
End Do
Write (nout,*)
End Do
99999 Format (1X,A,I0,A,I0,A,I0,A)
99998 Format (1X,I2,1X,20(1X,F10.4))
99997 Format (1X,3X,20(9X,I2))
End Subroutine display_results
End Program g05xefe