Program g01nbfe

!     G01NBF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g01nbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: abserr, beta, eps, y0
      Integer                          :: i, ifail, j, l1, l2, lda, ldb, ldc,  &
                                          ldsig, lmax, lwk, n
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), b(:,:), c(:,:), ela(:),      &
                                          emu(:), rmom(:), sigma(:,:), wk(:)
!     .. Executable Statements ..
      Write (nout,*) 'G01NBF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) beta, y0
      Read (nin,*) n, l1, l2

      lda = n
      ldb = n
      ldc = n
      ldsig = n
      lwk = 3*n*n + (8+l2)*n
      Allocate (a(lda,n),b(ldb,n),c(ldc,n),ela(n),emu(n),sigma(ldsig,n),       &
        wk(lwk),rmom(l2-l1+1))

!     Compute A, EMU, and SIGMA for simple autoregression
      Do i = 1, n
        Do j = i, n
          a(j,i) = 0.0E0_nag_wp
          b(j,i) = 0.0E0_nag_wp
        End Do
      End Do
      Do i = 1, n - 1
        a(i+1,i) = 0.5E0_nag_wp
        b(i,i) = 1.0E0_nag_wp
      End Do
      emu(1) = y0*beta
      Do i = 1, n - 1
        emu(i+1) = beta*emu(i)
      End Do
      sigma(1,1) = 1.0E0_nag_wp
      Do i = 2, n
        sigma(i,i) = beta*beta*sigma(i-1,i-1) + 1.0E0_nag_wp
      End Do
      Do i = 1, n
        Do j = i + 1, n
          sigma(j,i) = beta*sigma(j-1,i)
        End Do
      End Do

!     Use default accuracy
      eps = 0.0E0_nag_wp

!     Compute moments
      ifail = -1
      Call g01nbf('Ratio','Mean',n,a,lda,b,ldb,c,ldc,ela,emu,sigma,ldsig,l1,   &
        l2,lmax,rmom,abserr,eps,wk,ifail)
      If (ifail/=0) Then
        If (ifail<6) Then
          Go To 100
        End If
      End If

!     Display results
      Write (nout,99999) ' N = ', n, ' BETA = ', beta, ' Y0 = ', y0
      Write (nout,*)
      Write (nout,*) '      Moments'
      Write (nout,*)
      j = 0
      Do i = l1, lmax
        j = j + 1
        Write (nout,99998) i, rmom(j)
      End Do

100   Continue

99999 Format (A,I3,2(A,F6.3))
99998 Format (I3,E12.3)
    End Program g01nbfe