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

harville_test.f90
```!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! example matrix
!    2.000   1.000   1.000   0.000   0.000
!    1.000   3.000   0.000   1.000   1.000
!    1.000   0.000   3.000   0.000   1.000
!    0.000   1.000   0.000   1.000   1.000
!    0.000   1.000   1.000   1.000   2.000
! exact inverse
!    1.000  -0.500  -0.500   0.000   0.500
!   -0.500   0.750   0.250  -0.500  -0.250
!   -0.500   0.250   0.750   0.500  -0.750
!    0.000  -0.500   0.500   3.000  -1.500
!    0.500  -0.250  -0.750  -1.500   1.750
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! ======================
program test_harville
! ======================

use inv_harville

implicit none
integer, parameter                      :: maxneq = 10, maxzhz = 40
integer                                 :: ii, i, neqns, is1, is2

! set 'seeds' for random no. generator
is1 = 123456; is2 = 987650

maxrec = maxzhz; neqns = 5

! allocate vectors required
allocate(ifirst(neqns+1), inext(maxzhz), ivcol(maxzhz), dia(neqns), &
&           zhz(maxzhz), xvec(neqns), ivec(neqns), stat=ii)
if( ii /= 0) stop 'alloc lnklst'

! initialize
nrczhz = 0; ifirst = 0

! store example matrix (lower triangle)  in linked list format
dia(:neqns) = (/ 2.d0, 3.d0, 3.d0, 1.d0, 2.d0 /)
xvec = 1.d0; ivec(1) = 1
call lnkrow( 2, 1 ) ! row 2
call lnkrow( 3, 1 ) ! row 3
ivec(1) = 2; call lnkrow( 4, 1 ) ! row 4
ivec(:3)= (/ 2, 3, 4 /);  call lnkrow( 5, 3 ) ! row 5
write(*,*)'no. of non-zero elements: nrczhz=', nrczhz

! sort matrix
call lnksrt( neqns )

! don't need 'next' element any longer after sorting
deallocate(inext, stat=ii ); if( ii /= 0 ) stop 'deall inext'

! obtain approximate inverse
call harville( neqns, 55000, 5000, is1, is2 )

! write out diagonal elements of inverse
det = 0.d0
do i = 1, neqns
write(*,*) i, dia(i)
end do

end program test_harville```