Off[General::spell1];
Clear[SOR,a0,b0, r0, M0];
SOR[a0_,b0_, r0_, M0_] :=
Module[{a,b,k,n,x,u,r},
  a = a0; b = b0; M = M0; r = r0;
  n = Length[b];
  x = Table[0,{i, n}];
  u = Table[0, {i,n}];
  Clear[i,k];
  For[k=1,k<=M,k++,
    For[i=1, i<=n, i++,
       u[[i]] = (b[[i]] - Sum[a[[i,j]] x[[j]], {j,1,i-1}]
         - Sum[a[[i,j]] x[[j]], {j,i+1,n}])/a[[i,i]];
       x[[i]] = x[[i]] + r (u[[i]] - x[[i]]) ];
    Print[N[x,10]];  ]; 
  Return[x]; ]  
  On[General::spell1];

Clear[A,B];  
A = {{8,-1,1},{-1,8,1},{1,1,6}};
B = {8,8,8};  
StringForm["A= ``    B = ``", MatrixForm[A], MatrixForm[B]]

N[SOR[A,B,1,10],10]

Clear[n,Q,rho]; 
n=Length[A]; 
Q =  DiagonalMatrix[Table[1/A[[i,i]], {i,1,n}]];
rho=Max[N[Abs[Eigenvalues[IdentityMatrix[3] - Q. A]],10]];
omega= 2/(1+ (1- rho^2)^(1/2) )

SOR[A,B,omega,10]

Clear[A,B];  
A = {{1,0,-1/4,-1/4},{0,1,-1/4,-1/4},{-1/4,-1/4,1,0},
      {-1/4,-1/4,0,1}};
B = {1,1,1,1};   
StringForm["A= ``    B = ``", MatrixForm[A], MatrixForm[B]]

N[SOR[A,B,1,10],10]

Clear[n,Q,rho]; 
n=Length[A];
Q =  DiagonalMatrix[Table[1/A[[i,i]], {i,1,n}]];
rho=Max[N[Abs[Eigenvalues[IdentityMatrix[n] - Q. A]],10]];
omega= 2/(1+ (1- rho^2)^(1/2) )

N[SOR[A,B,omega,10],10]

Clear[Q]; 
L=Table[If[i >j, -A [[i,j]], 0], {i,n},{j,n}];   
Q =(DiagonalMatrix[Table[A[[i,i]],{i,1,n}]]-omega L)/omega;
Max[N[Abs[Eigenvalues[IdentityMatrix[n]-Inverse[Q]. A]],10]]

Clear[A,B];
A={ {1,1/2,1/3}, {1/3,1,1/2},{1/2,1/3,1} };   
B={11/18,11/18,11/18};

SOR[A,B,1,10]//N

Clear[n,Q,rho]; 
n=Length[A];
Q =  DiagonalMatrix[Table[1/A[[i,i]], {i,1,n}]];
rho=Max[N[Abs[Eigenvalues[IdentityMatrix[3] - Q. A]],10]];
omega= 2/(1+ (1- rho^2)^(1/2) )

N[SOR[A,B,omega,10],10]

Clear[Q]; 
L=Table[If[i >j, -A [[i,j]], 0], {i,n},{j,n}];   
Q =(DiagonalMatrix[Table[A[[i,i]],{i,1,n}]]-omega L)/omega;
Max[N[Abs[Eigenvalues[IdentityMatrix[n]-Inverse[Q]. A]],10]]