%%% Implementation of Gaussian elimination with Partial Pivoting, Jacobi,
%%% and Gauss-Seidel iteration methods for solving Ax = b; 
clear; clc
n = 10; 
%n = 1e6; % (2nd option)
e = ones(n,1);

if n < 1e6
    A = spdiags([-e 3*e -e], -1:1, n, n) + 1/2*fliplr(eye(n)); % only works for small n
else 
    
    
end 
x = ones(n,1); % true solution
b = A*x;


%% Gaussian Elimination with partial pivoting

% Augmented System: S 
S = [A b]; 
n = size(A,1); 
xPP = zeros(n,1); 

% Partially Pivot S down into an Upper Triangular Matrix
for j = 1:n-1 % only pivot first part of Augmented system
    [M,k] = max(S(:,j)); % biggest col value
    
    temp1 = S(k,:); 
    temp2 = S(j,:); 
    S(j,:) = temp1; % swap rows 
    S(k,:) = temp2;
  
    S(j,:) = S(j,:)./S(j,j); %normalize pivot
    % Eliminate Downwards
    for i = j+1:n 
        
    S(i,:) = -S(j+1,j).*S(i-1,:) + S(i,:);
    end 
end 

% Backsubstitute S
y = S(:,end); % Solve U*xPP = y 
U = S(1:end,1:end-1); 
for j=n:-1:1
    xPP(j)=y(j)/U(j,j);
    y(1:j-1)=y(1:j-1)-U(1:j-1,j)*xPP(j);
end

err = norm(A*xPP - b,Inf)/norm(b); 
String = ['GaussPP:',' error = ', num2str(err)];
disp(String)
%% Jacobi iteration 
tol = 1e-6; 
xJ = zeros(n,1); 
max_iter = 100; 
i = 1; 
err = Inf; 

S = diag(diag(A)); % use just the diagonal part of A 
T = S - A; 
tic
while i < max_iter && err > tol
    xJ=S \ (b+T*xJ);
    err = norm(A*xJ-b,Inf)/norm(b,Inf); 
    i = i + 1; 
    if i == max_iter
        fprintf('No Convergence')
    end 
end
t = toc; 


String = ['Jacobi: i = ',num2str(i),', error = ', num2str(err), ... 
    ', time = ', num2str(t), ' s'];
disp(String)

%% Gauss-Seidel iteration 
tol = 1e-6; 
xGS = zeros(n,1); 
max_iter = 100; 
i = 1; 
err = Inf; 


S = tril(A); % use lower triangular part of A 
T = S - A; 
tic
while i < max_iter && err > tol
    xGS=S \ (b+T*xGS);
    err = norm(A*xGS-b,Inf)/norm(b,Inf); 
    i = i + 1; 
    if i == max_iter
        fprintf('No Convergence')
    end 
end
t = toc; 

String = ['Gauss-Seidel: i = ',num2str(i),', error = ', num2str(err), ... 
    ', time = ', num2str(t), ' s'];
disp(String)