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