n = 4;
A = diag( 5 * ones( n, 1 ), 0 ) + diag( 1 * ones( n - 1, 1 ), 1 ) + ...
    diag( 1 * ones( n - 1, 1 ), -1 ) + diag( 1 * ones( n - 2, 1 ), 2 ) + ...
    diag( 1 * ones( n - 2, 1 ), -2 );
x_ex = ones( n, 1 );
b = A * x_ex;
tol = 1.0e-6;   kmax = 100;
x0 = zeros( n, 1 );
% gradient method (P=I)
P1 = eye( n );
[ x1, k1, res1 ] = preconditioned_gradient( A, b, P1, x0, tol, kmax );
err1 = norm( x_ex - x1 )
%   err1 =
%      1.6781e-08
k1, res1
%   k1 =
%        6
%   res1 =
%      1.2614e-07

% preconditioned gradient method (P=P_2)
P2 = diag( diag( A ) ) + diag( diag( A, 1 ), 1 ) + diag( diag( A, - 1 ), - 1 );
[ x2, k2, res2 ] = preconditioned_gradient( A, b, P2, x0, tol, kmax );
err2 = norm( x_ex - x2 )
%   err2 =
%      2.1285e-08
k2, res2
%   k2 =
%      4
%   res2 =
%      1.5999e-07

res_PG1_v = [];    res_PG2_v = [];
res_J_v = [];      res_GS_v = [];
klimit = 50;  tol = 1e-14;
k_vect = 1 : klimit;
for kmax = k_vect
    % gradient
    [ xPG1, kPG1, resPG1 ] = preconditioned_gradient( A, b, P1, x0, tol, kmax );
    res_PG1_v = [ res_PG1_v, resPG1 ];
    % preconditioned gardient
    [ xPG2, kPG2, resPG2 ] = preconditioned_gradient( A, b, P2, x0, tol, kmax );    
    res_PG2_v = [ res_PG2_v, resPG2 ];
    % Jacobi
    [ xJ, kJ, resJ ] = jacobi( A, b, x0, tol, kmax );
    res_J_v = [ res_J_v, resJ ];
    % Gauss-Seidel
    [ xGS, kGS, resGS ] = gauss_seidel( A, b, x0, tol, kmax );    
    res_GS_v = [ res_GS_v, resGS ];
end
semilogy( k_vect, res_PG1_v, '-b', k_vect, res_PG2_v, '-r', ...
          k_vect, res_J_v, '-m', k_vect, res_GS_v, '-k'  );
axis( [ 1 klimit 1e-13 10 ])
legend('Gradient', 'Prec.Gradient', 'Jacobi', 'Gauss-Seidel' );