%es 2 21.01.20
clear all
close all
clc

alphav=linspace(-12,-8);

b=[9 ; -16; -34];
[a,A] = Aa(alphav);

toll=sqrt(eps);
Nmax=100;
[xg] = GS(A,b,toll,Nmax);
[J,h]=Ja(A,b);
[xj] = JACOBI(J,h,toll,Nmax)

function [a,Ai] = Aa(alphav)
l0=0;
a=alphav(1);
Nmax=100;
toll=sqrt(eps);

for av=alphav
    A=[4 1 1; 2 -9 0; 0 av -6];
    B=A'*A;
    [l]=potenze(B,Nmax,toll);
    if l>l0
        l0=l;
        a=av;
        Ai=A;
    end
end
end

function [l] = potenze(Ai,Nmax,toll)
k=1;
m=1;
l1=1;
y=ones(size(Ai,1),1);
flag=1;

while k<Nmax && flag
    k=k+1;
    w=Ai*y;
    [~,m]=max(abs(w));
    l=w(m)/y(m);
    y=w/norm(w,inf);
    flag = (abs(l-l1))/abs(l)>toll;
    l1=l;
end
end

function [J,h] = Ja(A,b)

J= -diag(1./diag(A))*(A-diag(diag(A)));
h=diag(1./diag(A))*b;

end


function [xj] = JACOBI(J,h,toll,Nmax)
x0=zeros(length(h));
xj=h;
e=1;
k=1;

while k<Nmax && e>toll
    xj= J*xj +h;
    e=norm(xj-x0)/norm(xj);
    k=k+1;
end


end

function [xg] = GS(A,b,toll,Nmax)

L=tril(A);
U=A-L;
G=-L\U;
hg= L\b;

[esitoG]=RHO(G);

if ~esitoG
    disp('Gauss-Seidel non convergente!');
    return
end
k=1;
x0=zeros(length(b));
xg=hg;

e=1;
while k<Nmax && e>toll
    xg= G*xg +hg;
    e=norm(xg-x0)/norm(xg);
    k=k+1;
end

end

function [esito]=RHO(T)
Nmax=100;
toll=sqrt(eps);

k=1;
m=1;
l1=1;
l=1;
y=ones(size(T,1),1);
e=1;
esito=0;

while k<Nmax && e>toll
    l1=l;
    k=k+1;
    w=T*y;
    [~,m]=max(abs(w));
    l=w(m)/y(m);
    y=w/norm(w,inf);
    e = (abs(l-l1))/abs(l);
end

if abs(l1)<1
    esito=1;
end

end