    PROGRAM g08rafe

!      G08RAF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g08raf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: tol
       INTEGER                         :: i, idist, ifail, ip, j, ldprvr, ldx, &
                                          lparest, lvapvec, lwork, nmax, ns,   &
                                          nsum
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: eta(:), parest(:), prvr(:,:),        &
                                          vapvec(:), work(:), x(:,:), y(:),    &
                                          zin(:)
       INTEGER, ALLOCATABLE            :: irank(:), iwa(:), nv(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          maxval, sum
!      .. Executable Statements ..
       WRITE (nout,*) 'G08RAF Example Program Results'
       WRITE (nout,*)

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

!      Read number of samples, number of parameters to be fitted,
!      error distribution parameter and tolerance criterion for ties.
       READ (nin,*) ns, ip, idist, tol

       ALLOCATE (nv(ns))

!      Read the number of observations in each sample.
       READ (nin,*) nv(1:ns)

!      Calculate NSUM, NMAX and various array lengths
       nsum = sum(nv(1:ns))
       nmax = maxval(nv(1:ns))
       ldx = nsum
       ldprvr = ip + 1
       lvapvec = nmax*(nmax+1)/2
       lparest = 4*ip + 1
       lwork = nmax*(ip+1)
       ALLOCATE (y(nsum),x(ldx,ip),prvr(ldprvr,ip),irank(nmax),zin(nmax), &
          eta(nmax),vapvec(lvapvec),parest(lparest),work(lwork),iwa(nmax))

!      Read in observations and design matrices for each sample.
       READ (nin,*) (y(i),x(i,1:ip),i=1,nsum)

!      Display input information
       WRITE (nout,99999) 'Number of samples =', ns
       WRITE (nout,99999) 'Number of parameters fitted =', ip
       WRITE (nout,99999) 'Distribution =', idist
       WRITE (nout,99998) 'Tolerance for ties =', tol

       ifail = 0
       CALL g08raf(ns,nv,nsum,y,ip,x,ldx,idist,nmax,tol,prvr,ldprvr,irank,zin, &
          eta,vapvec,parest,work,lwork,iwa,ifail)

!      Display results
       WRITE (nout,*)
       WRITE (nout,*) 'Score statistic'
       WRITE (nout,99997) parest(1:ip)
       WRITE (nout,*)
       WRITE (nout,*) 'Covariance matrix of score statistic'
       DO j = 1, ip
          WRITE (nout,99997) prvr(1:j,j)
       END DO
       WRITE (nout,*)
       WRITE (nout,*) 'Parameter estimates'
       WRITE (nout,99997) parest((ip+1):(2*ip))
       WRITE (nout,*)
       WRITE (nout,*) 'Covariance matrix of parameter estimates'
       DO i = 1, ip
          WRITE (nout,99997) prvr(i+1,1:i)
       END DO
       WRITE (nout,*)
       WRITE (nout,99996) 'Chi-squared statistic =', parest(2*ip+1), ' with', &
          ip, ' d.f.'
       WRITE (nout,*)
       WRITE (nout,*) 'Standard errors of estimates and'
       WRITE (nout,*) 'approximate z-statistics'
       WRITE (nout,99995) (parest(2*ip+1+i),parest(3*ip+1+i),i=1,ip)

99999  FORMAT (1X,A,I2)
99998  FORMAT (1X,A,F8.5)
99997  FORMAT (1X,2F9.3)
99996  FORMAT (1X,A,F9.3,A,I2,A)
99995  FORMAT (1X,F9.3,F14.3)
    END PROGRAM g08rafe
