> restart; 
> with(LinearAlgebra): 
> Relax := proc (A, b, tol, r) 
> local n, k, vectornorm, residuum, RelaxIt; 
> global res, vnorm, x; 
> RelaxIt := proc (A, x, b, r) 
> local j, n, k, z, y; 
> n := Dimension(b); 
> y := x; 
> for j to n do z := b[j]; 
>     for k 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]*y[k] end do; 
>     y[j] := r*z/A[j, j]+(1-r)*y[j] 
> end do; 
> for j to n do x[j] := y[j] end do;
> end proc; 
> vectornorm := proc (x) 
> local vnorm, j, n; 
> vnorm := 0; 
> n := Dimension(x); 
> for j to n do if vnorm < abs(x[j]) then vnorm := abs(x[j]) end if 
> end do; 
> eval(vnorm) 
> end proc; 
> residuum := proc (A, x, b) 
> local j, k, n; global res; 
> n := Dimension(b); 
> res := Vector(n); 
> for j to n do 
>     for k 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); 
> residuum(A, x, b); 
> vnorm := vectornorm(res); 
> k := 0; 
> while tol < vnorm do 
>      RelaxIt(A, x, b, r); 
>      residuum(A, x, b); 
>      vnorm := 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;  
> Relax(A,b,tol, 0.9);
> 
>