function [ P, L, U ] = gauss_plu ( A, wait )

%*****************************************************************************80
%
%% gauss_plu() uses Gauss eliminiation to find the PLU factors of a matrix.
%
%  Discussion:
%
%    The desired factorization is A = P' * L * U
%    * P is a permutation matrix;
%    * L is a unit lower triangular matrix;
%    * U is an upper triangular matrix.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    07 January 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A(m,n), the system matrix.
%
%    logical wait:
%    0: no extra printout
%    1: extra printout, but no pause.
%    2: extra printout, and pausing.
%
%  Output:
%
%    global real P(m,m), L(m,m), U(m,n), the PLU factors.
%
  if ( nargin < 2 )
    wait = 0;
  end

  [ m, n ] = size ( A );
%
%  PLU initialize.
%
  P = eye(m);
  L = eye(m);
  U = A;

  info = 0;

  for j = 1 : n - 1
%
%  Choose the pivot row P for variable J.
%
    [ v, p ] = max ( abs ( U(j:n,j) ) );
    p = p + j - 1;

    if ( v == 0.0 )
      continue;
    end
%
%  Swap rows P and J.
%
    if ( p ~= j )
      T      = P(j,:);
      P(j,:) = P(p,:);
      P(p,:) = T;

      T      = U(j,:);
      U(j,:) = U(p,:);
      U(p,:) = T;

      T          = L(j,1:j-1);
      L(j,1:j-1) = L(p,1:j-1);
      L(p,1:j-1) = T;
    end
%
%  Eliminate U(i,j), and store multiplier in L(i,j).
%
    for i = j + 1 : n
      s = U(i,j) / U(j,j);
      U(i,:) = U(i,:) - s * U(j,:);
      U(i,j) = 0.0;
      L(i,j) = s;
    end
%
%  Show that it's still true that P' * L * U = A.
%
    if ( 0 < wait )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  After processing column %d:\n', j );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  P:\n' );
      fprintf ( 1, '\n' );
      disp ( P );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  L:\n' );
      fprintf ( 1, '\n' );
      disp ( L );
      fprintf ( 1, '\n' );
      fprintf ( 1, '  U:\n' );
      fprintf ( 1, '\n' );
      disp ( U );
      fprintf ( 1, '  ||A-P''*L*U|| = %g\n', norm ( A - P' * L * U ) );
      if ( 1 < wait )
        pause ( );
      end
    end

  end

  return
end

