%es1 05.09.22

clear all
close all
clc

% PRIMO PUNTO

% 1. Scrivo funzione che mi costruisce A in base ad alpha
% 2. Essendo alpha solo in posizione A(1,1) verifico la simmetria una sola
%    volta all'inizio del ciclo;
% 3. Chiamo fattorizzazione e la rendo tridiagonale;
% 4. Costruisco vettore polinomi di STURM;
% 5. Pongo gamma = 0 e verifico # var. di segno --> trovo 1 alpha ed esco;

b=[3 14 9 5 3]';
alphaV=linspace(0,3);
[a,A] = Ai(alphaV);

% SECONDO PUNTO

% 1. Metodo grad e grad coniugato
Nmax=100;
toll=sqrt(eps);
[x0,e0,k0] = Grad(A,b,Nmax,toll);
[x1,e1,k1] = Gradcon(A,b,Nmax,toll);

% TERZO PUNTO 

% 1. Grafico errore/residuo entrambi (in scala logaritmica)
K0=linspace(1,k0,k0);
K1=linspace(1,k1,k1);

% figure 
% plot(K0,log10(e0),'r')
% hold on
% plot(K1,log10(e1),'b');
% tilte ('Convegenza in base all''andamento dell''errore');

figure 
% plot(K0,log10(e0),'r')
% hold on
plot(K1,e1,'b');
title ('Convegenza in base all''andamento dell''errore');

function [a,A] = Ai(alphaV)
q=0;
for av=alphaV
    A=[av 0 1 0 0; 0 10 2 1 1; 1 2 6 0 0; 0 1 0 4 0; 0 1 0 0 2];
    n=length(A);
    if q==0
        k=0;
        for i=2:n
            for j=i-1:n
                if A(i,j)==A(j,i)
                    k=k+1;
                end
            end
        end

        if k~=(n^2/2 +n/2 -1)
            disp('matrice non simmetrica');
            return
        end
        q=1;
    end

    [H] = QRH(A);
    [esito] = PSturm(H);

    if esito
        a=av;
        return
    end

end
end

function [H] = QRH(A)
%matrice di hess
n=length(A);
Q=eye(n);

for i=1:n-2
    v=A(i+1:n,i);
    u=zeros(n,1);

    if v(1)==0
        u(i+1:n)= v + norm(v)*eye(n-i,1);
    else
        u(i+1:n)= v + sign(v(1))*norm(v)*eye(n-i,1);
    end

    Qi=eye(n) -2*(u*u')/(u'*u);
    Q=Qi*Q;
    A=Qi*A*Qi;
end
H=A;
end

function [esito] = PSturm(H)

n=length(H);
C=diag(H);
D=diag(H,1);
p(1)=1;
p(2)=-C(1);
esito=0;
r=0;
for i=2:n
    if p(i)==0
        return
    else
        p(i+1)= -C(i)*p(i) - D(i-1)^2*p(i-1);
        if p(i)*p(i-1)<0
            r=r+1;
        end
    end
end
if p(n+1)*p(n)<0 && p(n+1)~=0
    r=r+1;
end
if r==n
    esito=1;
end
end

function [x0,e0,k0] = Grad(A,b,Nmax,toll)
n=length(b);
x0=zeros(n,1);
r0=A*x0 -b;
a0=r0'*r0/(r0'*A*r0);
k0=1;
e0(k0)=norm(r0);

while k0<Nmax && e0(k0)>toll
    x=x0 -a0*r0;
    r0=A*x-b;
    k0=k0+1;
    e0(k0)=norm(r0);
    x0=x;
end

end

function [x1,e1,k1] = Gradcon(A,b,Nmax,toll)
n=length(b);
x1=zeros(n,1);
r1=A*x1 -b;
d1=-r1;
k1=1;
e1(k1)=norm(r1);

while k1<Nmax && e1(k1)>toll
    a1=r1'*r1/(r1'*A*r1);
    x=x1 +a1*d1;
    r1=A*x1-b;
    b1=(d1'*A*r1)/(d1'*A*d1);
    d=-r1+b1*d1;
    k1=k1+1;
    e1(k1)=norm(r1);
    x1=x;
    d1=d;
end
end