% GAUSS-SEIDEL EXAMPLES
%  (code by Benjamin Walter, METU-NCC, benjamin@metu.edu.tr)
%
% The Gauss-Seidel method is an iterative algorithm used to solve systems
% of equations.  It works as long as the matrix A has big elements on the
% diagonal and small elements elsewhere (it also works for symmetric, 
% positive definite matrices, but we haven't discussed them yet).  Matrices
% like this commonly arise when solving systems of differential equations,
% for example: linked spring systems and heat diffusion systems.
%
%
%% EXAMPLE 1: Simple iteration, by hand
%
clc

fprintf('Solution via simple iteration, by hand')

A=[3 .2 -.1; .3 4 .2; .1 -.4 2]
b=[1 1 -1]'    % Recall: x' is x transpose. Fastest way to make col vector.

%%
% Initial guess:
%
fprintf('Initial guess:')

x = [1 1 -1]'  % Initial guess doesn't matter much, as we'll see later.

%%
% Check our guess (if we guessed the right answer, we don't need to bother
% with any GS iteration) :)
%
fprintf('The error of our guess is b-Ax:')
Error = b-A*x

SizeOfError = norm(Error)

fprintf(' -->That wasn''t a very good guess... now we''ve got some work to do.\n')


%%
% Now we will refine x using info from Ax=b. 
%
% Let's print out the theory in the command window, so you can see what's
% going on clearly.
%
fprintf('-----\nAugmented matrix A|b:')

Ab = [A b]

fprintf('The equation Ax=b tells us that:\n')
fprintf(' 3*x(1) + .2*x(2) - .1*x(3) =  1 --> x(1) = ( 1 - .2*x(2) + .1*x(3))/3\n')
fprintf('.3*x(1) +  4*x(2) + .2*x(3) =  1 --> x(2) = ( 1 - .3*x(1) - .2*x(3))/4\n')
fprintf('.1*x(1) - .4*x(2) +  2*x(3) = -1 --> x(3) = (-1 - .1*x(1) + .4*x(2))/2\n')
%
%%
% Refinement 1:
%
fprintf('-------------\nREFINEMENT 1:')

x(1) = ( 1 - .2*x(2) + .1*x(3))/3;
x(2) = ( 1 - .3*x(1) - .2*x(3))/4;
x(3) = (-1 - .1*x(1) + .4*x(2))/2;

x
Error = b-A*x
SizeOfError = norm(Error)

fprintf('  We''re still not very close to the right answer...\n-------------\n')

%%
% We need to iterate this a lot more times to get closer to the solution!
%
%  Just for kicks, let's do this by defining INLINE FUNCTIONS.  Inline
%  functions are simple, short functions which are defined on the fly in
%  matlab rather than inside of .m files.  The format for them is
% <function name> = inline('<expression to be evaluated>', '<variable name>') 
%
fprintf('Inline functions for refinement:')

RefineX1 = inline('( 1 - .2*X2 + .1*X3)/3', 'X2', 'X3')  % This is not 
RefineX2 = inline('( 1 - .3*X1 - .2*X3)/4', 'X1', 'X3')  % necessary, but
RefineX3 = inline('(-1 - .1*X1 + .4*X2)/2', 'X1', 'X2')  % it is cute.

%%
% Now let's iterate a few times
%
for (i = 1:4)
    fprintf('-------------\nREFINEMENT %i:', i+1)
    
    x(1) = RefineX1(x(2),x(3));  % Note:  we could also have just typed in
    x(2) = RefineX2(x(1),x(3));  %  (or cut and pasted) the formulas from
    x(3) = RefineX3(x(1),x(2));  %  above.  But, inline functions are fun!
    x
    Error = b-A*x
    SizeOfError = norm(Error)
end

fprintf('Pretty good.  (We probably could have stopped at refinement 3.)\n')



%% EXAMPLE 2: Simple iteration, automatic (for the people)
% It's not hard to make this method fully automated.
%
clc

fprintf('Solution via simple iteration, automated')

A = [3  .1  .3 -.4  .3  
    .1   2 -.1  .2 -.1
    .2 -.1   4  .3 -.2
    .1  .3  .1   5  .2
    .1  .2 -.2 -.3  10]

b = [3 1 4 1 5]'   % pi

x = [2 7 1 8 2]'   % e

%%
%  Time to refine our guess:
%    We'll refine our answer until the error is really, really small.
%
Error = b-A*x;
SizeOfError = norm(Error);
[numrows,numcols] = size(A);

NumberOfIterations = 0; % just for kicks, we'll count our iterations.

while (SizeOfError > .000001)
    for (i=1:numcols)
        x(i) = ( b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]) )/A(i,i);
    end
    NumberOfIterations = NumberOfIterations + 1;
    
    Error = b-A*x;
    fprintf(' ---> Refinement %i:', NumberOfIterations);
    x
    SizeOfError = norm(Error)
end

%%
% Slight improvement on this code:
%  In the code above, every time we refine x(i), we divide by A(i,i).
%  we can speed things up slightly dividing row i of A and b by A(i,i) 
%  once at the beginning, so that we don't have to repeat this division 
%  over and over.
%
fprintf('------\nWe can speed up our code a little by normalizing Ax=b before solving.\n')
b = diag(A).\b       % divide b by diagonal elements of A
A = diag(diag(A))\A  % divide each row by it's diagonal element


Error = b-A*x;
SizeOfError = norm(Error);
[numrows,numcols] = size(A);

NumberOfIterations = 0; % just for kicks, we'll count our iterations.

while (SizeOfError > .000001)
    for (i=1:numcols)
        x(i) = b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]); % no division here
    end
    NumberOfIterations = NumberOfIterations + 1;
    
    Error = b-A*x;
    fprintf(' ---> Refinement %i:', NumberOfIterations);
    x
    SizeOfError = norm(Error)
end

%% COMMENT on vectorization and good algorithms
% Note that the above code uses the refined x(1) when refining x(2), so 
% without doing something clever, we cannot vectorize it to refine x(1) 
% and x(2) at the same time...  If we allow ourselves to be sloppy, though, 
% we can vectorize by using the unrefined x(1) when computing x(2):
% The formula for x(i) is: 
%  x(i) = b(i) - (row i of A without element i,i) * (x without element i)
% If we make a matrix that looks like A but with 0 on the diagonal then we 
% can refine x via the vector command:
%  x = b - (A with zero diagonal)*x 
% This won't converge as fast -- since it refines x(2) and x(1) at the same
% time, it cannot use the refined x(1) in the calculation of x(2).....
% Still, maybe the speed gained by being vectorized will make up for the 
% extra iterations required.
%
% Let's check!
%
clc
fprintf('A is a random 500 x 500 matrix with small off-diagonal entries\n')
fprintf('b is a random 500 vector\n')


% We may as well start out with A normalized so that it's diagonal elements
% are all 1 -- otherwise our first step would be to normalize everything
% anyway.  
% -->The command "tril(A,-1)" returns a matrix which equals A everywhere 
%   below the diagonal, and is zero on the diagonal and above.
%    (it is "tril" because it returns the Lower TRIangle of A)
% -->The command "triu(A,1)" returns a matrix which equals A everywhere
%   above the diagonal, and is zero on the diagonal and below.
%    (it is "triu" because it returns the Upper TRIangle of A)

A = randn(500)/10; A = tril(A,-1) + triu(A,1) + eye(500); 
b = rand(500,1)*10;

%%
%  first we'll check the gauss-seidel method using the refinement
%  x(i) = b(i) - (A with zero diagonal)*x 

x = zeros(500,1);

fprintf('Time required for Gauss-Seidel to be closer than .0000001:\n')
start = now;

SizeOfError = norm(b-A*x);
[numrows,numcols] = size(A);
NumberOfIterations = 0; 
while (SizeOfError > .0000001)
    for (i=1:numcols)
        x(i) = b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]);
    end
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end

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



%%
%  now we'll check a brute-force, vectorized method
%
x = zeros(500,1);

fprintf('-----\nTime required for vectorized iteration to be closer than .0000001:\n')
start = now;
AMinusDiag = tril(A, -1) + triu(A, 1);

SizeOfError = norm(b-A*x);
[numrows,numcols] = size(A);
NumberOfIterations = 0; 
while (SizeOfError > .000001)
    x = b - AMinusDiag*x;
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end
time=[minute(now-start), second(now-start)]
NumberOfIterations


%%
% If we are clever, there is another way to vectorize Gauss-Seidel.
% Looking carefully at the Gauss-Seidel code, you can see it solves:
%  tril(A)*NewX = b - triu(A,1)*OldX
% --> the for loop is just computing NewX using back substitution, while
%     reusing the vector x to store the solution.
% So we could try to vectorize by doing:
%  x = inv(tril(A))*( b - triu(A,1)*x )
% The issue with this is that you must first compute the inverse of
% tril(A,-1)... which is a bit of work by itself.

x = zeros(500,1);

fprintf('-----\nTime required for vectorized GS to be closer than .0000001:\n')
start = now;
TrilAInv = inv(tril(A));  TriuA = triu(A,1);

SizeOfError = norm(b-A*x);
[numrows,numcols] = size(A);
NumberOfIterations = 0; 
while (SizeOfError > .000001)
    x = TrilAInv*( b - TriuA*x);
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end
time=[minute(now-start), second(now-start)]
NumberOfIterations

fprintf('This comparison isn''t really fair, because we used matlab to compute inv.\n')
fprintf('If we used a script to compute inverse, it would probably take much longer.\n')

%%
%  let's check this against matlab's black magic
%
fprintf('-----\nTime required for matlab''s (hard-coded) algorithm:\n')
start=now;
A\b;
time=[minute(now-start), second(now-start)]

fprintf('This comparison isn''t really fair, since matlab''s algorithm is \n')
fprintf('written in C or fortran and already compiled to assembly code, \n')
fprintf('ready to be fed into a computer''s processor.  So of course it \n')
fprintf('will be faster than a script which must be compiled and then run. \n')




%% COMMENT on choosing initial guess:
% If the matrix A is dominated by it's diagonal elements, a good initial
% guess is  x(i) = b(i)/A(i,i).  If A consisted of only the diagonal 
% elements, this x would be the answer.
%
clc
fprintf('A is a random 500 x 500 matrix with small off-diagonal entries\n')
fprintf('b is a random 500 vector\n')

A = randn(500)/2; A = A + diag(diag(A)*30);
b = rand(500,1)*90;

%%
%
fprintf('Initial guess of all zeros')

x = zeros(500,1);

Error = b-A*x;
SizeOfError = norm(Error);
[numrows,numcols] = size(A);
NumberOfIterations = 0;
while (SizeOfError > .000001)
    for (i=1:numcols)
        x(i) = ( b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]) )/A(i,i);
    end
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end

NumberOfIterations

%%
%
fprintf('--------\nx = b/diag(A)')

x = b./diag(A);  % this would be the solution if A was purely diagonal.

Error = b-A*x;
SizeOfError = norm(Error);
[numrows,numcols] = size(A);
NumberOfIterations = 0;
while (SizeOfError > .000001)
    for (i=1:numcols)
        x(i) = ( b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]) )/A(i,i);
    end
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end

NumberOfIterations
%%
%
fprintf('--------\nA initial guess of all 1000:')

x = ones(500,1)*1000;

Error = b-A*x;
SizeOfError = norm(Error);
[numrows,numcols] = size(A);
NumberOfIterations = 0;
while (SizeOfError > .000001)
    for (i=1:numcols)
        x(i) = ( b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]) )/A(i,i);
    end
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end

NumberOfIterations

%%
%
fprintf('--------\nA initial guess of random 3 digit numbers:')

x = rand(500,1)*1000;

Error = b-A*x;
SizeOfError = norm(Error);
[numrows,numcols] = size(A);
NumberOfIterations = 0;
while (SizeOfError > .000001)
    for (i=1:numcols)
        x(i) = ( b(i) - A(i,[1:i-1 i+1:end])*x([1:i-1 i+1:end]) )/A(i,i);
    end
    NumberOfIterations = NumberOfIterations + 1;
    SizeOfError = norm(b-A*x);
end

NumberOfIterations

fprintf('------->Moral of the story:  Initial guess doesn''t matter much!\n')
