Clear[a,x,xx]
n=8;rt2=1./Sqrt[2.0];rt3=1./Sqrt[3.0];
a=Table[0.0,{i,1,n},{j,1,n}];
a[[1,1]]=-1.; a[[1,4]]=rt2;a[[1,5]]=1.;
a[[2,2]]=-1.;a[[2,4]]=rt2;
a[[3,3]]=-1.;a[[3,7]]=0.5;
a[[4,4]]=-rt2;a[[4,6]]=-1.;a[[4,7]]=0.5;
a[[5,5]]=-1.;a[[5,8]]=1.;
a[[6,6]]=1.;
a[[7,4]]=-rt2;a[[7,7]]=rt3;
a[[8,7]]=-rt3;a[[8,8]]=-1.;
b=Table[0.0,{i,1,n}];
b[[6]]=10000.0
MatrixForm[a]
MatrixForm[b]

x=Table[xx[i],{i,1,n}];
Det[a];
sys=Thread[Dot[a,x]==b];
xsol=x/.Solve[sys,x][[1]]//N

d=Diagonal[a]*IdentityMatrix[n];
r=a-d;
l=u=Table[0,{i,1,n},{j,1,n}];
Do[l[[i,j]]=-a[[i,j]],{i,2,n},{j,1,i-1}];
Do[u[[i,j]]=-a[[i,j]],{i,1,n-1},{j,i+1,n}];
MatrixForm[d];
MatrixForm[l];
MatrixForm[u];
tJ=-Dot[Inverse[d],r];
cJ=Dot[Inverse[d],b];
λJ=Eigenvalues[tJ]//N;
ρJ=Max[Abs[λJ]]
tGS=Dot[Inverse[d-l],u];
cGS=Dot[Inverse[d-l],b];
λGS=Eigenvalues[tGS]//N;
ρGS=Max[Abs[λGS]]
ω0=2.0/(1.0+Sqrt[1.0-ρJ^2])

ω=1.25;
tSOR=Dot[Inverse[d-ω*l],ω*u+(1-ω)*d];
cSOR=Dot[ω*Inverse[d-ω*l],b];
λSOR=Eigenvalues[tSOR]//N;
ρSOR=Max[Abs[λSOR]]

t=tJ;
c=cJ;
kmax=20000;tol=1.0 10^(-16);
k=0;xk[0]=Table[0.0,{i,1,n}];
While[k<kmax,
xk[k+1]=Dot[t,xk[k]]+c;
norm=Max[Abs[xk[k+1]-xk[k]]];If[norm≤tol,Break[]];
k++]
Print[k,"  ",xk[k]]
xsol

t=tGS;
c=cGS;
kmax=20000;
k=0;xk[0]=Table[0.0,{i,1,n}];
While[k<kmax,
xk[k+1]=Dot[t,xk[k]]+c;
norm=Max[Abs[xk[k+1]-xk[k]]];If[norm≤tol,Break[]];
k++]
Print[k,"  ",xk[k]]
xsol

t=tSOR;
c=cSOR;
kmax=10000;
k=0;xk[0]=Table[0.0,{i,1,n}];
While[k<kmax,
xk[k+1]=Dot[t,xk[k]]+c;
norm=Max[Abs[xk[k+1]-xk[k]]];If[norm≤tol,Break[]];
k++]
Print[k,"  ",xk[k]]
xsol

(* Stesso con Gauss-Jordan *)

Clear[ab,x];
ab=Table[Join[a[[i]],{b[[i]]}],{i,1,n}];
MatrixForm[ab]

Do[piv=ab[[i]][[i]];Print["i = ", i, "  piv = ", piv];Do[m[j]=ab[[j]][[i]]/piv;ab[[j]]=ab[[j]]-m[j]*ab[[i]],{j,i+1,n}],{i,1,n-1}];
MatrixForm[ab]

xx[n]=ab[[n,n+1]]/ab[[n,n]];
Do[xx[i]=(ab[[i,n+1]]-Sum[xx[j]*ab[[i,j]],{j,i+1,n}])/ab[[i,i]],{i,n,1,-1}];
x=Table[xx[i],{i,1,n}];
MatrixForm[x]