!   G13AEF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module g13aefe_mod

!     G13AEF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: nin = 5, nout = 6
    Contains
      Subroutine piv(mr,par,npar,c,kfc,icount,s,g,h,ldh,igh,itc,zsp)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: c, s
        Integer, Intent (In)                 :: igh, itc, kfc, ldh, npar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: g(igh), h(ldh,igh), par(npar), &
                                                zsp(4)
        Integer, Intent (In)                 :: icount(6), mr(7)
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,99999) 'Iteration', itc, '  residual sum of squares = ', s
        Return

99999   Format (1X,A,I3,A,E11.4)
      End Subroutine piv
    End Module g13aefe_mod
    Program g13aefe

!     G13AEF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: g13aef, nag_wp, x04caf
      Use g13aefe_mod, Only: nin, nout, piv
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: c, s
      Integer                              :: iex, ifail, igh, ist, itc, iwa,  &
                                              kfc, kpiv, kzsp, ldh, ndf, nex,  &
                                              ngh, nit, npar, nst, nx
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: al(:), ex(:), exr(:), g(:),      &
                                              h(:,:), hc(:,:), par(:), sd(:),  &
                                              st(:), wa(:), x(:)
      Real (Kind=nag_wp)                   :: zsp(4)
      Integer                              :: icount(6), isf(4), mr(7)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: max
!     .. Executable Statements ..
      Write (nout,*) 'G13AEF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size etc
      Read (nin,*) nx, kfc, c

!     Read in the orders
      Read (nin,*) mr(1:7)

!     Calculate NPAR and various array lengths
      npar = mr(1) + mr(3) + mr(4) + mr(6)
      iex = mr(3) + (mr(6)*mr(7)) + nx
      igh = mr(3) + (mr(6)*mr(7)) + npar + kfc
      ist = (mr(4)*mr(7)) + mr(2) + (mr(5)*mr(7)) + mr(3) + &
        max(mr(1),(mr(6)*mr(7)))
      iwa = ((nx+1+mr(1)+(mr(4)*mr(7))+mr(3)+(mr(6)*mr(7)))*8) + (9*npar)

      ldh = igh + 1
      Allocate (x(nx),par(npar),ex(iex),exr(iex),al(iex),g(igh),sd(igh), &
        h(ldh,igh),st(ist),wa(iwa),hc(ldh,igh))

!     Read in data
      Read (nin,*) x(1:nx)

!     Read in initial values
      Read (nin,*) par(1:npar)

!     Read in control parameters
      Read (nin,*) kpiv, nit, kzsp
      If (kzsp==1) Then
        Read (nin,*) zsp(1:4)
      End If

!     Fit ARIMA model
      ifail = -1
      Call g13aef(mr,par,npar,c,kfc,x,nx,icount,ex,exr,al,iex,s,g,igh,sd,h, &
        ldh,st,ist,nst,piv,kpiv,nit,itc,zsp,kzsp,isf,wa,iwa,hc,ifail)
      If (ifail/=0) Then
        If (ifail<7) Then
          Go To 100
        End If
      End If

!     Display results
      nex = icount(4)
      ndf = icount(5)
      ngh = icount(6)
      If (ifail==0) Then
        Write (nout,99998) 'Convergence was achieved after', itc, ' cycles'
      Else
        Write (nout,99998) 'Iterative process ran for', itc, ' cycles'
      End If
      Write (nout,*)
      Write (nout,*) &
        'Final values of the PAR parameters and the constant are as follows'
      Write (nout,99997) par(1:npar), c
      Write (nout,*)
      Write (nout,99996) 'Residual sum of squares is', s, '  with', ndf, &
        ' degrees of freedom'
      Write (nout,*)
      Write (nout,*) 'The final values of ZSP were'
      Write (nout,99995) zsp(1:4)
      Write (nout,*)
      Write (nout,99999) 'The number of parameters estimated was', ngh
      Write (nout,*) '( backward forecasts, PAR and C, in that order )'
      Write (nout,*)
      Write (nout,*) 'The corresponding G array holds'
      Write (nout,99994) g(1:ngh)
      If ((ifail==0 .Or. ifail==9) .And. itc>0) Then
        Write (nout,*)
        Write (nout,*) 'The corresponding SD array holds'
        Write (nout,99994) sd(1:ngh)
        Write (nout,*)
        Flush (nout)
        ifail = 0
        Call x04caf('General',' ',ngh,ngh,h,ldh,'Corresponding H matrix', &
          ifail)
        Write (nout,*) 'Holds second derivatives in the upper half &
          &(including the main diagonal)'
        Write (nout,*) 'and correlation coefficients in the lower triangle'
      End If
      Write (nout,*)
      Write (nout,99993) 'EX, EXR, and AL each hold', nex, &
        ' values made up of', icount(1), ' back forecast(s),'
      Write (nout,99992) icount(2), ' differenced values, and'
      Write (nout,99992) icount(3), ' element(s) of reconstituted information'
      Write (nout,*)
      Write (nout,*) '  EX'
      Write (nout,99991) ex(1:nex)
      If (ifail==0 .Or. ifail==9) Then
        Write (nout,*)
        Write (nout,*) '  EXR'
        Write (nout,99991) exr(1:nex)
      End If
      If (ifail==0) Then
        Write (nout,*)
        Write (nout,*) '  AL'
        Write (nout,99991) al(1:nex)
      End If
      If (ifail==0 .Or. ifail==9) Then
        Write (nout,*)
        Write (nout,99998) 'The state set consists of', nst, ' values'
        Write (nout,99991) st(1:nst)
      End If

100   Continue

99999 Format (1X,A,I5)
99998 Format (1X,A,I3,A)
99997 Format (1X,4F10.4)
99996 Format (1X,A,F10.3,A,I4,A)
99995 Format (1X,4E15.4)
99994 Format (1X,10F9.4)
99993 Format (1X,A,I5,A,I5,A)
99992 Format (1X,I5,A)
99991 Format (1X,5F11.4)
    End Program g13aefe
