Classical Iterative Solvers

The idea behind all these schemes is an iterative scheme based on an additive decomposition of a matrix.  The easiest decomposition to understand is
	A=L+U+Dia
where L is strictly lower triangular,  U is strictly uper triangular, and Dia is diagonal.  If 
	A.x=b
then 
	Dia.x=b-(L+U).x.
Dia is easily invertible (as long as it has non-zero entries) and
	x=(Dia^(-1)).(b-(L+U).x)
Turning this into an interative scheme
	(Subscript[x, n+1])=(Dia^(-1)).(b-(L+U).(Subscript[x, n]))
gives Jacobi’s method.  Slightly more complicated is the Gauss-Seidel iteration.
	(Subscript[x, n+1])=(((L+Dia))^(-1)).(b-U.(Subscript[x, n])).
The final iteration in the trio of classical methods is Successive Over Relaxation (SOR) which starts from the decomposition
	ω A=(ω L+Dia)+(ω U+(ω-1)Dia)
so 
	ω A.x=ω b
becomes
	(Subscript[x, n+1])=(((ω L+Dia))^(-1)).(ω b-(ω U + (ω-1)Dia).(Subscript[x, n])).
Note: this is  Gauss-Seidel for ω=1.

All of these converge for matrics with large diagonal entries.  A matrix is Diagonally Dominant if the diagonal entry is larger in magnitude than all the other entries on each row.

Diagonaly Dominant Test Matrices

Laplacian1D[n_]:= SparseArray[{Band[{1,1}]:>-2.0, Band[{1,2}]->1.0, Band[{2,1}]->1.0},{n,n}]

SparseDDominant[n_,nzr_]:= Module[
{A},
(*nzr is the number of non-zeroes per row.*)
A=SparseArray[Band[{1,1}]:> RandomChoice[{-1,1}]RandomReal[{1.1nzr,5 nzr}],{n,n}];
Do[
Do[
While[
j=RandomInteger[{1,n}];
A[[i,j]]≠0.0];
A[[i,j]]=RandomReal[{-1,1}],
{j,nzr}
],
{i,1,n}];
A
]

MatrixPlot[SparseDDominant[6,3]]
MatrixPlot[Laplacian1D[12]]

Jacobi, Gauss-Seidel, SOR in Matrix form

If  and the matrix A is already splot up we have the following iterations.

JaM[{L_,U_,Dia_},x_,b_]:= LinearSolve[Dia,(b-(L+U).x)]
GSM[{L_,U_,Dia_},x_,b_]:= LinearSolve[ L+Dia,(b-U.x)]
SORM[{L_,U_,Dia_},x_,b_,ω_]:= LinearSolve[L+ (Dia/ω) ,b-(U-((1-ω)/ω)Dia ).x]

Testing that these are

n=123;
A=Laplacian1D[n];
{L,U,Dia}={LowerTriangularize[A,-1],UpperTriangularize[A,1], DiagonalMatrix[Diagonal[A]]};
Norm[A-L-U-Dia];
xSol=RandomReal[{-1,1},n];
b=A.xSol;
ω=1.2;
Norm[JaM[{L,U,Dia},xSol,b]-xSol]
Norm[GSM[{L,U,Dia},xSol,b]-xSol]
Norm[SORM[{L,U,Dia},xSol,b,ω]-xSol]

x0=RandomReal[{-1,1},n];

MaxIter=200;
x=x0;
JData=Table[
x=JaM[{L,U,Dia},x,b];
Norm[x-xSol],
{MaxIter}];
x=x0;
GSData=Table[
x=GSM[{L,U,Dia},x,b];
Norm[x-xSol],
{MaxIter}];
ListLogLogPlot[ {JData,GSData},
PlotRange->All,
PlotStyle->{Green, Blue},
Joined->True,
PlotMarkers->Automatic]

SOR is a little more complicated because of the parameter ω

{ωMin,ωMax}={0.8,1.9};
MaxIter=30;
SORData=Table[
x=x0;
Table[
x=SORM[{L,U,Dia},x,b,ω];
Norm[x-xSol],
{MaxIter}],
{ω,ωMin,ωMax,(ωMax-ωMin)/24}];
ListContourPlot[Log[10,SORData],
FrameLabel->{"n","ω"},
ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),
DataRange->{{1,MaxIter},{ωMax,ωMin}}]

{ωMin,ωMax}={1.9,1.99};
MaxIter=800;
SORData=Table[
x=x0;
Table[
x=SORM[{L,U,Dia},x,b,ω];
Norm[x-xSol],
{MaxIter}],
{ω,ωMin,ωMax,(ωMax-ωMin)/12}];
ListContourPlot[Log[10,SORData],
FrameLabel->{"n","ω"},
(*ContourLabels->(Text[Framed[#3],{#1,#2},Background->White]&),*)
DataRange->{{1,MaxIter},{ωMax,ωMin}}]

So the value of ω matters!

MaxIter=5000;
x=x0;
ω=1.948;
SORDataGood=Table[
x=SORM[{L,U,Dia},x,b,ω];
Norm[x-xSol],
{MaxIter}];
ω=1.8;
x=x0;
SORDataBad=Table[
x=SORM[{L,U,Dia},x,b,ω];
Norm[x-xSol],
{MaxIter}];
x=x0;
GSData=Table[
x=GSM[{L,U,Dia},x,b];
Norm[x-xSol],
{MaxIter}];
ListLogLogPlot[{GSData, SORDataGood,SORDataBad},
PlotRange->All]