function x=matrixGS(A,b,n)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
BG=(D-L)\U;
f=(D-L)\b;
x=zeros(n,1);
rx=BG*x+f;
while norm(rx-x)>=1.0e-6
    x=rx;
    rx=BG*x+f;
end
x=rx;