subroutine gauss_lu(A)
! performs an LU factorization of the matrix A using
! Gaussian reduction.
! A is the matrix to be factored.
! L is a lower triangular factor with 1's on the diagonal
! U is an upper triangular factor.
! A = L * U
! L is stored as the lower triangle without diagonal of A, U is
!      stored as the upper triangle of A
! from the Matlab
! will abort if pivoting is needed

! mms 3/4/13
real(PREC) :: A(:,:)
integer :: n,Irow,Jcol,i,i4max(1)
real(PREC) :: mult


n=size(A,1);

do Jcol=1,n
  ! check if pivoting needed.  Hope not!
  i4max = maxloc ( abs ( A(Jcol:n, Jcol ) ) );
  if (i4max(1) /= 1) then
    if (abs(A(Jcol-1+i4max(1),Jcol)) > abs(A(Jcol,Jcol))*10._PREC) &
      & stop 'gauss_lu: Pivoting needed!'
  endif

  do Irow=Jcol+1,n
    ! compute Irow multiplier
    mult=A(Irow,Jcol)/A(Jcol,Jcol);
  
    ! multiply row "Jcol" by L(Irow,Jcol) and subtract from row "Irow"
    A(Irow,Jcol:n)=A(Irow,Jcol:n)-mult*A(Jcol,Jcol:n);

    ! put multiplier in A(Irow,Jcol)
    A(Irow,Jcol)=mult;
  enddo
enddo

end subroutine gauss_lu