WOMBAT – A program for Mixed Model Analyses by Restricted Maximum Likelihood

harville.f90
``` module inv_harville                                       !© karin meyer
integer, save                            :: nrczhz, maxrec
integer, dimension(:),allocatable, save  :: ifirst, ivcol, inext
real(8), dimension(:),allocatable, save  :: dia, zhz
integer, dimension(:), allocatable, save :: ivec
real(8), dimension(:), allocatable, save :: xvec

contains

! ==========================================================
subroutine harville( nrow, ngibbs, nburn, iseed1, iseed2 )
! ==========================================================

! Purpose : Obtain inverse of large, sparse matrix (symmetric & positive
!           definite) using Gibbs sampling

!           Assumes non-zero elements of the matrix are stored
!           in linked list format, which has been sorted
!               dia: diagonal elements
!               zhz: off-triagonal elements (lower triangle only)
!               ifirst: position of first element in a row
!               inext: position of next off-diagonal element
!               ivcol: column no. for off-diagonal elements

!          ** Version to give diagonal elements of inverse only **

! Reference: Harville, D.A. (1999) "Use of the Gibbs sampler to
!            invert large, possibly sparse, positive definite matrices"
!            Linear Algebra and its Applications 289 : 203-224

! Arguments:
!            nrow - no. of rows in matrix
!            ngibbs - no. of Gibbs samples to be drawn
!            nburn -  no. of samples to be discarded in "burn-in" phase
!            is1, is2 - seed values to  intialize random no. generator
!------------------------------------------------------------- km 11/2002 --------

implicit none
integer, intent(in)       :: nrow, ngibbs, nburn, iseed1, iseed2
real(8), dimension(nrow)  :: xvec, dd, dinv, adj, sdev
real(8)                   :: aa, xbar, xx
integer                   :: irow, igibbs, ij, is1, is2

! function to generate random normal deviate from the ranlibf library
!                           (replace with similiar function as desired)
real, external            :: snorm

! initialize random number generator from ranlibf library (replace if desired)
call setall( iseed1, iseed2 )

! initialize
dinv = 0.d0; xvec = 0.d0
do irow = 1, nrow
if( dia(irow) <= 0 ) then
write(*,*)'Error in subroutine HARVILLE: zero/negative diagonal element'
stop 'Programmed stop'
end if
dd(irow) = 1.d0/dia(irow); sdev(irow) = sqrt( dd(irow) )
end do

! start of gibbs sampler
do igibbs = 1, ngibbs
!     ... adjustments for previous samples (upper triangle of matrix)
do irow = 1, nrow
do ij = ifirst(irow)+1, ifirst(irow+1)
end do
end do
do irow = 1, nrow
!        ... adjustments for samples from this iterate (lower triangle)
aa = 0.d0
do ij = ifirst(irow)+1, ifirst(irow+1)
aa = aa + zhz(ij)*xvec(ivcol(ij))
end do
!        ... conditional mean
xbar = -( aa + adj(irow) ) * dd(irow)
!        ... new sample
xvec(irow) = xbar + sdev(irow) * snorm()
!        ... accumulate
if( igibbs > nburn) dinv(irow) = dinv(irow) + xbar * xbar
end do ! irow
end do ! igibbs

! inverse
xx = 1.d0 / dble( ngibbs - nburn )
do irow = 1, nrow
dia(irow) = dd(irow) + dinv(irow) * xx
end do

return
end subroutine harville

! ==============================
subroutine lnkrow ( irow, nn )
! ==============================

! routine to add elements to a row in a sparse matrix held in linked list
!     arguments:
!       irow - row no.
!       nn   - no. of elements to be added
!     information passed in module harville:
!       ivec(nn) - vector of column no.s (must be in ascending order)
!       xvec(nn) - vector of elements to be added
!------------------------------------------------------ km  ------------

implicit none
integer, intent(in) :: irow, nn
logical, intent(in) :: opt
integer             :: kk, ipre, ipcol, iplace, icol, i

if( nn < 1 ) return

! check that column no.s are in ascending sequence
ipcol = 0
do kk = 1, nn
icol = ivec(kk)
if( icol <= ipcol ) then
write(*,*)'subroutine "lnkrow" : col. no.s out of order !'
write(*,*)' row no. =',irow
do  i = 1, nn
write(*,*)' column   ',i,ivec(i),xvec(i)
end do
stop 'Programmed stop'
end if
ipcol = icol
end do

ipre = 0; iplace = ifirst(irow)
do kk = 1, nn
icol = ivec(kk)
!    go to next element in this row
10   if( iplace > 0 .and. ivcol(iplace) < icol )then
ipre = iplace; iplace = inext(iplace); go to 10
!    element found
else if( iplace > 0 .and. icol == ivcol(iplace) )then
zhz(iplace) = zhz(iplace) + xvec(kk)
ipre = iplace; iplace = inext(iplace)
!    element (irow,icol) has not been accessed before
else
nrczhz = nrczhz + 1
!        check that array size is sufficient
if( nrczhz > maxrec ) then
write(*,*)'Subroutine "LNKROW" : DIMENSIONS EXCEEDED !!'
write(*,*)'max. no. of elements =', maxrec
write(*,*)'row no.              =', irow
write(*,*)'no. elements in row  =', nn
write(*,*)'column no.           =', icol
write(*,*)'element in row no.   =', kk
stop 'reset parameter "maxrec" '
end if

zhz(nrczhz) = xvec(kk); ivcol(nrczhz) = icol
if( iplace == 0 .and. ifirst(irow) == 0 )then
!            first element in this row
ifirst(irow) = nrczhz; inext(nrczhz) = 0
else if( iplace == 0 .and. ifirst(irow) /= 0 ) then
!            all previous elements in the row have lower column no.
inext(ipre) = nrczhz; inext(nrczhz)=0
else if(iplace /= 0 .and. ipre == 0 ) then
!            fit element in "right" place, swapping indicators
inext(nrczhz) = ifirst(irow); ifirst(irow) = nrczhz
else
inext(nrczhz) = inext(ipre); inext(ipre) = nrczhz
end if
ipre = nrczhz; iplace = inext(nrczhz)
end if
end do

return
end subroutine lnkrow

! ===========================
subroutine lnksrt( neqns )
! ===========================

! routine to sort a sparse matrix held in linked list
!     arguments:
!       newns - no. of rows
!------------------------------------------------------ km -----------

implicit none
integer, intent(in) :: neqns
integer             :: newplc, irow, iplace, jplace, nrow, icol
real(8)             :: zave

newplc = 0; nrow = 0

do irow = 1, neqns
iplace = ifirst(irow); ifirst(irow) = nrow
do while( iplace > 0 )
newplc = newplc + 1
do while( iplace < newplc .and. iplace > 0 )
iplace = inext(iplace)
end do
jplace = inext(iplace)
if( iplace /= newplc ) then
zave = zhz(newplc)
zhz(newplc) = zhz(iplace)
zhz(iplace) = zave
inext(iplace) = inext(newplc)
inext(newplc) = iplace
icol = ivcol(iplace)
ivcol(iplace) = ivcol(newplc)
ivcol(newplc) = icol
end if
nrow = nrow+1; iplace = jplace
end do
end do ! irow
ifirst(neqns+1) = nrczhz

if( nrow /= nrczhz ) then
write(*,*)'Subroutine LNKSRT: Error in sorting',nrow,nrczhz,neqns
stop 'Programmed stop'
end if

return
end subroutine lnksrt

end module inv_harville```