!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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