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

dkmxhf.f90
```! ================================================
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
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```