Module g02jgfe_mod
! Mark 28.6 Release. NAG Copyright 2022.
! 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 ..
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