> restart: 
> with(LinearAlgebra): 
> printlevel:=0:
> #Gauss-Seidel-Verfahren 
> GaussSeidel:=proc(A,b,tol) 
> local n,k, vectornorm, residuum, GaussSeidelIt;
> global res, vnorm ,x;
> # option trace;
> GaussSeidelIt:=proc(A,x,b) #Gauss-Seidel-Iterationsschritt
> local j,n,k,z,y;
> # option trace;
> n:=Dimension(b); 
> y:=Vector(n);
> for j from 1 to n do
>   z:=b[j];
>   for k from 1 to j-1 do z:=z-A[j,k]*y[k] end do;
>   for k from j+1 to n do z :=z-A[j,k]*x[k] end do;
>   y[j]:=z/A[j,j]
> end do;
> for j from 1 to n do x[j]:= y[j] end do;
> end proc:
> vectornorm:=proc(x)  #Maximumnorm
> global  vnorm;
> local j, n;
> # option trace;
> vnorm:=0; 
> n:=Dimension(x); 
> for j from 1 to n do
>   if abs(x[j]) >vnorm then vnorm :=abs(x[j]) end if
> end do;
> end proc;
> residuum:=proc(A,x,b)   #Residuum Ax-b
> global res;
> local j,k,n;
> # option trace;
> n:=Dimension(b); 
> res:=Vector(n);
> for j from 1 to n do
>   for k from 1 to n do res[j]:= res[j]+A[j,k]*x[k] end do;
>   res[j]:=res[j]-b[j];
> end do;
> end proc;
> n:=Dimension(b); 
> x:=Vector(n); 
> res:=Vector(n); 
> residuum(A,x,b);
> vectornorm(res); 
> k:=0;
> while (vnorm > tol) do
>    GaussSeidelIt(A,x,b);
>    residuum(A,x,b);
>    vectornorm(res);    
>    k:=k+1;
>   print(k, x); 
> end do; 
> end proc:
> A:=Matrix([[10.,-1.,0.,2.],[1.,12.,-1.,2.],[-2.,1.,15.,0.],[1.,-2.,0.,20.]]): b:=Vector([11,14,14,19]):  tol:=0.001:  GaussSeidel(A,b,tol);
> 
>