!   G02JDF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    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
!      .. Parameters ..
       INTEGER, PARAMETER                  :: nin = 5, nout = 6
    CONTAINS
!dec$ attributes offload : mic :: print_results
       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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 non-zero 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))
!dec$ attributes offload : mic :: print_results
       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 Functions ..
       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 + 1
       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 initialisation 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 20
       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)

20     CONTINUE

    END PROGRAM g02jdfe
