    PROGRAM g05pjfe

!      G05PJF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g05kff, g05pjf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: lseed = 1, nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                         :: genid, i, ifail, ii, ip, iq, j, k,   &
                                          k2, l, ldvar, ldx, lphi, lr, lstate, &
                                          ltheta, mode, n, nreal, rn, subid
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: phi(:), r(:), theta(:), var(:,:),    &
                                          x(:,:), xmean(:)
       INTEGER                         :: seed(lseed)
       INTEGER, ALLOCATABLE            :: state(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          max
!      .. Executable Statements ..
       WRITE (nout,*) 'G05PJF Example Program Results'
       WRITE (nout,*)

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

!      Read in the base generator information and seed
       READ (nin,*) genid, subid, seed(1)

!      Initial call to initialiser 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)

!      Read in the sample size and number of realizations
       READ (nin,*) n, nreal

!      Read in the number of coefficients
       READ (nin,*) k, ip, iq

       k2 = k**2
       rn = max(ip,iq)
       l = k*(k+1)/2
       IF (ip>0) THEN
          l = l + (ip-1)*k2
       END IF
       IF (k>=6) THEN
          lr = (5*rn**2+1)*k2 + (4*rn+3) + 4
       ELSE
          lr = ((ip+iq)**2+1)*k2 + (4*(ip+iq)+3)*k + &
             max(k*rn*(k*rn+2),k2*(ip+iq)**2+l*(l+3)+k2*(iq+1)) + 4
       END IF
       lphi = ip*k*k
       ltheta = iq*k*k
       ldvar = k
       ldx = k
       ALLOCATE (phi(lphi),theta(ltheta),var(ldvar,k),r(lr),x(ldx,n),xmean(k))

!      Read in the AR parameters
       DO l = 1, ip
          DO i = 1, k
             ii = (l-1)*k*k + i
             READ (nin,*) (phi(ii+k*(j-1)),j=1,k)
          END DO
       END DO

!      Read in the MA parameters
       DO l = 1, iq
          DO i = 1, k
             ii = (l-1)*k*k + i
             READ (nin,*) (theta(ii+k*(j-1)),j=1,k)
          END DO
       END DO

!      Read in the means
       READ (nin,*) xmean(1:k)

!      Read in the variance / covariance matrix
       READ (nin,*) (var(i,1:i),i=1,k)

!      For the first realization we need to set up the reference vector
!      as well as generate the series
       mode = 2

!      Generate NREAL realizations
D_LP:  DO rn = 1, nreal

          ifail = 0
          CALL g05pjf(mode,n,k,xmean,ip,phi,iq,theta,var,ldvar,r,lr,state,x, &
             ldx,ifail)

!         Display the results
          WRITE (nout,99999) ' Realization Number ', rn
          DO i = 1, k
             WRITE (nout,*)
             WRITE (nout,99999) '  Series number ', i
             WRITE (nout,*) '  -------------'
             WRITE (nout,*)
             WRITE (nout,99998) x(i,1:n)
          END DO
          WRITE (nout,*)

!         For subsequent realizations we use previous reference vector
          mode = 3

       END DO D_LP

99999  FORMAT (1X,A,I0)
99998  FORMAT (8(2X,F8.3))
    END PROGRAM g05pjfe
