    PROGRAM g13ebfe

!      G13EBF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

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

20     CONTINUE

99999  FORMAT (6F12.4)
99998  FORMAT (A,E13.4)
    END PROGRAM g13ebfe
