!   G02JDF Example Program Text
!   Mark 26.1 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
          Write (nout,99999,Advance='NO') gamma(k)
          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
                Write (nout,99988,Advance='NO')(rndm(3+nv+l,tb),l=1,ns)
              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 g02jdfe_mod, Only: nin, nout, print_results
      Use nag_library, Only: g02jcf, g02jdf, nag_wp
!     .. 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 (nin,*)

!     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 (nin,*) levels(1:ncol)

!     Read in the fixed part of the model
      Read (nin,*)

!     Number of variables
      Read (nin,*) fixed(1)
      nv = fixed(1)

!     Intercept
      Read (nin,*) fixed(2)

!     Variable IDs
      If (nv>0) Then
        Read (nin,*) fixed(3:(nv+2))
      End If

!     Read in the random part of the model
      lvpr = 0
      Do j = 1, nrndm
!       Skip header
        Read (nin,*)

!       Number of variables and intercept
        Read (nin,*) rndm(1,j)
        Read (nin,*) rndm(2,j)
        nv = rndm(1,j)

!       Variable IDs
        If (nv>0) Then
          Read (nin,*)(rndm(i,j),i=3,nv+2)
        End If

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

!       Subject variable IDs
        If (nl>0) Then
          Read (nin,*)(rndm(i,j),i=nv+4,nv+nl+3)
        End If
        lvpr = lvpr + rndm(2,j) + nv
      End Do

!     Read in the dependent and independent data
      If (lwt>0) Then
        Read (nin,*)(y(i),dat(i,1:ncol),wt(i),i=1,n)
      Else
        Read (nin,*)(y(i),dat(i,1:ncol),i=1,n)
      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 (nin,*) vpr(1:lvpr)

!     Read in GAMMA
      Read (nin,*) gamma(1:nvpr)

!     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