    PROGRAM g13eafe

!      G13EAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : daxpy, ddot, dgemv, dpotrf, dtrmv, dtrsv,       &
                               g13eaf, 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, is_const, read_matrix, stq
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), ax(:), b(:,:), c(:,:),       &
                                          h(:,:), k(:,:), p(:,:), q(:,:),      &
                                          r(:,:), s(:,:), wk(:), x(:), y(:),   &
                                          ymean(:)
       INTEGER, ALLOCATABLE            :: iwk(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          log
!      .. Executable Statements ..
       WRITE (nout,*) 'G13EAF Example Program Results'
       WRITE (nout,*)

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

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

       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),iwk(m),wk(lwk),x(n),ymean(m),y(m),ax(n),p(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 (full), perform the 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,*) x(1:n)

!      Read in mean of the series
       READ (nin,*) ymean(1:m)

!      Read in control parameter
       READ (nin,*) ncall, tol

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

!      Initialise variables
       dev = zero
       read_matrix = .TRUE.

!      Loop through data
       DO istep = 1, ncall
!         Read in the various matrices. If the series is constant
!         then this only happens at the first call
          IF (read_matrix) THEN

!            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 (full), perform the 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

!            If series is constant set flag to false
             read_matrix = .NOT. is_const
          END IF

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

!         Call G13EAF
          ifail = 0
          CALL g13eaf(n,m,l,a,lds,b,stq,q,ldq,c,ldm,r,s,k,h,tol,iwk,wk,ifail)

!         Subtract the mean y:= y-ymean
!         The NAG name equivalent of daxpy is f06ecf
          CALL daxpy(m,-one,ymean,inc1,y,inc1)

!         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,1)
          CALL dgemv('N',n,n,one,a,lds,x,inc1,zero,ax,1)
          CALL dgemv('N',n,m,one,k,lds,y,inc1,one,ax,1)
          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,1)
!         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

!      Compute P from S
!      The NAG name equivalent of dtrmv is f06pff
       DO i = 1, n
          p(1:i,i) = s(i,1:i)
          CALL dtrmv('L','N','N',i,s,lds,p(1,i),inc1)
          p(i,1:i-1) = p(1:i-1,i)
       END DO

!      Display final results
       WRITE (nout,*)
       WRITE (nout,*) ' Final X(I+1:I) '
       WRITE (nout,*)
       WRITE (nout,99999) x(1:n)
       WRITE (nout,*)
       FLUSH (nout)
       ifail = 0
       CALL x04caf('Lower','Non-Diag',n,n,p,lds,'Final Value of P',ifail)
       WRITE (nout,*)
       WRITE (nout,99998) ' Deviance = ', dev

20     CONTINUE

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