close all
clear all
m_max = 100; eps = 1e-8;

legend_cell = {};  % for the legend of the plot
k = 1;


for N = [10,20,40]
    % get A
    B_cell = diff_advec(N);
    n = N^2;
    tau = (1/N+1)^2;
    A = speye(n) + tau*B_cell{1}{4};
    H = 0.5*(A+A');
    gamma_low = tau;
    gamma_high = max(sum(abs(H)));  %Gershgorin
    %evals = eig(full(H));   % works only for small matrices
    %gamma_low = min(evals); gamma_high = max(evals);
    b = rand(n,1);
    x0 = zeros(n,1);
    alpha = sqrt(gamma_low*gamma_high)
 
    [x,resvec,iter] = hss(A,b,alpha,x0,m_max,eps);
    semilogy( (0:iter),[1;resvec(1:iter)]);
    hold on
    str = [' N = ', num2str(N)];
    legend_cell{k} = str;
    k = k+1;
end

hold off
title('residual norms');
legend(legend_cell);

