% LU FACTORIZATION EXAMPLES
%   (Code by Ben Walter, METU-NCC, benjamin@metu.edu.tr)
%
% In this lecture, we give a series of different algorithms for computing
% the LU factorization of a matrix.  Our goal is to write A=LU where 
% L is lower triangular and U is upper triangular.  In the next section
% we will discuss what we can do with such a factorization, for now
% we are concerned only with how to find it.
%
% We give two different kinds of algorithm -- one based on reducing A to
% get U, storing the row operations at each step in L; and the other based
% on directly solving A=LU to get L and U.  Our algorithms are:
%
%   ALGORITHM 1.0: Using row reduction by hand
%   ALGORITHM 1.1: Using row reduction by hand with one for loop
%   ALGORITHM 1.2: Using row reduction and nested for loops
%   ALGORITHM 1.3: Using row reduction and a loop with matrix
%     multiplication
%
%   ALGORITHM 2.0: Direct solution (Doolittle's method) by hand
%   ALGORITHM 2.1: Direct solution (Doolittle) using nested for loops
%   ALGORITHM 2.2: Direct solution (Doolittle) using a single for loop and 
%     matrix ops
%
% Doolittle's method is faster than row reduction, because the row 
% reduction method will sometimes require you to do operations like 
%    a*x + b*x 
% In Doolittle's method you will instead do the operation
%    (a+b)*x
% which takes one less multiplication.  The fastest algorithm of all is
% 2.2 where we replace the inner for loop by clever matrix multiplication. 
% Matlab is built to do matrix multiplication very quickly, so using matrix
% multiplication is faster than using a for loop.  The difference in speed 
% is quite noticible already for 1000 x 1000 matrices. 
%
%% ALGORITHM 1.0:  Using row reduction, by hand
%
clc    % this command clears the command window, so we don't get confused
fprintf('LU decomposition via row reduction (by hand)\n')   % print description

%%
% Enter Matrix
%
A = [2 2 3; 2 3 4; 3 2 6]

%%
% Prepare L for filling with entries
% We start with a 3x3 identity matrix for L, we will fill this with the
% rest of L's entries as we do our row reduction
%
fprintf('-----\nGet the matrix L ready to be filled with row operations')

L = eye(3)   

%%
% U will be ref(A)
% We begin with U = A and will reduce step by step, recording each row
% operation in the appropriate position of L
%
fprintf('-----\nGet the matrix U ready to be reduced to upper triangular')

U = A

%%
% Do first reduction step with U
% replace row 2 by row 2 - row 1
%
fprintf(' ---> Step 1:  (row 2) = (row 2) - (row 1) in U')

U(2,:) = U(2,:) - 1*U(1,:)

%%
% Record the reduction step in L
%
fprintf('      Record (row 2) = (row 2) - 1*(row 1) in L(2,1)')

L(2,1) = 1

%%
% Check that L recorded the row operation correctly:
%  LU should be the same as A
%
fprintf('-----\nCheck our work so far: Does A = LU?')

A
LU = L*U

%%
% Do second reduction step with U
% replace row 3 by row 3 - 3/2 row 1
%
fprintf(' ---> Step 2:  (row 3) = (row 3) - (3/2)*(row 1) in U')

U(3,:) = U(3,:) - 3/2*U(1,:)

%%
% Record the reduction step in L
%
fprintf('      Record (row 3) = (row 3) - (3/2)*(row 1) in L(3,1)')
L(3,1) = 3/2

%%
% Check that L recorded the row operation correctly:
%  LU should be the same as A
%
fprintf('-----\nCheck our work so far: Does A = LU?')

A
LU = L*U

%%
% Do third reduction step with U
% replace row 3 by row 3 + row 2
%
fprintf(' ---> Step 3:  (row 3) = (row 3) + (row 2)')

U(3,:) = U(3,:) + U(2,:)

%%
% Record the reduction step in L
%
fprintf('      Record: (row 3) = (row 3) - (-1)*(row 2) in L(3,2)')

L(3,2) = -1

%%
% LU factorization is complete!  Check answer
% LU should be the same as A
%
fprintf('-----\nFinished!  Check our work:  Does A = LU?')

A
LU = L*U




%% ALGORITHM 1.1: Using row reduction by hand with one for loop 
%  
clc

fprintf('LU decomposition via row reduction with a for loop')

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

%%
%  Set things up
%
fprintf('Set up L and U:')

L = eye(3)
U = A

%%
% Our for loop will reduce the first column
%  we will loop over the rows that need to be reduced (rows 2 and 3)
%
fprintf(' --->Reduce the first column:')

for (i=2:3)
	L(i,1) = U(i,1)/U(1,1)   % record the row op we'll do
	U(i,:) = U(i,:) - U(1,:)*L(i,1) % do the row op
end

%%
% Now finish the reduction by taking care of the second column
%
fprintf(' --->Reduce the second column:')

L(3,2) = U(3,2)/U(2,2)     % record the row op we'll do
U(3,:) = U(3,:) - U(2,:)*L(3,2)   % do the row op

%%
% Check the factorization: LU should equal A 
%
fprintf('-----\nFinished!  Check our work: Does A = LU?')

LU = L*U
A




%% ALGORITHM 1.2: Using row reduction and nested for loops
%  
clc

fprintf('LU decomposition via row reduction with two for loops')

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

%%
% We will set things up so these commands will work for bigger matrices
%
fprintf('Set up L and U:')

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

%%
% Look one column at a time, reducing all rows below row j.
%
fprintf('-----\nDoing row reduction...\n')
for (j = 1:numcols)        % j will be column we're reducing 
	for (i = j+1:numrows)    % i will be row we're reducing 
        
        % Please ignore the following fprintf statements
        fprintf(' ---> killing element %i,%i with the row operation:\n',[i,j])
        fprintf('       (row %i) = (row %i) - %g*(row %i)', [i,i,U(i,j)/U(j,j),j])

        L(i,j) = U(i,j)/U(j,j) % record the row op we'll do
		U(i,:) = U(i,:) - U(j,:)*L(i,j) % do the row op
	end
end
fprintf('Done!\n')

%%
% Check the answer:  LU == A?
%
fprintf('-----\nCheck our work:  Does A = LU?')

LU = L*U
A



%% ALGORITHM 1.3: Using row reduction and a loop with matrix multiplication
% We can speed up our script considerably by removing the inner for loop
% and replacing it by matrix multiplication.  It is a little more
% complicated to work out, but it runs much faster.
%  
clc

fprintf('LU decomposition via row reduction with one for loop and matrix mult')

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

%%
% We will set things up so it will work for bigger matrices
%
fprintf('Set up L and U:')

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

%%
% We will reduce the matrix U one column at a time, but now we will use 
% matrix multiplication instead of a second for loop
%  NOTE: matrix multiplication is faster than for loops, because matlab is 
%   specificially programmed to do it really, really fast.
%
fprintf('-----\nDoing row reduction...\n')

for (j = 1:numcols-1)        % j will be column we're reducing 
    
% instead of inserting an inner for loop which runs over all rows below j
% we will reduce all these rows at once, using a matrix multiplication

    fprintf(' --->Reducing column %i', j) 

% The column vector L(:,j) stores the row ops we will do on all rows of U
% when reducing column j of U (note that we only change rows below row j)
    L((j+1):end,j) = U((j+1):end,j)/U(j,j) 
    
% We can use matrix multiplication to do all our row ops at once!
% Notice that we multiply L(:,j)*U(j,:) (column * row)
%  the result of this multiplication is many copies of U(j,:) each scaled
%  by the appropriate L(i,j) element.
	U((j+1):end,j:end) = U((j+1):end,j:end) - L((j+1):end,j)*U(j,j:end) 
    
end

fprintf('Done!\n')

%%
% Check the answer:  LU == A?
%
fprintf('-----\nCheck our work:  Does A = LU?')
LU = L*U
A





%% ALGORITHM 2.0: Direct solution by hand
% Instead of finding L and U using row reduction, we can find them by 
% directly solving A = LU, using the fact that L will be lower triangular
% with 1's on the diagonal, and U will be upper triangular
%
clc

fprintf('LU decomposition using Doolittle''s method by hand\n')

%%
% Enter matrix, set up L and U
%
fprintf('Set up A, L, and U:')

A = [2 2 3; 2 3 4; 3 2 6]
L = eye(3)
U = zeros(3) 

%%
% We get row 1 of U by multiplying (columns of U) * (row 1 of L)
% row 1 of L is [1 0 0] so this gives: 
%   A(1,j) = 1*U(1,j) + 0*U(2,j) + 0*U(3,j)  <-- solve for U(1,j)
% (In other words Row 1 of U is Row 1 of A)
%
fprintf('-----\nFill row 1 of U:  (row 1 of U) is (row 1 of A)')

U(1,:) = A(1,:)

%%
% We get column 1 of L by multiplying (rows of L) * (column 1 of U)
% column 1 of U is [A(1,1); 0; 0] so this gives:
%   A(1,j) = L(1,j)*A(1,1) + L(2,j)*0 + L(3,j)*0  <-- solve for L(i,j)
% (In other words, Column 1 of L is (Column 1 of A)/A(1,1))
%
fprintf('-----\nFill column 1 of L: (column 1 of L) is (column 1 of A)/A(1,1)')

L(:,1) = A(:,1)/A(1,1)

%%
% Check our work:  
%  At this point LU should correctly give the first row and column of A
%
fprintf('-----\nCheck work so far:  LU should be correct on first row and column')

A
LU = L*U

%%
% Fill row 2 of U:
% We figure out these elements by multiplying (row 2 of L) * (cols of U)
% row 2 of L is [L(2,1) 1 0]

fprintf('-----\nFill row 2 of U')

% multiplying (row 2 of L) * (column 2 of U) gives: 
%  A(2,2) = L(2,1)*U(1,2) + 1*U(2,2) + 0*U(3,2) <-- solve for U(2,2)
U(2,2) = A(2,2) - L(2,1)*U(1,2)

% multiplying (row 2 of L) * (column 3 of U) gives:
%  A(2,3) = L(2,1)*U(1,3) + 1*U(2,3) + 0*U(3,3) <-- solve for U(2,3)
U(2,3) = A(2,3) - L(2,1)*U(1,3)

%%
% Fill column 2 of L:
% We figure out these elements by multiplying (rows of L) * (col 2 of U)
% col 2 of U is [U(2,1); U(2,2); 0]

fprintf('-----\nFill column 2 of L')

% multiplying (row 3 of L) * (column 2 of U) gives:
%  A(3,2) = L(3,1)*U(1,2) + L(3,2)*U(2,2) + L(3,3)*0 <-- solve for L(3,2)
L(3,2) = ( A(3,2) - L(3,1)*U(1,2) )/U(2,2)

%%
% Check our work:
%  At this point LU should correctly give the first two rows and cols of A
%
fprintf('-----\nCheck work so far:  LU should be correct on first two rows and columns of A')

A
LU=L*U

%%
% Fill row 3 of U:
% We can figure these out by multiplying (row 3 of L)* (columns of U)
% (row 3 of L is [L(3,1) L(3,2) 1])

fprintf('-----\nFill row 3 of U')

% Multiplying (row 3 of L) * (column 3 of U) gives:
%  A(3,3) = L(3,1)*U(1,3) + L(3,2)*U(2,3) + 1*U(3,3) <-- solve for U(3,3)
U(3,3) = A(3,3) - ( L(3,1)*U(1,3) + L(3,2)*U(2,3) )

%%
% Check A = LU for all rows and columns
%
fprintf('-----\nDone!  Check our work: Does A = LU?')
A
LU = L*U



%% ALGORITHM 2.1: Direct solution using nested for loops
% Our previous work required us to do a number of solving steps by hand 
% and then tell matlab specificially what operations to do.  Now we'll
% automate the process.
%
clc

fprintf('LU decomposition using Doolittles method with for loops')

A = [2 2 3; 2 3 4; 3 2 6]

%%
% We'll make this work for nxn matrices
%
fprintf('-----\nSet up L and U')

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

%%
%
fprintf('-----\nComputing LU decomposition...\n')

for (i=1:numrows)  % i will be the row of U/column of L we are filling 

% First we fill row i of U from column i to the end of the row 
%  (the other elements in row i of U are 0)
% We get U(i,j) by multiplying (row i of L) * (column j of U)

    fprintf(' --->Filling row %i of U', i)

	for (j=i:numcols)  % j is the col of the element we're adding to U     
		U(i,j) = A(i,j) - L(i,1:(i-1))*U(1:(i-1),j)
    end

    
% Now we fill column i of L from row i+1 to the end of the column
%  (the other elements of column i of L are either 0 or 1 (the diagonal))
% We get L(j,i) by miltiplying (row j of L) * (column i of U)

    fprintf(' --->Filling column %i of L', i)

	for (j=(i+1):numrows) % j is the row of the element we're adding to L
		L(j,i) = ( A(j,i) - L(j,1:(i-1))*U(1:(i-1),i) )/U(i,i)
	end

end

fprintf('\nDone!\n')

%%
% Test that we succeeded
%
fprintf('-----\nCheck our work: Does A = LU?')
A
LU = L*U




%% ALGORITHM 2.2: Direct solution using a single for loop and matrix ops
% We can speed up our script considerably by removing the inner for loop
% and replacing it by matrix multiplication.  It is a little more
% complicated to work out, but it runs much faster.
%
clc

fprintf('LU decomposition using Doolittles method with loop and matrix mult')

A = [2 2 3; 2 3 4; 3 2 6]

%%
%
fprintf('-----\nSet up L and U')

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

%%
%
fprintf('-----\nComputing LU decomposition...\n')

for (i=1:numrows)  % i will be the row of U/col of L we are filling 

% instead of using a second loop, we will cleverly use matrix operations
% to fill the entire row with one command:
%  to get row U(i,:) we want to multiply (row i of L) * (columns of U)
%  we will stick all of the columns of U together and multiply them by 
%  row i at once.  This is faster than using a for loop, because matlab is 
%  optimized to do matrix multiplication quickly.

    fprintf(' --->Filling row %i of U', i)

	U(i,i:end) = A(i,i:end) - L(i,1:(i-1))*U(1:(i-1),i:end)
 
% similarly we will fill the entire column of L with one command:
%  instead of multiplying (row i of U) * (columns of L) one at a time, we
%  gather the columns of L together and do a single matrix multiplication

    fprintf(' --->Filling column %i of L', i)

	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('\nDone!\n')

%%
% test that we succeeded
%
fprintf('-----\nCheck our work:  Does A = LU?')

A
LU = L*U




%% TIME COMPARISON
% Let's compare how long it takes to run these different algorithms.
%
clc
fprintf('Time comparison of LU decomposition methods.\n')

fprintf('   A is a random 1000x1000 matrix\n')
A = randn(1000);   % A is a random 1000x1000 matrix

%%

fprintf('Time for LU decomposition by reduction with nested for loops [min, sec]:')
start=now;                   % record the starting time
[numrows,numcols] = size(A);
L = eye(numrows);
U = A;
for(j = 1:numcols)       
	for (i = j+1:numrows)   
		L(i,j) = U(i,j)/U(j,j); 
		U(i,j:end) = U(i,j:end) - U(j,j:end)*L(i,j); 
	end
end
time=[minute(now-start), second(now-start)]

%%

fprintf('Time for LU decomposition by reduction with matrix mult [min, sec]:')
start=now;                  % record starting time
[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
time=[minute(now-start), second(now-start)]

%%

fprintf('Time for LU decomposition by Doolittle with nested for loops [min, sec]:')
start=now;
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows); 
for (i=1:numrows)
	for (j=i:numcols)     
		U(i,j) = A(i,j) - L(i,1:(i-1))*U(1:(i-1),j);
    end
	for (j=(i+1):numrows)
		L(j,i) = ( A(j,i) - L(j,1:(i-1))*U(1:(i-1),i) )/U(i,i);
    end
end
time=[minute(now-start), second(now-start)]

%%

fprintf('Time for LU decomposition by Doolittle with matrix mult [min, sec]:')
start=now;
[numrows,numcols] = size(A);
L = eye(numrows);
U = zeros(numrows,numcols);
for (i=1:numrows)  
	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
time=[minute(now-start), second(now-start)]

fprintf('WOW!\n')
