Program g05pjfe

!     G05PJF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. 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 Procedures ..
      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 initializer 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