Example description
!   G13AEF Example Program Text
!   Mark 27.0 Release. NAG Copyright 2019.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: piv
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine piv(mr,par,npar,c,kfc,icount,s,g,h,ldh,igh,itc,zsp)

!       .. Parameters ..
        Integer, Parameter             :: param_inds(4) = (/1,3,4,6/)
!       .. 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)
!       .. Local Scalars ..
        Integer                        :: i, param, start_ind
!       .. Executable Statements ..
        Write (nout,*)

        If (itc==0) Then
          Write (nout,*) 'Begin monitoring output'
          Write (nout,99999) 'Degrees of freedom: ', icount(5)
        End If

        Write (nout,99999) 'Iteration: ', itc

param_loop: Do i = 1, 4
          param = mr(param_inds(i))
          If (param<=0) Then
            Cycle param_loop
          End If
          Select Case (param_inds(i))
          Case (1)
            Write (nout,*) 'Autoregressive parameters:'
            start_ind = 1
          Case (3)
            Write (nout,*) 'Moving-average parameters:'
          Case (4)
            Write (nout,*) 'Seasonal autoregressive parameters:'
          Case (6)
            Write (nout,*) 'Seasonal moving-average parameters:'
          End Select
          Write (nout,99998) par(start_ind:start_ind+param-1)
          start_ind = start_ind + param
        End Do param_loop

        If (kfc==1) Then
          Write (nout,*) 'Constant term:'
          Write (nout,99998) c
        End If

        Write (nout,*) 'Residual sum of squares:'
        Write (nout,99998) s

        Write (nout,*)                                                         &
          'Parameters (alpha, beta, delta, gamma) controlling the search ',    &
          'procedure:'
        Write (nout,99998) zsp

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

!     G13AEF Example Main Program

!     .. Use Statements ..
      Use g13aefe_mod, Only: nin, nout, piv
      Use nag_library, Only: g13aef, nag_wp, x04caf
!     .. 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