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