NAG Library Manual, Mark 27.3
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Module g02jgfe_mod
!     Mark 27.3 Release. NAG Copyright 2021.
!     G02JGF Example Program Module:
!            Parameters and Utility Routines

!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: read_line
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
    Contains
      Subroutine read_line(v1)
!       Read in a line from NIN and remove any comments

!       .. Scalar Arguments ..
        Character (*), Intent (Out)    :: v1
!       .. Local Scalars ..
        Integer                        :: pend
!       .. Intrinsic Procedures ..
        Intrinsic                      :: adjustl, index
!       .. Executable Statements ..
        Continue

        Read (nin,'(A200)') v1
        pend = index(v1,'::')
        If (pend/=0) Then
          v1 = v1(1:pend-1)
        End If
        v1 = adjustl(v1)

        Return
      End Subroutine read_line
    End Module g02jgfe_mod
    Program g02jgfe

!     G02JGF Example Program Text

!     .. Use Statements ..
      Use g02jgfe_mod, Only: nin, nout, read_line
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: g02jff, g02jgf, g02jhf, g02zlf, g22yaf, g22ybf,   &
                             g22ydf, g22zaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (c_ptr)                     :: hddesc, hfixed, hlmm, hsform
      Real (Kind=nag_wp)               :: lnlike, rvalue
      Integer                          :: blicomm, blk, blrcomm, effn, fnlsv,  &
                                          i, ifail, ip, lb, ldcxx, ldcxz,      &
                                          ldczz, lddat, licomm, lisx, lrcomm,  &
                                          lvinfo, lvnames, lwt, n, nblk, ncov, &
                                          nff, nrf, nrndm, nvar, nvpr, nxx,    &
                                          nzz, optype, rnkx, rnlsv, sddat,     &
                                          tlicomm, tlrcomm, totaln
      Character (20)                   :: cvalue
      Character (200)                  :: formula, line
      Character (1)                    :: intcpt, weight
!     .. Local Arrays ..
      Type (c_ptr), Allocatable        :: hrndm(:)
      Real (Kind=nag_wp), Allocatable  :: b(:), brcomm(:), cxx(:,:), cxz(:,:), &
                                          czz(:,:), dat(:,:), gamma(:),        &
                                          rcomm(:), se(:), trcomm(:), wt(:),   &
                                          y(:)
      Integer, Allocatable             :: bicomm(:), icomm(:), isx(:),         &
                                          levels(:), ticomm(:), vinfo(:)
      Character (50), Allocatable      :: vnames(:), vplab(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: adjustl, allocated, move_alloc,      &
                                          size, trim
!     .. Executable Statements ..
      Write (nout,*) 'G02JGF Example Program Results'
      Write (nout,*)
      Flush (nout)

      hfixed = c_null_ptr
      hddesc = c_null_ptr
      hlmm = c_null_ptr
      hsform = c_null_ptr

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

!     Read in size of the data matrix and number of variable labels supplied
      Read (nin,*) weight, n, nvar, lvnames

!     Read in number of levels and names for the variables
      Allocate (levels(nvar),vnames(lvnames))
      Read (nin,*) levels(1:nvar)
      If (lvnames>0) Then
        Read (nin,*) vnames(1:lvnames)
      End If

!     Create a description of the data matrix
      ifail = 0
      Call g22ybf(hddesc,n,nvar,levels,lvnames,vnames,ifail)

!     Read in the formula for the fixed part of the model,
!     remove comments and parse it
      Call read_line(formula)
      ifail = 0
      Call g22yaf(hfixed,formula,ifail)

!     Read in the number of random statements
      Read (nin,*) nrndm

      Allocate (hrndm(nrndm))
      hrndm(:) = c_null_ptr

!     Read in the formula for the random parts of the model,
!     remove comments and parse them
      Do i = 1, nrndm
        Call read_line(formula)
        ifail = 0
        Call g22yaf(hrndm(i),formula,ifail)
      End Do

!     Skip heading
      Read (nin,*)

!     Read in the number of blocks of data
      Read (nin,*) nblk

!     Loop over each block of data
      totaln = 0
      Do blk = 1, nblk
!       Skip header
        Read (nin,*)

!       Read in the number of observations in this block
        Read (nin,*) n

!       Keep a running total of the number of observations processed
        totaln = totaln + n

        If (weight=='W' .Or. weight=='w') Then
          lwt = n
        Else
          lwt = 0
        End If
        lddat = n
        sddat = nvar

        Allocate (dat(lddat,sddat),y(n),wt(lwt))

!       Read in the dependent and independent data for this block
        If (lwt>0) Then
          Read (nin,*)(dat(i,1:nvar),wt(i),y(i),i=1,n)
        Else
          Read (nin,*)(dat(i,1:nvar),y(i),i=1,n)
        End If

        blicomm = 2
        blrcomm = 0
        If (allocated(bicomm)) Then
          Deallocate (bicomm,brcomm)
        End If
        Allocate (bicomm(blicomm),brcomm(blrcomm))

!       Prepare the data
        ifail = 1
        Call g02jff(hlmm,hddesc,hfixed,nrndm,hrndm,weight,n,y,wt,dat,lddat,    &
          sddat,fnlsv,nff,rnlsv,nrf,nvpr,brcomm,blrcomm,bicomm,blicomm,ifail)
        If (ifail/=0) Then
          If (ifail==191) Then
!           COMM arrays were not large enough, call the routine a second
!           time with the correct array sizes
            blicomm = bicomm(1)
            blrcomm = bicomm(2)
            Deallocate (bicomm,brcomm)
            Allocate (bicomm(blicomm),brcomm(blrcomm))
          End If

!         Re-run the routine, either to make use of the larger communication
!         arrays or print out information relevant to the error message
!         and terminate
          ifail = -1
          Call g02jff(hlmm,hddesc,hfixed,nrndm,hrndm,weight,n,y,wt,dat,lddat,  &
            sddat,fnlsv,nff,rnlsv,nrf,nvpr,brcomm,blrcomm,bicomm,blicomm,      &
            ifail)
          If (ifail/=0 .And. ifail/=34 .And. ifail/=102) Then
!           terminate on an error
            Stop
          End If
        End If

!       Combine the blocks of data
        If (blk==1) Then
!         This is the first block of data, so move the communication arrays
!         into ICOMM and RCOMM
          Call move_alloc(bicomm,icomm)
          Call move_alloc(brcomm,rcomm)
          licomm = blicomm
          lrcomm = blrcomm

        Else
!         Combine the current block with the previous blocks
          ifail = 1
          Call g02jgf(hlmm,rcomm,lrcomm,icomm,licomm,brcomm,bicomm,ifail)
          If (ifail/=0) Then
            If (ifail==52) Then
!             ICOMM and / or RCOMM are not sufficiently large to hold the
!             combined communication array
              If (licomm<icomm(1)) Then
!               ICOMM is not large enough, reallocate, preserving the contents
                tlicomm = licomm
                licomm = icomm(1)
                Call move_alloc(icomm,ticomm)
                Allocate (icomm(licomm))
                icomm(1:tlicomm) = ticomm(1:tlicomm)
                Deallocate (ticomm)
              End If
              If (lrcomm<icomm(2)) Then
!               RCOMM is not large enough, reallocate, preserving the contents
                tlrcomm = lrcomm
                lrcomm = icomm(2)
                Call move_alloc(rcomm,trcomm)
                Allocate (rcomm(lrcomm))
                rcomm(1:tlrcomm) = trcomm(1:tlrcomm)
                Deallocate (trcomm)
              End If
            End If
            ifail = 0

!           Re-run the routine. If IFAIL=32 this will re-run with the larger
!           communication sizes, otherwise it will print out information
!           relevant to the error message and terminate
            Call g02jgf(hlmm,rcomm,lrcomm,icomm,licomm,brcomm,bicomm,ifail)
          End If
        End If

!       Free up the data arrays
        Deallocate (dat,y,wt)

!       Free up the communication arrays for the current block
        If (allocated(bicomm)) Then
          Deallocate (bicomm,brcomm)
        End If
      End Do

!     Clear up the handles that are no longer required
      ifail = 0
      Call g22zaf(hddesc,ifail)
      Call g22zaf(hfixed,ifail)
      Do i = 1, nrndm
        Call g22zaf(hrndm(i),ifail)
      End Do

!     Get the number of fixed and random factors and the number of
!     variance components for the combined dataset
      ifail = 0
      Call g02zlf('nff',nff,rvalue,cvalue,optype,icomm,rcomm,ifail)
      Call g02zlf('nrf',nrf,rvalue,cvalue,optype,icomm,rcomm,ifail)
      Call g02zlf('rnlsv',rnlsv,rvalue,cvalue,optype,icomm,rcomm,ifail)
      Call g02zlf('fnlsv',fnlsv,rvalue,cvalue,optype,icomm,rcomm,ifail)
      Call g02zlf('nvpr',nvpr,rvalue,cvalue,optype,icomm,rcomm,ifail)

      nzz = nrf*rnlsv
      nxx = nff*fnlsv
      lb = nxx + nzz

!     We don't want the output from CXX, CZZ and CXZ
      Allocate (czz(0,0),cxx(0,0),cxz(0,0))
      ldczz = size(czz,1)
      ldcxx = size(cxx,1)
      ldcxz = size(cxz,1)

!     Allocate the rest of the output arrays
      Allocate (gamma(nvpr+1),b(lb),se(lb))

!     Use MIVQUE for the initial values for gamma
      gamma(:) = -1.0_nag_wp

!     Fit the model
      ifail = -1
      Call g02jhf(hlmm,nvpr,gamma,effn,rnkx,ncov,lnlike,lb,b,se,czz,ldczz,cxx, &
        ldcxx,cxz,ldcxz,rcomm,icomm,ifail)
      If (ifail/=0 .And. ifail/=1001 .And. ifail/=1002 .And. ifail/=1003 .And. &
        ifail/=1004) Then
!       terminate if there is an error
        Stop
      End If

!     Print some results
      Write (nout,99998) 'Number of observations           = ', totaln
      Write (nout,99998) 'Total number of random factors   = ', nzz
      Write (nout,99998) 'Total number of fixed factors    = ', nxx
      Write (nout,99998) 'Effective N                      = ', effn
      Write (line,99996) lnlike
      line = adjustl(line)
      Write (nout,99999) '-2Log Likelihood                 = ', trim(line)
      Write (line,99996) gamma(nvpr+1)
      line = adjustl(line)
      Write (nout,99999) 'Sigma**2                         = ', trim(line)

      If (nvpr>0) Then
!       Get labels for the variance component estimates
        lvinfo = 0
        lisx = 0
        Allocate (vinfo(lvinfo),vplab(nvpr),isx(lisx))
        Call g22ydf('V',hsform,hlmm,intcpt,ip,lisx,isx,nvpr,vplab,lvinfo,      &
          vinfo,ifail)

        Write (nout,*)
        Write (nout,99999) 'Variance Components'
        Write (nout,99995) 'Component', 'Estimate'
        Do i = 1, nvpr
          Write (*,99997) adjustl(vplab(i)), gamma(i)
        End Do
      End If

!     Clear up the remaining handles
      Call g22zaf(hlmm,ifail)

99999 Format (1X,A,A)
99998 Format (1X,A,I0)
99997 Format (1X,A25,2(5X,F10.4))
99996 Format (F10.4)
99995 Format (A10,22X,A10)
    End Program g02jgfe