```    Program g02mbfe

!     G02MBF Example Program Text

!     Mark 26 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 in the problem size

!     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
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))

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
```