Program g12zafe

!     G12ZAF Example Program Text.

!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: g11caf, g12zaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: dev, tol
      Integer                          :: cm, i, ifail, ip, iprint, ldz, lisi, &
                                          lwk, m, maxit, mxn, n, nd, ns, num,  &
                                          nxs
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:), cov(:), sc(:), se(:), t(:),    &
                                          tp(:), wk(:), x(:,:), z(:,:)
      Integer, Allocatable             :: cnt(:), ic(:), id(:), irs(:),        &
                                          isi(:), isz(:), ixs(:), nca(:),      &
                                          nct(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: count, maxval
!     .. Executable Statements ..
      Write (nout,*) 'G12ZAF Example Program Results'
      Write (nout,*)

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

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

      If (ns>0) Then
        lisi = n
      Else
        lisi = 0
      End If

      ldz = n
      Allocate (z(ldz,m),isz(m),t(n),ic(n),isi(lisi),tp(n),irs(n))

!     Read in the data
      If (ns>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),i=1,n)
      End If

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

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

!     Call the routine once to calculate size of MXN ...

!     Dummy allocation
      mxn = 0
      Allocate (x(mxn,ip),id(mxn),ixs(mxn))

!     Call G12ZAF to calculate MXN
      ifail = 1
      Call g12zaf(n,m,ns,z,ldz,isz,ip,t,ic,isi,num,ixs,nxs,x,mxn,id,nd,tp,irs, &
        ifail)
      If (ifail/=0 .And. ifail/=3) Then
        Go To 100
      End If

!     Required size for MXN is returned in NUM, so reallocate memory
      mxn = num
      Deallocate (x,id,ixs)
      Allocate (x(mxn,ip),id(mxn),ixs(mxn))

!     Create risk set
      ifail = 0
      Call g12zaf(n,m,ns,z,ldz,isz,ip,t,ic,isi,num,ixs,nxs,x,mxn,id,nd,tp,irs, &
        ifail)

      Allocate (cnt(nxs),b(ip),se(ip),sc(ip),nca(nxs),nct(nxs),cov(ip*(ip+     &
        1)/2))

!     Set tolerance
      tol = 1.0E-5_nag_wp

!     Read in initial parameter estimates
      Read (nin,*) b(1:ip)

!     Count the number of observations in each stratum
      cnt(1:nxs) = 0
      Do i = 1, num
        cnt(ixs(i)) = cnt(ixs(i)) + 1
      End Do
      cm = maxval(cnt(1:nxs))

      lwk = ip*num + (cm+1)*(ip+1)*(ip+2)/2 + cm
      Allocate (wk(lwk))

!     Get parameter estimates from conditional logistic analysis
      ifail = 0
      Call g11caf(num,ip,nxs,x,mxn,isz,ip,id,ixs,dev,b,se,sc,cov,nca,nct,tol,  &
        maxit,iprint,wk,lwk,ifail)

!     Display results
      Write (nout,*) ' Parameter      Estimate', '       Standard Error'
      Write (nout,*)
      Write (nout,99999)(i,b(i),se(i),i=1,ip)

100   Continue

99999 Format (I6,10X,F8.4,10X,F8.4)
    End Program g12zafe