%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% PA=LU decomposition and 2-step Gaussian elimination in MATLAB
% Application of Sherman-Morrison formula 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
i=1;
A = [2 4 -2; 4 9 -3; -2 -1 7]; b = [2 8 10]'; len=3;

sizes = 10:100:910; 
for i = sizes
        
A = rand(i); 
len = size(A,1); b = rand(len,1);
I = eye(len);
u = 2*I(:,1); v = I(:,2);
B = A - u*v';

%fprintf('\n Gaussian elimination for Ax=b\n');
% PA=LU, where P=permulation matrix 

[L,U,P]=lu(A);
z = L\(P*b);
x1 = U\z;


%fprintf('\n Gaussian elimination for Bx=b\n');
t0=clock;
x2 = B\b;
t1(i) = etime(clock,t0);

%fprintf('\n Sherman - Morrison formula result:\n'); 
% x2 = x1 + (v'*x1/(1-v'*y))*y

t0=clock;
z = L\(P*u);
y = U\z;

x2m = x1 + y*((v'*x1)/(1-v'*y));
t2(i) = etime(clock,t0);

norm(x2-x2m)

end

hold on;
xlabel('N');
ylabel('Time required');
plot(sizes,t1(sizes),'ro-',sizes,t2(sizes),'b.-');
legend('direct inversion','Sherman-Morrison update');


