Program g01nafe

!     G01NAF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g01naf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: beta, con
      Integer                          :: i, ifail, j, l, lda, ldsig, lwk, n
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), emu(:), rkum(:), rmom(:),    &
                                          sigma(:,:), wk(:)
!     .. Executable Statements ..
      Write (nout,*) 'G01NAF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) beta, con
      Read (nin,*) n, l

      ldsig = n
      lda = n
      lwk = 3*n*(n+1)/2 + n
      Allocate (a(lda,n),emu(n),sigma(ldsig,n),rkum(l),rmom(l),wk(lwk))

!     Compute A, EMU, and SIGMA for simple autoregression
      Do i = 1, n
        Do j = i, n
          a(j,i) = 0.0E0_nag_wp
        End Do
      End Do
      Do i = 1, n - 1
        a(i+1,i) = 0.5E0_nag_wp
      End Do
      emu(1) = con*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

!     Compute cumulants
      ifail = 0
      Call g01naf('M','M',n,a,lda,emu,sigma,ldsig,l,rkum,rmom,wk,ifail)

!     Display results
      Write (nout,99999) ' N = ', n, ' BETA = ', beta, ' CON = ', con
      Write (nout,*)
      Write (nout,*) '    Cumulants       Moments'
      Write (nout,*)
      Write (nout,99998)(i,rkum(i),rmom(i),i=1,l)

99999 Format (A,I3,2(A,F6.3))
99998 Format (I3,E12.4,4X,E12.4)
    End Program g01nafe