!........................................................................
!
! Subroutine :
!
! Description:
!
! Arguments  :          
!
! Mod/Commons:    
!
! Calls      :      
!
! Author     : 
!
! Modified   :     
!             
!........................................................................
!
! $Id$
!
! $Log$
! Revision 1.4  2004/08/09 13:41:21  dcs
! Minor change to reflect changes in precision organization.
!
! Revision 1.3  2004/03/08 14:03:43  cesrulib
! Update for recent DCS changes to recipes_f-90_LEPP
!
! Revision 1.2  2003/04/30 17:14:51  cesrulib
! dlr's changes since last import
!
! Revision 1.1.1.1  2002/12/13 19:23:28  cesrulib
! import bmadz
!
!
!........................................................................
!
  SUBROUTINE ludcmp(a,indx,d)
  USE precision_def, except => dp
  USE nrtype; USE nrutil, ONLY : assert_eq,imaxloc,nrerror,outerprod,swap
  IMPLICIT NONE
  REAL(rp), DIMENSION(:,:), INTENT(INOUT) :: a
  INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx
  REAL(DP), INTENT(OUT) :: d
  REAL(rp), DIMENSION(size(a,1)) :: vv
  REAL(DP), PARAMETER :: TINY=1.0e-20_dp
  INTEGER(I4B) :: j,n,imax
  n=assert_eq(size(a,1),size(a,2),size(indx),'ludcmp')
  d=1.0
  vv=maxval(abs(a),dim=2)
  if (any(vv == 0.0)) call nrerror('singular matrix in ludcmp')
  vv=1.0/vv
  do j=1,n
    imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j)))
    if (j /= imax) then
      call swap(a(imax,:),a(j,:))
      d=-d
      vv(imax)=vv(j)
    end if
    indx(j)=imax
    if (a(j,j) == 0.0) a(j,j)=TINY
    a(j+1:n,j)=a(j+1:n,j)/a(j,j)
    a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n))
  end do
  END SUBROUTINE ludcmp