```    Program g12bafe

!     G12BAF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. Use Statements ..
Use nag_library, Only: g02gcf, g12baf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: dev, tol
Integer                          :: i, idf, ifail, ip, ip1, iprint,      &
irank, ldv, ldz, lisi, lomega, m,    &
maxit, n, nd, ndmax, ns
Character (1)                    :: offset
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: b(:), cov(:), omega(:), res(:),      &
sc(:), se(:), sur(:,:), t(:), tp(:), &
v(:,:), wk(:), y(:), z(:,:)
Integer, Allocatable             :: ic(:), isi(:), isz(:), iwk(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: count, eoshift, log, max, real
!     .. Executable Statements ..
Write (nout,*) 'G12BAF Example Program Results'
Write (nout,*)

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

!     Read in problem size
Read (nin,*) n, m, ns, maxit, iprint, offset

If (offset=='Y' .Or. offset=='y') Then
lomega = n
Else
lomega = 0
End If
If (ns>0) Then
lisi = n
Else
lisi = 0
End If
ldz = n
ndmax = n
ldv = n
Allocate (z(ldz,m),isz(m),t(n),ic(n),omega(lomega),isi(lisi),res(n),     &
tp(ndmax),sur(ndmax,max(ns,1)),iwk(2*n),y(n))

!     Read in the data
If (ns>0) Then
If (lomega==0) Then
Read (nin,*)(t(i),z(i,1:m),ic(i),isi(i),i=1,n)
Else
Read (nin,*)(t(i),z(i,1:m),ic(i),isi(i),omega(i),i=1,n)
End If
Else
If (lomega==0) Then
Read (nin,*)(t(i),z(i,1:m),ic(i),i=1,n)
Else
Read (nin,*)(t(i),z(i,1:m),ic(i),omega(i),i=1,n)
End If
End If

!     Read in the variable indication
Read (nin,*) isz(1:m)

!     Calculate number of parameters in the model
ip = count(isz(1:m)>0)

!     We are fitting a mean in the exponential model, so arrays used by G02GCF
!     need to be based on IP + 1
ip1 = ip + 1
Allocate (b(ip1),se(ip1),sc(ip),cov(ip1*(ip1+1)/2),wk(ip1*(ip1+          &
9)/2+n),v(ldv,ip1+7))

!     Specify tolerance to use
tol = 5.0E-5_nag_wp

!     Get initial estimates by fitting an exponential model
Do i = 1, n
y(i) = 1.0E0_nag_wp - real(ic(i),kind=nag_wp)
v(i,7) = log(t(i))
End Do

!     Fit exponential model
ifail = -1
Call g02gcf('L','M','Y','U',n,z,ldz,m,isz,ip1,y,res,0.0E0_nag_wp,dev,    &
idf,b,irank,se,cov,v,ldv,tol,maxit,0,0.0E0_nag_wp,wk,ifail)
If (ifail/=0) Then
If (ifail<5) Then
Go To 100
End If
End If

!     Check exponential model was of full rank
If (irank/=ip1) Then
Write (nout,*) ' WARNING: covariates not of full rank'
End If

!     Move all parameter estimates down one so as to drop the parameter
!     estimate for the mean.
b(1:ip1) = eoshift(b(1:ip1),1)

!     Fit Cox proportional hazards model
ifail = -1
Call g12baf('No-offset',n,m,ns,z,ldz,isz,ip,t,ic,omega,isi,dev,b,se,sc,  &
cov,res,nd,tp,sur,ndmax,tol,maxit,iprint,wk,iwk,ifail)
If (ifail/=0) Then
If (ifail<5) Then
Go To 100
End If
End If

!     Display results
Write (nout,*) ' Parameter      Estimate', '       Standard Error'
Write (nout,*)
Write (nout,99999)(i,b(i),se(i),i=1,ip)
Write (nout,*)
Write (nout,99998) ' Deviance = ', dev
Write (nout,*)
Write (nout,*) '    Time     Survivor Function'
Write (nout,*)
ns = max(ns,1)
Write (nout,99997)(tp(i),sur(i,1:ns),i=1,nd)

100   Continue

99999 Format (I6,10X,F8.4,10X,F8.4)
99998 Format (A,E13.4)
99997 Format (F10.0,5X,F8.4)
End Program g12bafe
```