Program g05xbfe ! G05XBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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 and allocate B 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 + 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 non-free Wiener process npaths = 2 rcord = 2 d = 3 a = 1 Allocate (start(d),term(d),c(d,d)) start(1:d) = 0.0_nag_wp term(1:d) = (/1.0_nag_wp,0.5_nag_wp,0.0_nag_wp/) ! 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,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,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 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,*) 'G05XBF 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 g05xbfe