```!   G02JDF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

Module g02jdfe_mod

!     G02JDF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Accessibility Statements ..
Private
Public                           :: print_results
!     .. Parameters ..
Integer, Parameter, Public       :: nin = 5, nout = 6
Contains
Subroutine print_results(n,nff,nlsv,nrf,fixed,lfixed,nrndm,rndm,ldrndm,  &
nvpr,vpr,lvpr,gamma,effn,rnkx,ncov,lnlike,lb,id,ldid,b,se)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: lnlike
Integer, Intent (In)           :: effn, lb, ldid, ldrndm, lfixed,      &
lvpr, n, ncov, nff, nlsv, nrf,       &
nrndm, nvpr, rnkx
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: b(lb), gamma(nvpr+1), se(lb)
Integer, Intent (In)           :: fixed(lfixed), id(ldid,lb),          &
rndm(ldrndm,nrndm), vpr(lvpr)
!       .. Local Scalars ..
Integer                        :: aid, i, k, l, ns, nv, p, pb, tb,     &
tdid, vid
Character (120)                :: pfmt, tfmt
!       .. Executable Statements ..
!       Display the output
Write (nout,*) 'Number of observations (N)                    = ', n
Write (nout,*) 'Number of random factors (NRF)                = ', nrf
Write (nout,*) 'Number of fixed factors (NFF)                 = ', nff
Write (nout,*) 'Number of subject levels (NLSV)               = ',     &
nlsv
Write (nout,*) 'Rank of X (RNKX)                              = ',     &
rnkx
Write (nout,*) 'Effective N (EFFN)                            = ',     &
effn
Write (nout,*) 'Number of nonzero variance components (NCOV) = ', ncov

Write (nout,99990) 'Parameter Estimates'
tdid = nff + nrf*nlsv

If (nrf>0) Then
Write (nout,*)
Write (nout,99990) 'Random Effects'
End If
pb = -999
pfmt = ' '
Do k = 1, nrf*nlsv
tb = id(1,k)
If (tb/=-999) Then
vid = id(2,k)
nv = rndm(1,tb)
ns = rndm(3+nv,tb)
Write (tfmt,*)(id(3+l,k),l=1,ns)
If (pb/=tb .Or. tfmt/=pfmt) Then
If (k/=1) Then
Write (nout,*)
End If
Write (nout,99991) ' Subject: ', ('Variable ',rndm(3+nv+l,tb),   &
' (Level ',id(3+l,k),')',l=1,ns)
End If
If (vid==0) Then
!             Intercept
Write (nout,99994) b(k), se(k)
Else
!             VID'th variable specified in RNDM
aid = rndm(2+vid,tb)
If (id(3,k)==0) Then
Write (nout,99992) aid, b(k), se(k)
Else
Write (nout,99993) aid, id(3,k), b(k), se(k)
End If
End If
pfmt = tfmt
End If
pb = tb
End Do

If (nff>0) Then
Write (nout,*)
Write (nout,99990) 'Fixed Effects'
End If
Do k = nrf*nlsv + 1, tdid
If (vid/=-999) Then
vid = id(2,k)
If (vid==0) Then
!             Intercept
Write (nout,99997) b(k), se(k)
Else
!             VID'th variable specified in FIXED
aid = fixed(2+vid)
If (id(3,k)==0) Then
Write (nout,99995) aid, b(k), se(k)
Else
Write (nout,99996) aid, id(3,k), b(k), se(k)
End If
End If
End If
End Do

Write (nout,*)
Write (nout,*) 'Variance Components'
Write (nout,*) ' Estimate      Parameter        Subject'
Do k = 1, nvpr
p = 0
Do tb = 1, nrndm
nv = rndm(1,tb)
ns = rndm(3+nv,tb)
If (rndm(2,tb)==1) Then
p = p + 1
If (vpr(p)==k) Then
End If
End If
Do i = 1, nv
p = p + 1
If (vpr(p)==k) Then
Write (nout,99989,Advance='NO') rndm(2+i,tb),                  &
(rndm(3+nv+l,tb),l=1,ns)
End If
End Do
End Do
Write (nout,*)
End Do
Write (nout,*)
Write (nout,99998) 'SIGMA**2         = ', gamma(nvpr+1)
Write (nout,99998) '-2LOG LIKELIHOOD = ', lnlike

Return
99999   Format (1X,F10.5,5X)
99998   Format (1X,A,F15.5)
99997   Format (3X,'Intercept',20X,F10.4,1X,F10.4)
99996   Format (3X,'Variable ',I2,' (Level ',I2,')',7X,F10.4,1X,F10.4)
99995   Format (3X,'Variable ',I2,18X,F10.4,1X,F10.4)
99994   Format (5X,'Intercept',18X,F10.4,1X,F10.4)
99993   Format (5X,'Variable ',I2,' (Level ',I2,')',5X,F10.4,1X,F10.4)
99992   Format (5X,'Variable ',I2,16X,F10.4,1X,F10.4)
99991   Format (1X,A,4(A,I2,A,I2,A,1X))
99990   Format (1X,A)
99989   Format (1X,'Variable',1X,I2,5X,'Variables',1X,100(I2,1X))
99988   Format (1X,'Intercept',7X,'Variables',1X,100(I2,1X))
End Subroutine print_results
End Module g02jdfe_mod
Program g02jdfe

!     G02JDF Example Main Program

!     .. Use Statements ..
Use nag_library, Only: g02jcf, g02jdf, nag_wp
Use g02jdfe_mod, Only: nin, nout, print_results
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: lnlike
Integer                          :: effn, i, ifail, j, lb, ldcxx, ldcxz, &
ldczz, lddat, ldid, ldrndm, lfixed,  &
licomm, liopt, lrcomm, lropt, lvpr,  &
lwt, n, ncol, ncov, nff, nl, nlsv,   &
nrf, nrndm, nv, nvpr, nzz, rnkx
Character (1)                    :: weight
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: b(:), cxx(:,:), cxz(:,:), czz(:,:),  &
dat(:,:), gamma(:), rcomm(:),        &
ropt(:), se(:), wt(:), y(:)
Integer, Allocatable             :: fixed(:), icomm(:), id(:,:),         &
iopt(:), levels(:), rndm(:,:),       &
vpr(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: max
!     .. Executable Statements ..
Write (nout,*) 'G02JDF Example Program Results'
Write (nout,*)

!     Skip the heading in data file

!     Read in the problem size
Read (nin,*) weight, n, ncol, nrndm, nvpr

!     Set LFIXED and LDRNDM to maximum value they could
!     be for this dataset
lfixed = ncol + 2
ldrndm = 3 + 2*ncol

If (weight=='W' .Or. weight=='w') Then
lwt = n
Else
lwt = 0
End If
lddat = n
Allocate (dat(lddat,ncol),levels(ncol),y(n),wt(lwt),fixed(lfixed),       &
rndm(ldrndm,nrndm))

!     Read in the number of levels associated with each of the
!     independent variables

!     Read in the fixed part of the model

!     Number of variables
nv = fixed(1)

!     Intercept

!     Variable IDs
If (nv>0) Then
End If

!     Read in the random part of the model
lvpr = 0
Do j = 1, nrndm

!       Number of variables and intercept
nv = rndm(1,j)

!       Variable IDs
If (nv>0) Then
End If

!       Number of subject variables
nl = rndm(nv+3,j)

!       Subject variable IDs
If (nl>0) Then
End If
lvpr = lvpr + rndm(2,j) + nv
End Do

!     Read in the dependent and independent data
If (lwt>0) Then
Else
End If

licomm = 2
lrcomm = 0
Allocate (icomm(licomm),rcomm(lrcomm))

!     Call the initialization routine once to get LRCOMM and LICOMM
ifail = 0
Call g02jcf(weight,n,ncol,dat,lddat,levels,y,wt,fixed,lfixed,nrndm,rndm, &
ldrndm,nff,nlsv,nrf,rcomm,lrcomm,icomm,licomm,ifail)

!     Reallocate ICOMM and RCOMM
licomm = icomm(1)
lrcomm = icomm(2)
Deallocate (icomm,rcomm)
Allocate (icomm(licomm),rcomm(lrcomm))

!     Pre-process the data
ifail = 0
Call g02jcf(weight,n,ncol,dat,lddat,levels,y,wt,fixed,lfixed,nrndm,rndm, &
ldrndm,nff,nlsv,nrf,rcomm,lrcomm,icomm,licomm,ifail)

!     Use the default options
liopt = 0
lropt = 0

!     Calculate LDID
ldid = 0
Do i = 1, nrndm
nv = rndm(1,i)
ldid = max(rndm(3+nv,i),ldid)
End Do
ldid = ldid + 3

lb = nff + nrf*nlsv
nzz = nrf*nlsv
ldczz = nzz
ldcxx = nff
ldcxz = nff
Allocate (vpr(lvpr),gamma(nvpr+1),id(ldid,lb),b(lb),se(lb),              &
czz(ldczz,nzz),cxx(ldcxx,nff),cxz(ldcxz,nzz),iopt(liopt),ropt(lropt))

!     Read in VPR

!     Read in GAMMA

!     Perform the analysis
ifail = -1
Call g02jdf(lvpr,vpr,nvpr,gamma,effn,rnkx,ncov,lnlike,lb,id,ldid,b,se,   &
czz,ldczz,cxx,ldcxx,cxz,ldcxz,rcomm,icomm,iopt,liopt,ropt,lropt,ifail)
If (ifail/=0 .And. ifail<100) Then
Go To 100
End If

!     Display results
Call print_results(n,nff,nlsv,nrf,fixed,lfixed,nrndm,rndm,ldrndm,nvpr,   &
vpr,lvpr,gamma,effn,rnkx,ncov,lnlike,lb,id,ldid,b,se)

100   Continue

End Program g02jdfe
```