Program g13befe

!     G13BEF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g13bef, nag_wp, x04abf, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: iset = 1, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: d, s
      Integer                          :: dp, i, ifail, imwa, inc, isttf, itc, &
                                          iwa, kef, kfc, kpriv, kzef, kzsp,    &
                                          ldcm, ldxxy, mx, nadv, ncd, nce,     &
                                          ncf, ncg, ndf, ndv, nis, nit, npara, &
                                          nser, nsttf, nxxy, qp, qx, smx
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: cm(:,:), para(:), res(:), sd(:),     &
                                          sttf(:), wa(:), xxy(:,:)
      Real (Kind=nag_wp)               :: zsp(4)
      Integer                          :: mr(7)
      Integer, Allocatable             :: mt(:,:), mwa(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max, sum
!     .. Executable Statements ..
      Write (nout,*) 'G13BEF Example Program Results'
      Write (nout,*)

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

!     Read in problem size
      Read (nin,*) kzef, kfc, nxxy, nser, kef, nit, kzsp, kpriv
      If (kzsp/=0) Then
        Read (nin,*) zsp
      End If

!     Number of input series
      nis = nser - 1

!     Set the advisory channel to NOUT for monitoring information
      If (kpriv/=0) Then
        nadv = nout
        Call x04abf(iset,nadv)
      End If

      Allocate (mt(4,nser))

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

!     Read in transfer function
      Do i = 1, nis
        Read (nin,*) mt(1:4,i)
      End Do

!     Calculate NPARA and various other quantities required
!     for calculate array sizes
      npara = 0
      ncg = 0
      qx = 0
      smx = 0
      ncf = nser
      inc = 1
      Do i = 1, nis
        npara = npara + mt(2,i) + mt(3,i)
        If (mt(4,i)>1) Then
          ncg = ncg + sum(mt(1:3,i))
          If (mt(4,i)==3) Then
            mx = max(mt(1,i)+mt(2,i),mt(3,i))
            qx = max(qx,mx)
            smx = smx + mx
          End If
        Else If (mt(4,i)==1 .And. kef==3) Then
          If (mt(3,i)>0) Then
            ncf = ncf + 1
          End If
          inc = inc + 1
        End If
      End Do
      npara = npara + mr(1) + mr(3) + mr(4) + mr(6) + nser

!     Calculate size of arrays
      isttf = mr(4)*mr(7) + mr(2) + mr(5)*mr(7) + mr(3) +                      &
        max(mr(1),mr(6)*mr(7)) + ncg
      ldxxy = nxxy
      ldcm = npara
      qp = mr(3) + mr(6)*mr(7)
      dp = mr(2) + mr(5)*mr(7)
      If (mr(3)>0 .And. kef>1) Then
        inc = inc + 1
      End If
      If (kfc>0 .And. kef==3) Then
        inc = inc + 1
      End If
      qx = qp
      ncd = npara + kfc + smx
      If (mr(1)>0) Then
        ncf = ncf + inc
      End If
      If (mr(3)>0) Then
        ncf = ncf + inc
      End If
      If (mr(4)>0) Then
        ncf = ncf + inc
      End If
      If (mr(6)>0) Then
        ncf = ncf + inc
      End If
      If (qx>0) Then
        ncf = ncf + 1
      End If
      If (kfc>0) Then
        ncf = ncf + 1
      End If
      ncd = ncd + qx
      nce = nxxy + dp + 6*qx
      iwa = 2*ncd**2 + nce*(ncf+4)
      iwa = 2*iwa
      imwa = 16*nser + 7*ncd + 3*npara + 3*kfc + 27

      Allocate (xxy(ldxxy,nser),para(npara),sd(npara),cm(ldcm,npara),          &
        res(nxxy),sttf(isttf),wa(iwa),mwa(imwa))

!     Read in rest of data
      Read (nin,*) para(1:npara)
      Read (nin,*)(xxy(i,1:nser),i=1,nxxy)

      ifail = -1
      Call g13bef(mr,nser,mt,para,npara,kfc,nxxy,xxy,ldxxy,kef,nit,kzsp,zsp,   &
        itc,sd,cm,ldcm,s,d,ndf,kzef,res,sttf,isttf,nsttf,wa,iwa,mwa,imwa,      &
        kpriv,ifail)
      If (ifail/=0) Then
        If (ifail/=8 .And. ifail/=9) Then
          Go To 100
        End If
      End If

!     Display results
      Write (nout,99999) 'The number of iterations carried out is', itc
      Write (nout,*)
      If (ifail/=8 .And. ifail/=9) Then
        Write (nout,*)                                                         &
          'Final values of the parameters and their standard deviations'
        Write (nout,*)
        Write (nout,*) '   I            PARA(I)                 SD'
        Write (nout,*)
        Write (nout,99998)(i,para(i),sd(i),i=1,npara)
        Write (nout,*)
        Flush (nout)
        ifail = 0
        Call x04caf('General',' ',npara,npara,cm,ldcm,                         &
          'The correlation matrix is',ifail)
      End If
      Write (nout,*)
      Write (nout,*) 'The residuals and the z and n values are'
      Write (nout,*)
      Write (nout,*) '   I         RES(I)          z(t)           n(t)'
      Write (nout,*)
      ndv = nxxy - mr(2) - mr(5)*mr(7)
      Do i = 1, nxxy
        If (i<=ndv) Then
          Write (nout,99997) i, res(i), xxy(i,1:nser)
        Else
          Write (nout,99996) i, xxy(i,1:nser)
        End If
      End Do
      If (mr(2)/=0 .Or. mr(5)/=0) Then
        Write (nout,*)
        Write (nout,*)                                                         &
          '** Note that the residuals relate to differenced values **'
      End If
      Write (nout,*)
      Write (nout,99995) 'The state set consists of', nsttf, ' values'
      Write (nout,*)
      Write (nout,99994) sttf(1:nsttf)
      Write (nout,*)
      Write (nout,99999) 'The number of degrees of freedom is', ndf

100   Continue

99999 Format (1X,A,I4)
99998 Format (1X,I4,2F20.6)
99997 Format (1X,I4,3F15.3)
99996 Format (1X,I4,F30.3,F15.3)
99995 Format (1X,A,I4,A)
99994 Format (1X,6F10.4)
    End Program g13befe