% Examples of using the LU decomposition and Notes about the LU
% decomposition.
% (Code by Ben Walter, METU-NCC, benjamin@metu.edu.tr)
%
%% EXAMPLES of using the LU decomposition to solve problems
%
% We will give three examples of using the LU decomposition to solve
% the equation Ax=b in matlab. 
%  Example 1 uses matlab like a calculator: without for loops or any 
%   programming. 
%  Example 2 uses matlab like a simple programming language: computing the 
%   answer with nested for loops.
%  Example 3 takes advantage of matlab's ability to do matrix computations
%   quickly: it computes the answer with a single for loop and a matrix
%   multiplication operation.
%
%% EXAMPLE 1: solving Ax=b using LU decomposition (by hand)
% Once you know the LU decomposition of A, solving Ax=b may be done by 
% solving LUx=b.  The benefit of this is that solving Tx=b where T is 
% triangular is easy -- just use back-substitution.  Since we have two 
% triangular matrices, our solution will take two steps:  
%  First we deal with the outer matrix -- solve Ly=b
%  Then we deal with the inner matrix -- solve Ux=y
% The x we find will make Ax = LUx = L(Ux) = Ly = b.
%
clc

fprintf('Example 1: LU decomposition to solve Ax=b (by hand).')
A = [2 2 3; 2 3 4; 3 2 6]
b = [2; 3; 1]

%%
% Now we'll get the LU decomposition using our last algorithm from the 
% previous part -- we'll add in semicolons so that this runs silently
%
fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')
%% 
% Now we'll solve Ly=b. Let's write the augmented matrix L|b so our work
% is more clear
%
fprintf('----------\nThe augmented matrix L|b:')
Lb = [L,b]
clear y    % this command ensures that y isn't already storing data

%%
% Step 1: Row 1 of Ly=b tells us 1*y(1) + 0*y(2) + 0*y(3) = 2
%
fprintf('-----\n From row 1:   1*y(1)  +   0*y(2) + 0*y(3) = 2.')

y(1) = 2

%%
% Step 2: Row 2 of Ly=b tells us 1*y(1) + 1*y(2) + 0*y(3) = 3
%
fprintf('-----\n From row 2:   1*y(1)  +   1*y(2) + 0*y(3) = 3.')

y(2) = 3 - 1*y(1)

%%
% Step 3: Row 3 of Ly=b tells us 3/2*y(1) + (-1)*y(2) + 1*y(3) = 1
%
fprintf('-----\n From row 3: 3/2*y(1) + (-1)*y(2) + 1*y(3) = 1.')

y(3) = 1 - 3/2*y(1) - (-1)*y(2)

%%
% We've now computed y so that Ly = b -- let's check our work
%
fprintf('-----\nDoes Ly = b?')
y = y';    % we made y be a row vector, but it should be a column vector
Ly = L*y
b

%%
% Now we can solve Ux=y. Let's write the augmented matrix U|y so our work
% is more clear 
%
fprintf('----------\nThe augmented matrix U|y:')
Uy = [U,y]
clear x     % this command ensures that x isn't already storing data

%%
% Step 1: Row 3 of Ux=y tells us 0*x(1) + 0*x(2) + 5/3*x(3) = -1
%
fprintf('-----\n From row 3: 0*x(1) + 0*x(2) + 5/3*x(3) = -1.')

x(3) = (-1)/(5/2)

%%
% Step 2: Row 2 of Ux=y tells us 0*x(1) + 1*x(2) + 1*x(3) = 1
%
fprintf('-----\n From row 2: 0*x(1) + 1*x(2) + 1*x(3) = 1.')

x(2) = ( 1 - 1*x(3) )/1

%%
% Step 3: Row 1 of Ux=y says 2*x(1) + 2*x(2) + 3*x(3) = 2
%
fprintf('-----\n From row 1: 2*x(1) + 2*x(2) + 3*x(3) = 2.')

x(1) = ( 2 - 2*x(2) - 3*x(3) )/2

%% 
% We've now computed x so that Ux=y, this should also make
% Ax = LUx = Ly = b.  Let's check our work:
%
fprintf('-----\nDoes Ax = b?')
x = x';  % we made x a row vector, but it should be a column vector
Ax = A*x
b





%% EXAMPLE 2: Solving Ax=b using LU decomposition (using two for loops)
% We'll automate the solving process using some simple programming.
%
clc

fprintf('Example 2: LU decomposition to solve Ax=b (with for loops).')
A = [3 1 2; 2 1 3; 3 3 3]
b = [1; -2; 2]

%%
% Get LU decomposition
%
fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')
%%
% We'll compute y with a for loop: row i of Ly=b says  
%  b(i) = L(i,1)*y(1) + ... + L(i,i-1)*y(i-1) + 1*y(i) + 0*y(i+1) + 0*... 
% solving for y(i) we get:
%  y(i) = b(i) - L(i,1)*y(1) - L(i,2)*y(2) - ... - L(i,i-1)*y(i-1)
%
y=0;
fprintf('-----\nComputing y:')

for (i=1:numrows) 
    y(i) = b(i);
    for (j=1:i-1)
        y(i) = y(i) - L(i,j)*y(j);
    end
end

y = y'

%%
% Check our answer
%
fprintf('-----\nDoes Ly = b?')

Ly = L*y
b

%% 
% We'll compute x with a for loop: row i of Ux=y says
%  y(i) = 0*x(1) + ... + 0*x(i-1) + U(i,i)*x(i) + U(i,i+1)*x(i+1) + ...
% solving for x(i) we get:
%  x(i) = ( y(i) - U(i,i+1)*x(i+1) - U(i,i+2)*x(i+2) - ... )/U(i,i)
%
% Note that our loop needs to run backwards -- from 3 to 1
%
x=0;
fprintf('-----\nComputing x:')

for (i=numcols+1-(1:numcols)) % 4 - [1 2 3] = [(4-1) (4-2) (4-3)] = [3 2 1]
    x(i) = y(i);
    for (j=i+1:numcols)
        x(i) = x(i) - U(i,j)*x(j);
    end
    x(i) = x(i)/U(i,i);
end

x = x'

%%
% Check our answer
%
fprintf('-----\nDoes Ax = b?')

Ax = A*x
b





%% EXAMPLE 3: Solving Ax=b using LU decomposition (using a single for loop)
% As usual, we can replace the inner 'for' loop by matrix multiplication.
%
clc

fprintf('Example 3: LU decomposition to solve Ax=b (vectorized).')
A = [2 1 3; 4 1 2; 1 2 3]
b = [0; -1; 3]

%%
% Get LU decomposition
%
fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')
%%
% Looking at row i of Ly=b gave us the formula:
%  y(i) = b(i) - L(i,1)*y(1) - L(i,2)*y(2) - ... - L(i,i-1)*y(i-1)
%       = b(i) - (row i of L from column 1 to column i-1)*(y from 1 to i-1)
%
y=0;
fprintf('-----\nComputing y:')

for (i=1:numrows) 
    y(i) = b(i) - L(i,1:(i-1))*y(1:(i-1))'; 
end

y = y'

%%
% Check our answer
%
fprintf('-----\nDoes Ly = b?')

Ly = L*y
b

%% 
% Looking at row i of Ux=y gave us the formula:
%  x(i) = ( y(i) - U(i,i+1)*x(i+1) - U(i,i+2)*x(i+2) - ... )/U(i,i)
%     = ( y(i) - (row i of U from i+1 to end)*(x from i+1 to end) )/U(i,i)
%
x=0;
fprintf('-----\nComputing x:')

for (i=numcols+1-(1:numcols))      % i=[... 4 3 2 1]
    x(i) = ( y(i) - U(i,i+1:end)*x(i+1:end)' )/U(i,i);
end

x = x'

%%
% Check our answer
%
fprintf('-----\nDoes Ax = b?')

Ax = A*x
b



%% EX 4: Time comparison of vectorized method
% Since we went to the trouble of figuring out a vectorized version of our
% algorithm to solve Ax=b using the LU decomposition.  Let's check how much
% faster it is!
%
clc

fprintf('Example 4:  Time comparison of methods!\n')
fprintf(' -->A will be a random 1000 x 1000 matrix\n')
fprintf(' -->b will be a random 1000 vector\n')

A = randn(1000);
b = randn(1000,1);

fprintf(' (be patient in the next step... LU decomposition should take a few seconds)\n')
%%
fprintf('Computing LU decomposition...')

[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')

%%
fprintf('---------\nTime for solution using nested for loops [min, sec]:\n')
start=now;                   % record the starting time

y=0;
for (i=1:numrows) 
    y(i) = b(i);
    for (j=1:i-1)
        y(i) = y(i) - L(i,j)*y(j);
    end
end
y = y';

x=0;
for (i=numcols+1-(1:numcols)) % 4 - [1 2 3] = [(4-1) (4-2) (4-3)] = [3 2 1]
    x(i) = y(i);
    for (j=i+1:numcols)
        x(i) = x(i) - U(i,j)*x(j);
    end
    x(i) = x(i)/U(i,i);
end
x = x';

time=[minute(now-start), second(now-start)]

%%
fprintf('---------\nTime for solution using vectorized algorithm [min, sec]:\n')
start=now;                   % record the starting time

y=0;
for (i=1:numrows) 
    y(i) = b(i) - L(i,1:(i-1))*y(1:(i-1))'; 
end
y = y';

x=0;
for (i=numcols+1-(1:numcols))      % i=[... 4 3 2 1]
    x(i) = ( y(i) - U(i,i+1:end)*x(i+1:end)' )/U(i,i);
end
x = x';

time=[minute(now-start), second(now-start)]



%% EX 5: Getting inv(A) from the LU decomposition
% To be added....  
%   NOTE: matlab computes inv(A) by first finding LU decomposition
%  

%% EX 6: Getting det(A) from the LU decomposition
% Once you know A = LU, you get det(A) = det(LU) = det(L)*det(U).
% L and U are both triangular, so their determinants are merely the
% products of their diagonal elements.  Thus:
%   det(A) = det(L)*det(U) = 1*U(1,1)*U(2,2)*U(3,3)
% If you check the matlab manual, it says this is how matlab computes det.
%
clc
fprintf('Calculating determinant using LU decomposition')
A = [2 1 3; 4 1 2; 1 2 3]

%%
% Get LU decomposition
%
fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')

%%
% Check determinant
%
fprintf('-----\nDeterminant should be product of diagonal elements of U')

U
MatlabDetA = det(A)
ComputedDetA = prod(diag(U))



%% NOTE 1: Missing pivots.
% One way of interpreting A=LU is to say U is ref(A) and L is the matrix of
% row operations.  This is only true if you reduce A without any scaling or
% row swap operations.  For many matrices this can be done, but sometimes
% there aren't enough pivots, or aren't pivots in the right places for this
% to work.
%
%% NOTE 1.1: LU decomposition fails due to required row swaps
% 
clc
fprintf('In this example, LU decomposition will FAIL due to required row swap.')

A = [0 1 2; 1 2 3; 2 3 5]

fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')

L
U

fprintf('The LU decomposition failed because we could not divide by A(1,1)=0.\n')
%%
% If we swap rows 1 and 2 then it will now work
%
fprintf('-----\nIf we swap rows 1 and 2, we can get an LU decomposition:\n')
PA = A([2 1 3],:)

fprintf('Computing LU decomposition...')
[numrows,numcols] = size(PA);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = PA(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( PA((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')

L
U

%%
% This L and U is not a decomposition of the original matrix A, it is a
% decomposition of A with the first two rows swapped.  In matrix language,
% it is a decomposition of PA where P is the "permutation matrix" which
% swaps rows 1 and 2:
%
fprintf('-----\nThe permutation matrix which swaps rows 1 and 2 is:')
P = [0 1 0; 1 0 0; 0 0 1] 

%%
% You can make P by operating on eye(3) with the row swap you want to encode
%
fprintf('-----\nWe can build P by doing this row operation to the identity matrix')
P = eye(3)
P = P([2 1 3],:)

%%
% Check our work:

fprintf('-----\nCheck: P*A = L*U')
PA = P*A
LU = L*U

fprintf('--->Now let''s try to add row swaps to our LU decomposition code...\n')

%% LU decomposition with simple partial pivoting
% Clearly we need to introduce some sort of partial pivoting into our code.
% It is unclear how to add partial pivoting into Doolittle's method 
%(algorithm 2), so we'll use the LU decomposition method from algorithm 1.3.
% We can just copy the method directly, adding row swaps on L and U into 
% the middle.

fprintf('Computing LU decomposition (with simple partial pivoting)...')
[numrows,numcols] = size(A);
L = zeros(numrows); 
U = A;

P = eye(numrows);   % This matrix will record row swaps

for (j=1:numcols) 
	[maxvalue,maxindex] = max(abs(U(j:end,j)));  % find the pivot for column j
	maxindex = maxindex+j-1;                   % correct maxindex to be index in U
	P([j maxindex],:) = P([maxindex j],:);     % record row swap
	U([j maxindex],:) = U([maxindex j],:);     % do row swap
	L([j maxindex],:) = L([maxindex j],:);     % do row swap
	L((j+1):end,j) = U((j+1):end,j)/U(j,j);              % record row ops
	U((j+1):end,:) = U((j+1):end,:) - L((j+1):end,j)*U(j,:); % do row ops
end

L = L + eye(numrows);

fprintf('finished.\n')


%%
% Check our partial pivoting code
%

fprintf('-----\nCheck LU decomposition with partial pivoting: PA = LU?')
P
A
PA = P*A
L
U
LU = L*U

%%
% This is how matlab does LU decompositions:
%
fprintf('-----\nresults of the matlab command lu(A):')
[L,U,P] = lu(A)




%% NOTE 1.2: LU decomposition fails due to no pivot
% Adding in row swaps helps out a lot, but it doesn't solve all problems.
% Sometimes a column doesn't have a pivot anywhere, so row swaps won't save
% you.

clc
fprintf('In this example, LU decomposition will FAIL due lack of rank.')

A = [1 2 3; 1 2 3; 1 2 3]


fprintf('Computing LU decomposition (with partial pivoting)...')
[numrows,numcols] = size(A);
L = zeros(numrows); 
U = A;
P = eye(numrows);  
for (j=1:numcols) 
	[maxvalue,maxindex] = max(abs(U(j:end,j)));  % find the pivot for column j
	maxindex = maxindex+j-1;                   % correct maxindex to be index in U
	P([j maxindex],:) = P([maxindex j],:);     % record row swap
	U([j maxindex],:) = U([maxindex j],:);     % do row swap
	L([j maxindex],:) = L([maxindex j],:);     % do row swap
	L((j+1):end,j) = U((j+1):end,j)/U(j,j);    % record row ops
                                               % do row ops:
	U((j+1):end,j:end) = U((j+1):end,j:end) - L((j+1):end,j)*U(j,j:end);
end
L = L + eye(numrows);
fprintf('finished.\n')

L
U

fprintf(' The LU decomposition failed because even partial pivoting couldn''t\n');
fprintf(' find a pivot in the second column.\n-----\n')

%%  LU decomposition with good partial pivoting
% After reducing the first column, everything left was zero.  We should not
% have continued trying to reduce, since that just introduced NaN in row 3.
% We can fix this by modifying our code so that it skips columns with no
% pivot.

fprintf('Computing LU decomposition (with good partial pivoting)...')
[numrows,numcols] = size(A);
L = zeros(numrows); 
U = A;
P = eye(numrows);

i=1; j=1;    % i,j will be the row,col of the next pivot
while ((j<=numcols) & (i<numrows))
	[maxvalue,maxindex] = max(abs(U(i:end,j))); % find the pivot for column j
	if (maxvalue ~= 0)
		maxindex = maxindex+i-1;        % correct maxindex to be index in U
		P([i maxindex],:) = P([maxindex i],:);     % record row swap
		U([i maxindex],:) = U([maxindex i],:);     % do row swap
		L([i maxindex],:) = L([maxindex i],:);     % do row swap
		L((i+1):end,j) = U((i+1):end,j)/U(i,j);    % record row ops
                                                   % do row ops:
		U((i+1):end,j:end) = U((i+1):end,j:end) - L((i+1):end,j)*U(i,j:end); 
        i=i+1;                          % next pivot is one row down
    end
    j=j+1;                              % next pivot is one column over
end
L = L + eye(numrows);
fprintf('finished.\n')

%%
% Check our work:

fprintf('-----\nCheck our work PA = LU?')
PA = P*A
L
U
LU = L*U




%% NOTE 2: rectangular matrices
% So far we've been dealing only with square matrices, but there's no
% reason not to try this with rectangular matrices.
%
clc
fprintf('Doolittle''s method already works if #cols > #rows:')

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


fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows) 
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')

L
U
LU = L*U

%%
% 
%
clc
fprintf('The row reduction method already works if #rows > #cols:')

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

fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = A;
for (j = 1:numcols) 
    L((j+1):end,j) = U((j+1):end,j)/U(j,j); 
	U((j+1):end,j:end) = U((j+1):end,j:end) - L((j+1):end,j)*U(j,j:end); 
end
fprintf('finished.\n')

L
U
LU = L*U

%%
%  If we are more careful with our indexing, Doolittle can work for all 
%  rectangular matrices 
%
clc
fprintf('If #cols < #rows, Doolittle''s method FAILS (unless we modify it)')

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

fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);

for (i=1:numrows)  % TO FIX, CHANGE THIS TO BE  "i=1:min(numrows,numcols)"
	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end);
	L((i+1):end,i) = ( A((i+1):end,i) - L((i+1):end,1:(i-1))*U(1:(i-1),i) )/U(i,i);
end
fprintf('finished.\n')

L
U
LU = L*U

%%
% If we are more careful with our indexing, the row reduction method can
% also work for all rectangular matrices.
%
clc
fprintf('If #rows < #cols, the row reduction method FAILS (unless we modify it)')

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

fprintf('Computing LU decomposition...')
[numrows,numcols] = size(A);
L = eye(numrows);
U = A;
for (j = 1:numcols) % TO FIX, CHANGE THIS TO BE  "j=1:min(numcols,numrows)"
    L((j+1):end,j) = U((j+1):end,j)/U(j,j); 
	U((j+1):end,j:end) = U((j+1):end,j:end) - L((j+1):end,j)*U(j,j:end); 
end
fprintf('finished.\n')

L
U
LU = L*U
