function lu_demo
% LU factorization demo.
% This is a function rather than a script so that the helper function at
% the end of the file is visible in both MATLAB and Octave.

% Define the matrix and row permutation.
A = [2, 1, 1, 0; 4, 3, 3, 1; 8, 7, 9, 5; 6, 7, 9, 8];

P = eye(4);
pivoting = true;
if pivoting
    row_order = [3, 4, 2, 1];
    P = P(row_order, :);
end
PA = P * A;

fprintf('Matrix A:\n')
disp(A)
fprintf('Permuted matrix P*A:\n')
disp(PA)

% Apply three elimination steps.
fprintf('Elimination matrices and intermediate products:\n')
M1 = mkM(4, -PA(2:4, 1) / PA(1, 1))
M1PA = M1 * PA

M2 = mkM(4, -M1PA(3:4, 2) / M1PA(2, 2))
M2M1PA = M2 * M1PA

M3 = mkM(4, -M2M1PA(4, 3) / M2M1PA(3, 3))
M3M2M1PA = M3 * M2M1PA

U = M3 * M2 * M1 * PA

fprintf('Each elimination matrix and its inverse:\n')
[M1, inv(M1)]
[M2, inv(M2)]
[M3, inv(M3)]

% Form L and check the factorization.
L = inv(M1) * inv(M2) * inv(M3)

fprintf('L*U and P*A (these should match):\n')
L * U
P * A

end

function M = mkM(n, m)
% Create an elimination matrix.

M = eye(n);
k = n - length(m);
M(k + 1:n, k) = m(:);
end
