NAG Library Manual, Mark 29.2
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program g05xefe

!     G05XEF Example Program Text

!     Mark 29.2 Release. NAG Copyright 2023.

!     .. 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