% Gauss-Seidelova & Jacobiho iteracni metoda
A = [20 1 2 4;
     5 30 1 8;
     2 8 25 1;
     1 2 8 42];

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

x_exact = A\b;
x = [0;0;0;0];
N = length(x);

% Jacobiova & Gauss-Seidelova metoda
iter = 100;

x_J  = x;
x_GS = x;

for k=1:iter
    xk = x_J;
    for i=1:N
        sum_J  = 0;
        sum_GS = 0;
        for j=1:N
            if(j ~= i)
                sum_J =  sum_J  + A(i,j)*xk(j);
                sum_GS = sum_GS + A(i,j)*x_GS(j);
            end
        end
        x_J(i)  = 1/A(i,i)*(b(i) - sum_J);
        x_GS(i) = 1/A(i,i)*(b(i) - sum_GS);
    end
    err_J(k) = norm(x_exact-x_J,2);
    err_GS(k) = norm(x_exact-x_GS,2);
end

x_exact
x_J
x_GS

hold on;
plot(1:iter, err_J,'-r')
plot(1:iter, err_GS,'-b')
grid on
yscale log
hold off;

% Pozn.: PIM
U = eye(4) - A;
x = [0;0;0;0];
for k=1:iter
    x = U*x + b
end

% Pozn.: PIM nekonverguje!! Proc?
% zkusme spektralni polomer:
max(abs(eig(U)))
