    PROGRAM g13befe

!      G13BEF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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 Functions ..
       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 20
          END IF
       END IF

!      Display results
       WRITE (nout,99999) 'The number of iterations carried out is', itc
       WRITE (nout,*)
       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)
       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

20     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
