n1 = 15;
A1 = hilb( n1 );
k2_1 = cond( A1 )
%   k2_1 =
%      4.4333e+17

x1_ex = ones( n1, 1 );
b1 = A1 * x1_ex;
x1_0 = zeros( n1, 1 );
% criterion based on RELATIVE residual
tol_rel = 1.0e-5;       kmax = 1000;
tol =  tol_rel * norm( b1 );
[ x1_gs, k1_gs, res1_gs ] = gauss_seidel( A1, b1, x1_0, tol, kmax );
k1_gs
%   k1_gs =
%      599
res_rel1_gs = res1_gs / norm( b1 )
%   res_rel1_gs =
%      9.9853e-06
err_rel1_gs = norm( x1_ex - x1_gs ) / norm( x1_ex )
%   err_rel1_gs =
%      0.0412

n2 = 100;
A2 = diag( ( 4 + 5e-5 ) * ones( n2, 1 ), 0 ) + ...
     diag( - 2 * ones( n2 - 1, 1 ), 1 ) + diag( - 2 * ones( n2 - 1, 1 ), -1 );
D2 = diag( diag( A2 ) );
E2 = - tril( A2, -1 );
B2_GS = eye( n2 ) - inv( D2 - E2 ) * A2;
format short e
rho2_GS = max( abs( eig( B2_GS ) ) )
%   rho2_GS =
%      9.9901e-01
format

x2_ex = ones( n2, 1 );
b2 = A2 * x2_ex;
tol = 1.0e-5;    kmax = 10000;
x2_0 = zeros( n2, 1 );
[ x2_gs, k2_gs, diff2_gs ] = gauss_seidel_difference_iterates( A2, b2, x2_0, ...
                                                               tol, kmax );
k2_gs, diff2_gs,
%   k2_gs =
%           6852
%   diff2_gs =
%      9.9964e-06
format short e
err2_gs = norm( x2_ex - x2_gs )
%   err2_gs =
%      1.0065e-02
format