
format compact
format short

clear
clc

% define L and b
L = [1,0;2,5]
b = [2;3]
% solve Ly=b
y = Lsolve(L,b)

disp('** check Lsolve, this should eval to zero')
ztest = L*y - b

% % test case: b has the wrong size
% bWrong = [2;3;4]
% Lsolve(L,bWrong)

% % test case: L is singular
% Lwrong = [1,0;2,0]
% Lsolve(Lwrong,b)

% % test case: L is not lower-triangular
% Lwrong = [1,2;2,5]
% y = Lsolve(Lwrong,b)
% disp('** check Lsolve, does this eval to zero?')
% ztest = Lwrong*y - b  % no it doesn't...
%                       % it only looks at the lower triangle of L

% % define U
% U = [2,1;0,3]
% % solve Ux=y, this is now the solution to LUx = b
% x = Usolve(U,y)
% 
% disp('** check Usolve, these should both eval to zero (up to roundoff)')
% ztest = U*x - y
% ztest = L*U*x - b
% 
% % do it in one step
% x = LUsolve(L,U,b)
% 
% disp('** check LUsolve, this should eval to zero (up to roundoff)')
% ztest = L*U*x - b
% 
% A = rand(5)
% [L,U] = LUdecomp(A)
% disp('** check LUdecomp, this should eval to zero (up to roundoff)')
% L*U-A
% 
% b = rand(5,1)
% x = Asolve(A,b)
% disp('** check Asolve, this should eval to zero (up to roundoff)')
% A*x - b

% 
% % this will throw an error...
% A = [0,1;1,0]
% [L,U] = LUdecomp(A)

% disp('** try it with a large matrix')
% nn=800
% A = rand(nn);
% b = zeros(nn);
% 
% tic
% [L,U] = LUdecomp(A);
% toc
% 
% disp('** top left block of L')
% disp( L( 1:5 , 1:5 ) )
% disp('** top left block of U')
% disp( U( 1:5 , 1:5 ) )
% checkLU = L*U;
% disp('** top left block of LU')
% disp( checkLU( 1:5 , 1:5 ) )
% disp('** top left block of A')
% disp( A( 1:5 , 1:5 ) )

