! ================================================
  subroutine dkmxhf( a, det, zero, nrank, n, iopt)
! ================================================
 
  implicit none
  integer, intent(in)                          :: n     ! no. of rows/columns
  logical, intent(in)                          :: iopt  ! .true./.false. : print warnings
  real(8), dimension(n*(n+1)/2), intent(inout) :: a     ! matrix; upper triangle, stored row-wise
  real(8), intent(in)                          :: zero  ! operational zero
  real(8), intent(out)                         :: det   ! log determinant
  integer, intent(out)                         :: nrank ! rank of matrix
  real(8), dimension(n)                        :: v, w
  integer, dimension(n)                        :: iflag
  real(8)                                      :: ww, xx, dmax, amax, bmax, dimax
  integer                                      :: ii, i1, n1, il, imax, &
&                                                 imaxm1,imaxp1,ij,j,i,neg
 
! matrix is a scalar
  if( n == 1 )then
      if(a(1) > zero)then
         det=dlog( a(1) ); a(1)=1.d0/a(1); nrank=1
      else
         a(1)=0.d0; nrank=0; det=0.d0
      end if
      return
   end if
 
!  initialize
   neg=0; det=0.d0; iflag=0;  n1=n+1
 
!  pick out diagonal elements
   ii=-n
   do i=1,n
      ii=ii+n1; w(i)=a(ii); ii=ii-i
   end do
 
!  gaussian elimination steps
   do ii=1,n
 
!     ... find diag. element with largest absolute value (pivot)
      dmax=0.d0;  amax=0.d0
      do i=1,n
         if(iflag(i) /= 0 ) cycle
         bmax=dabs(w(i))
         if(bmax.gt.amax)then
            dmax=w(i); amax=bmax; imax=i
         end if
      end do
 
!     ... check for singularity
      if(amax.le.zero)go to 11
      if(iopt .and. amax.lt.1.d-5) print *,'small pivot ',ii,dmax
!     ... set flag
      iflag(imax)=ii
!     accumulate log determinant
      det=det+dlog(amax)
      if(dmax.lt.0.d0)then
         neg=neg+1
         if(iopt)print *,'negative pivot, ignore sign for log det',ii,imax,dmax
      end if
      dimax=1.d0/dmax; imaxm1=imax-1; imaxp1=imax+1
 
!     ... pick out elements for row imax
      il=imax-n
      do i=1,imaxm1
         il=il+n1-i; v(i)=a(il)
      end do
      il=il+n1-imax
      do i=imaxp1,n
         il=il+1; v(i)=a(il)
      end do
 
!     ... transform matrix
      ij=0
      do i=1,imaxm1
         ww = v(i)
         if(ww /= 0.d0)then
            xx = ww*dimax
            ij = ij+1; w(i) = w(i) - ww*xx
            do j=i+1,imaxm1
               ij=ij+1; if( v(j) /= 0) a(ij) = a(ij)-xx*v(j)
            end do
!           element a(i,imax)
            ij=ij+1; a(ij) = xx
            do j = imaxp1,n
               ij=ij+1; if( v(j) /= 0) a(ij) = a(ij)-xx*v(j)
            end do
         else
            ij=ij+n1-i ! update address counter
         end if
      end do
 
!     row imax
      ij=ij+1; w(imax)=-dimax
      do j =imaxp1,n
         ij=ij+1; a(ij)=v(j)*dimax
      end do
      do i=imaxp1,n
         ww=v(i)
         if(ww /= 0.d0)then
            xx=ww*dimax; ij=ij+1; w(i)=w(i)-ww*xx
            do j=i+1,n
               ij=ij+1; a(ij)=a(ij)-xx*v(j)
            end do
         else
            ij=ij+n1-i
         end if
      end do
 
  end do ! loop over rows/columns
 
! store diagonals back & reverse signs
      ij=0
      do i=1,n
         ij=ij+1; ww=-w(i); w(i)=ww; a(ij)=ww
         do j=i+1,n
            ij=ij+1; a(ij)=-a(ij)
         end do
      end do
      nrank=n
 
300 if(iopt) then
       if(neg.gt.0) write(*,*)'no. of negative pivots =', neg
       if(mod(neg,2)/=0) write(*,*)'dkmxhf : uneven no. of negative pivots',&
&                          ' - log determinant not correct !'
   end if
   return
 
!  matrix not of full rank, return generalised inverse
11 nrank=ii-1; ij=0
   do i=1,N
      if(iflag(i) == 0)then ! ... zero out row/column
         w(i)=0.d0
         do j=i,n
            ij=ij+1;  a(ij)=0.d0
         end do
      else
         ij=ij+1; ww=-w(i); w(i)=ww; a(ij)=ww
         do j=i+1,n
            ij=ij+1
            if(iflag(j).ne.0)then
               a(ij)=-a(ij)
            else
               a(ij)=0.d0
            end if
         end do
      end if
  end do
  if(iopt) write(*,'(a,i5,a,i5)')'Generalised inverse of matrix with order =', &
&           n,'   and rank =',nrank
  go to 300
  end subroutine dkmxhf