function [ l, u ] = r8mat_lu ( a )

%*****************************************************************************80
%
%% r8mat_lu() computes the LU factorization of an R8MAT.
%
%  Discussion:
%
%    The routine is given an M by N matrix A, and produces
%
%      L, an M by M unit lower triangular matrix, and
%      U, an M by N upper triangular matrix
%
%    so that
%
%      A = L * U.
%
%    This procedure will fail if it encounters a pivot entry that is
%    exactly zero; it will perform badly if it encounters a pivot entry
%    that is very small and affected by roundoff error.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 November 2021
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(M,N), the M by N matrix to be factored.
%
%  Output:
%
%    real L(M,M), the M by M unit lower triangular factor.
%
%    real U(M,N), the M by N upper triangular factor.
%
  [ m, n ] = size ( a );
%
%  Initialize:
%
%    U:=A
%    L:=Identity
%    P:=Identity
%
  u(1:m,1:n) = a(1:m,1:n);

  l = r8mat_identity ( m, m );
%
%  On step J, find the pivot row, IPIV, and the pivot value PIVOT.
%
  for j = 1 : min ( m - 1, n )

    if ( u(j,j) == 0.0 )
      error ( 'r8mat_lu(): Fatal error, zero pivot.' );
    end
%
%  Zero out the entries in column J, from row J+1 to M.
%
    for i = j + 1 : m

      if ( u(i,j) ~= 0.0 )
        l(i,j) = u(i,j) / u(j,j);
        u(i,j) = 0.0;
        u(i,j+1:n) = u(i,j+1:n) - l(i,j) * u(j,j+1:n);
      end

    end

  end

  return
end
