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)
      adj = 0.d0
      do irow = 1, nrow
         do ij = ifirst(irow)+1, ifirst(irow+1)
            adj(ivcol(ij)) = adj(ivcol(ij)) + zhz(ij)*xvec(irow) 
         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