% Section 3.1

A=[5 -2 3; -3 9 1; 2 -1 -7]
D=diag(diag(A));
U=triu(A)-D;
L=tril(A)-D;

b=[-1;2;3];

% Jacobi iteration
xj0=[0;0;0];
xj1=D\(-(U+L)*xj0+b);
rj=norm(xj1-xj0)/norm(xj1);
count=0;
while rj>1e-6 && count<100
    count=count+1;
    xj0=xj1;
    xj1=D\(-(U+L)*xj0+b);
    rj=norm(xj1-xj0)/norm(xj1);
end
count
rj
max(abs(eig(D\(-(U+L))))) % spectral radius

% Gauss-Seidel iteration
xgs0=[0;0;0];
xgs1=(L+D)\(-U*xgs0+b);
rgs=norm(xgs1-xgs0)/norm(xgs1);
count=0;
while rgs>1e-6 && count<100
    count=count+1;
    xgs0=xgs1;
    xgs1=(L+D)\(-U*xgs0+b);
    rgs=norm(xgs1-xgs0)/norm(xgs1);
end
count
rgs
max(abs(eig((L+D)\(-U)))) % spectral radius

% SOR
w=1.1;
max(abs(eig((D+w*L)\((1-w)*D-w*U)))) % spectral radius
xgs0=[0;0;0];
xgs1=(w*L+D)\(((1-w)*D-w*U)*xgs0)+w*((D+w*L)\b);
rgs=norm(xgs1-xgs0)/norm(xgs1);
count=0;
while rgs>1e-6 && count<100
    count=count+1;
    xgs0=xgs1;
    xgs1=(w*L+D)\(((1-w)*D-w*U)*xgs0)+w*((D+w*L)\b);% or (w*L+D)\(w*b-(w*U+(w-1)*D)*xgs0);
    rgs=norm(xgs1-xgs0)/norm(xgs1);
end
count
rgs