Program g13ebfe

!     G13EBF Example Program Text

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: ddot, dgemv, dpotrf, dsyrk, dtrsv, g13ebf,        &
                             nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: inc1 = 1, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: dev, tol
      Integer                          :: i, ifail, info, istep, l, ldm, ldq,  &
                                          lds, lwk, m, n, ncall, tdq
      Logical                          :: full, stq
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), ax(:), b(:,:), c(:,:),       &
                                          h(:,:), k(:,:), p(:,:), q(:,:),      &
                                          r(:,:), s(:,:), u(:,:), us(:,:),     &
                                          wk(:), x(:), y(:)
      Integer, Allocatable             :: iwk(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: log
!     .. Executable Statements ..
      Write (nout,*) 'G13EBF Example Program Results'
      Write (nout,*)

!     Skip heading in data file
      Read (nin,*)

!     Read in problem size
      Read (nin,*) n, m, l, stq

      lds = n
      If (.Not. stq) Then
        ldq = l
        tdq = l
      Else
        ldq = 1
        tdq = 1
      End If
      ldm = m
      lwk = (n+m)*(n+m+l)
      Allocate (a(lds,n),b(lds,l),q(ldq,tdq),c(ldm,n),r(ldm,m),s(lds,n), &
        k(lds,m),h(ldm,m),u(lds,n),iwk(m),wk(lwk),ax(n),y(m),x(n),p(lds,n), &
        us(lds,n))

!     Read in the state covariance matrix, S
      Read (nin,*)(s(i,1:n),i=1,n)
!     Read in flag indicating whether S is the full matrix, or its 
!     Cholesky decomposition.
      Read (nin,*) full
!     If required, perform Cholesky decomposition on S       
      If (full) Then
!       The NAG name equivalent of dpotrf is f07fdf
        Call dpotrf('L',n,s,lds,info)
        If (info>0) Then
          Write (nout,*) ' S not positive definite'
          Go To 100
        End If
      End If

!     Read in initial state vector
      Read (nin,*) ax(1:n)
!     Read in transition matrix, A
      Read (nin,*)(a(i,1:n),i=1,n)
!     Read in noise coefficient matrix, B
      Read (nin,*)(b(i,1:l),i=1,n)
!     Read in measurement coefficient matrix, C
      Read (nin,*)(c(i,1:n),i=1,m)

!     Read in measurement noise covariance matrix, R
      Read (nin,*)(r(i,1:m),i=1,m)
!     Read in flag indicating whether R is the full matrix, or its Cholesky
!     decomposition
      Read (nin,*) full
!     If required, perform Cholesky decomposition on R
      If (full) Then
!       The NAG name equivalent of dpotrf is f07fdf
        Call dpotrf('L',m,r,ldm,info)
        If (info>0) Then
          Write (nout,*) ' R not positive definite'
          Go To 100
        End If
      End If

!     Read in state noise matrix Q, if not assume to be identity matrix
      If (.Not. stq) Then
        Read (nin,*)(q(i,1:l),i=1,l)
!       Read in flag indicating whether Q is the full matrix, or 
!       its Cholesky decomposition
        Read (nin,*) full
!       Perform cholesky factorisation on Q, if full matrix is supplied
        If (full) Then
!         The NAG name equivalent of dpotrf is f07fdf
          Call dpotrf('L',l,q,ldq,info)
          If (info>0) Then
            Write (nout,*) ' Q not positive definite'
            Go To 100
          End If
        End If
      End If

!     Read in control parameters
      Read (nin,*) ncall, tol

!     Display titles
      Write (nout,*) '         Residuals'
      Write (nout,*)

!     Loop through data       
      dev = 0.0E0_nag_wp
      Do istep = 1, ncall

!       Read in observed values
        Read (nin,*) y(1:m)

        If (istep==1) Then
!         Make first call to G13EBF
          ifail = 0
          Call g13ebf('T',n,m,l,a,lds,b,stq,q,ldq,c,ldm,r,s,k,h,u,tol,iwk,wk, &
            ifail)

!         The NAG name equivalent of dgemv is f06paf
          Call dgemv('N',n,n,one,u,lds,ax,inc1,zero,x,inc1)

        Else
!         Make remaining calls to G13EBF
          ifail = 0
          Call g13ebf('H',n,m,l,a,lds,b,stq,q,ldq,c,ldm,r,s,k,h,u,tol,iwk,wk, &
            ifail)
        End If

!       Perform time and measurement update x <= Ax + K(y-Cx)
!       The NAG name equivalent of dgemv is f06paf
        Call dgemv('N',m,n,-one,c,ldm,x,inc1,one,y,inc1)
        Call dgemv('N',n,n,one,a,lds,x,inc1,zero,ax,inc1)
        Call dgemv('N',n,m,one,k,lds,y,inc1,one,ax,inc1)
        x(1:n) = ax(1:n)

!       Display the residuals
        Write (nout,99999) y(1:m)

!       Update loglikelihood
!       The NAG name equivalent of dtrsv is f06pjf
        Call dtrsv('L','N','N',m,h,ldm,y,inc1)
!       The NAG name equivalent of ddot is f06eaf
        dev = dev + ddot(m,y,1,y,1)
        Do i = 1, m
          dev = dev + 2.0_nag_wp*log(h(i,i))
        End Do
      End Do

!     Calculate back-transformed x <- U^T x
!     The NAG name equivalent of dgemv is f06paf
      Call dgemv('T',n,n,one,u,lds,ax,inc1,zero,x,inc1)

!     Compute back-transformed P from S
      Do i = 1, n
        Call dgemv('T',n-i+1,n,one,u(i,1),lds,s(i,i),inc1,zero,us(1,i),inc1)
      End Do
!     The NAG name equivalent of dsyrk is f06ypf
      Call dsyrk('L','N',n,n,one,us,lds,zero,p,lds)

!     Display final results
      Write (nout,*)
      Write (nout,*) ' Final X(I+1:I) '
      Write (nout,99999) x(1:n)
      Write (nout,*)
      Flush (nout)
      ifail = 0
      Call x04caf('Lower','N',n,n,p,lds,'Final Value of P',ifail)
      Write (nout,99998) ' Deviance = ', dev

100   Continue

99999 Format (6F12.4)
99998 Format (A,E13.4)
    End Program g13ebfe