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