%% LU Factorization - Derivation and Application
% Here we show the connection between Gaussian elimination and LU factorizaion. Then, we solve a linear system using A=LU, backward and forward subsitution to solve A\;\mathbf{x}=\mathbf{b}\ldotp
%
%% Gaussian Elimination Row Operations
% Begin with square matrix A.

clear
A=[1 1 0 3;2 1 -1 1;3 -1 -1 2;-1 2 3 -1]

%
% Use the basis of column vectors of the identity matrix.

E=eye(4)

%
% Multiply row 1 by A_{21}/A_{11}=2 and subtract from row 2. This is equivalent to the new matrix (I-2\mathbf{e}_2\mathbf{e}_1^T)A:

L=E-2*E(:,2)*E(:,1)'
U=L*A

%
% Then, multiply row 1 by A_{31}/A_{11}=3 and subtract from row 2. This gives (I-3\mathbf{e}_3\mathbf{e}_1^T)(I-2\mathbf{e}_2\mathbf{e}_1^T)A = (I-2\mathbf{e}_2\mathbf{e}_1^T-3\mathbf{e}_3\mathbf{e}_1^T)A = \left(I-\left(2\mathbf{e}_2+3\mathbf{e}_3\right)\mathbf{e}_1^T\right)A:

L=E-3*E(:,3)*E(:,1)'
U=L*U

%
% Finishing the column, (I+\mathbf{e}_4\mathbf{e}_1^T)(I-3\mathbf{e}_3\mathbf{e}_1^T)(I-2\mathbf{e}_2\mathbf{e}_1^T)A =  \left(I-\left(2\mathbf{e}_2+3\mathbf{e}_3 - \mathbf{e}_4\right)\mathbf{e}_1^T\right)A:

L1=(E-(-1)*E(:,4)*E(:,1)')*(E-3*E(:,3)*E(:,1)')*(E-2*E(:,2)*E(:,1)')
U=L1*A

%
% Letting L_{ij}=(I+\alpha_{ij}\mathbf{e}_i\mathbf{e}_j^T), then so far we have the matrix L_{41}L_{31}L_{21}A.Notice that L_1=L_{41}L_{31}L_{21} is lower triangular.
% Continuing with the remaining columns, we compute L_2=L_{43}L_{42}L_{32}L_{41}L_{31}L_{21}.
% [Note: For this problem \alpha_{43}=0, so L_{43}=Iand is not shown.]

L2=(E-(-3)*E(:,4)*E(:,2)')*(E-4*E(:,3)*E(:,2)')*L1
U=L2*A

%
% The final matrix, U=L_2A, is upper triangular. Solving forA=L_2^{-1}U. We note that L=L_2^{-1}is also lower triangular. So, we have factored A=\textrm{LU} into the product of upper and lower triangular matrices.
% The function (below) lufact carries out this procedure and produces the same results:

[L,U]=lufact(A)
L2^(-1)

%
%% Solving Linear Systems
%
% Now we solve A\;\mathbf{x}=\mathbf{b} for \mathbf{x} given square matrix A and constant vector \mathbf{x}. We do this in three steps:
% Factor A=L\;U\ldotp [Use lufact.]
% Solve L\;\mathbf{z}=\mathbf{b}.    [Use forwardsub.]
% Solve U\;\mathbf{x}=\mathbf{z}\ldotp    [Use backsub.]
% This is shown below using three functions lufact, forwardsub, backsub which need to either show below or are in m-files in the current directory.

A=[1 1 0 3;2 1 -1 1;3 -1 -1 2;-1 2 3 -1]
b=[4 1 -3 4]';

[L,U]=lufact(A);
z=forwardsub(L,b);
x=backsub(U,z)

%
%% Failure of LU Factorization
% Consider the system A\;\mathbf{x}=\mathbf{b} with

A = [2 0 4 3 ; -4 5 -7 -10 ; 1 15 2 -4.5 ; -2 0 2 -13]
b = [ 4; 9; 29; 40 ];

%

[L,U]=lufact(A);
z=forwardsub(L,b);
x=backsub(U,z)

%
% Switch rows 2 and 4 and apply LU factorization.

A([4 2],:) = A([2 4],:)
[L,U]=lufact(A)

%
% Something went wrong! To see what happened, we begin the Gaussian elimination process on column one.

for k=2:4
    alpha(k)=A(k,1)/A(1,1);
end
L1=(E-alpha(4)*E(:,4)*E(:,1)')*(E-alpha(3)*E(:,3)*E(:,1)')*(E-alpha(2)*E(:,2)*E(:,1)');
L1*A

%
% Note the zero in the A(2,2) position. We cannot continue. This pivot element can be changed with row pivoting. We swap rows 2 and 4 to obtain a nonzero pivot element. Then, we can solve the linear system.
%
% MATLAB and Julia have built-in functions lu() which implement the process of pivoting when needed if A is nonsingular. Overall, lu() finds L, U, P such that AP=LU\ldotp Then, to solve A\;\mathbf{x}=\mathbf{b}, we instead solve LU\;=PA\;\mathbf{x}=P\;\mathbf{b}.

b([4 2],:) = b([2 4],:);
[L,U,P]=lu(A);
z=forwardsub(L,P*b);
x=backsub(U,z)

%
%% LU Factorization routine.
%

function [L,U] = lufact(A)

n = length(A);
L = eye(n);   % ones on diagonal

% Gaussian elimination
for j = 1:n-1
  for i = j+1:n
    L(i,j) = A(i,j) / A(j,j);   % row multiplier
    A(i,j:n) = A(i,j:n) - L(i,j)*A(j,j:n);
  end
end
U=A;
%U = triu(A);  % If the output doesn't have clean zeros, triu cleans it up.
end