Program g02mbfe

!     G02MBF Example Program Text

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g02buf, g02mbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: sw
      Integer                          :: i, ifail, intcpt, ip, k, ldb, lddtd, &
                                          lisx, lropt, m, mnstep, mtype, n,    &
                                          nstep, pm, pm2, pred
      Character (10)                   :: mean
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:,:), dtd(:,:), dy(:,:),           &
                                          fitsum(:,:), ropt(:), wmean(:)
      Real (Kind=nag_wp)               :: wt(1)
      Integer, Allocatable             :: isx(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: count, floor, max, repeat
!     .. Executable Statements ..

!     .. Executable Statements ..
      Write (nout,*) 'G02MBF Example Program Results'
      Write (nout,*)

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

!     Read in the problem size
      Read (nin,*) n, m

!     Read in the model specification
      Read (nin,*) mtype, pred, intcpt, mnstep, lisx

!     Read in variable inclusion flags (if specified) and calculate IP
      Allocate (isx(lisx))
      If (lisx==m) Then
        Read (nin,*) isx(1:lisx)
        ip = count(isx(1:m)==1)
      Else
        ip = m
      End If

!     Optional arguments (using defaults)
      lropt = 0
      Allocate (ropt(lropt))

!     Read in the augmented matrix [D y] and calculate cross-product matrices
!     (NB: Datasets with a large number of observations can be split into
!     blocks with the resulting cross-product matrices being combined
!     using G02BZF)
      Allocate (dy(n,m+1))
      Read (nin,*)(dy(i,1:m),dy(i,m+1),i=1,n)

      pm = m*(m+1)/2
      pm2 = (m+1)*(m+2)/2

!     We are calculating the cross-product matrix using G02BUF
!     which returns it in packed storage
      lddtd = 1

!     Calculate the cross-product matrices
      Allocate (wmean(m+1),dtd(1,pm2))
      If (intcpt==1) Then
        mean = 'Mean'
      Else
        mean = 'Zero'
      End If
      ifail = 0
      Call g02buf(mean,'UnWeighted',n,m+1,dy,n,wt,sw,wmean,dtd(1,:),ifail)
!     The first PM elements of DTD(1,:) contain the cross-products of D
!     elements DTD(1,PM+1:PM2-1) contains cross-product of D with y and
!     DTD(1,PM2) contains cross-product of y with itself

!     Allocate output arrays
      ldb = ip
      Allocate (b(ldb,mnstep+1),fitsum(6,mnstep+1))

!     Call the model fitting routine
      ifail = -1
      Call g02mbf(mtype,pred,intcpt,n,m,dtd,lddtd,isx,lisx,dtd(1,pm+1:pm2-1),  &
        dtd(1,pm2),mnstep,ip,nstep,b,ldb,fitsum,ropt,lropt,ifail)
      If (ifail/=0) Then
        If (ifail/=112 .And. ifail/=161 .And. ifail/=162 .And. ifail/=163)     &
          Then
!         IFAIL = 112, 161, 162 and 163 are warnings, so no need to terminate
!         if they occur
          Go To 100
        End If
      End If

!     Display the parameter estimates
      Write (nout,*) ' Step ', repeat(' ',max((ip-2),0)*5),                    &
        ' Parameter Estimate'
      Write (nout,*) repeat('-',5+ip*10)
      Do k = 1, nstep
        Write (nout,99998) k, b(1:ip,k)
      End Do
      Write (nout,*)
      Write (nout,99999) 'alpha: ', wmean(m+1)
      Write (nout,*)
      Write (nout,*)                                                           &
        ' Step     Sum      RSS       df       Cp       Ck     Step Size'
      Write (nout,*) repeat('-',64)
      Do k = 1, nstep
        Write (nout,99997) k, fitsum(1:2,k), floor(fitsum(3,k)+0.5_nag_wp),    &
          fitsum(4:6,k)
      End Do
      Write (nout,*)
      Write (nout,99999) 'sigma^2: ', fitsum(5,nstep+1)

100   Continue
99999 Format (1X,A,F9.3)
99998 Format (2X,I3,10(1X,F9.3))
99997 Format (2X,I3,2(1X,F9.3),1X,I6,1X,3(1X,F9.3))
    End Program g02mbfe