//
clear;
//
Percorso = "...";
exec(Percorso + "EGPP.sci");
exec(Percorso + "SA.sci");
exec(Percorso + "SI.sci");
//
T_J = 0;
T_EG = 0;
NE = [];
//
N = 100;
ripetizioni = 50;
ITER = 35;
//
// commentare la riga seguente per n > 3000
A = 3*eye(N,N); A(2:N,1) = ones(N-1,1); A(1:N-1,N) = ones(N-1,1);
b = ones(N,1);
//
c = b/3;
//
x_esatta = ones(N,1)/6; x_esatta(1) = 1/4; x_esatta(N) = 1/4;
tic();
for j=1:ripetizioni,
    x = zeros(N,1);
    xn = x;
    for k = 1:ITER,
        s = -(x(1)+x(N))/3;
        for i = 2:N-1, xn(i) = s + c(i); end;
            xn(1) = c(1) - x(N)/3;
            xn(N) = c(N) - x(1)/3;
            x = xn;
            NE(k) = norm(x-x_esatta,1);
    end;
end;
T_J = toc()/ripetizioni;
//
printf("\nerrore relativo J = %4.3e\n",norm(x - x_esatta,1)/norm(x_esatta,1));
printf("\ntempo per Jacobi (N = %d) = %4.3e\n",N, T_J);
//
ERR_ZERO = find(NE == 0);
if ERR_ZERO ~= [] then NE = NE(1:min(ERR_ZERO)-1); end;
K = [1:length(NE)]';
scf(0);clf();
plot2d(K, [log10(3*NE(1)) - K * log10(3), log10(NE)], style = [2,5]);
xgrid();
xlabel('k (iterazione)'); legend('C'' - k * log10(3)' , 'log10( || x(k) - x* || )');
//
tic();
//for j=1:ripetizioni,
    [S,D,P] = EGPP(A);
    [c] = SA(S,P*b);
    [x_EG] = SI(D,c);
//end;
T_EG = toc();
//
printf("\nerrore relativo EG = %4.3e\n",norm(x_EG - x_esatta,1)/norm(x_esatta,1));
printf("\ntempo per Gauss (N = %d)= %4.3e\n", N, T_EG);
//
