%esercizio 1 esame 28.04.22

clear all
close all
clc

A=[10 0 1 0; 0 5 2 -1; 1 2 8 2; 0 -1 2 7];
b=[11 6 13 8]';

n=length(A);
q=0;

for i=2:n 
    for j=1:i-1
        if A(i,j)==A(j,i)
            q=q+1;
        end
    end
end

if q~=n^2/2 +n/2 -n 
    disp('Matrice non simmetrica')
    return
end
[H] = QRA(A);
[esito] = Sturm(H);

if ~esito
    disp('Matrice non definita positiva');
    return
end
Nmax=100;
toll=sqrt(eps);
[x] = Gradiente(A,b,Nmax,toll);
[x] = Gradiente(A,b,Nmax,toll);

function [H] = QRA(A)
%porto a matrice tridiagonale
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] = Sturm(H)
%verifica che tutti gli autovalori siano maggiori di 0

D=diag(H);
C=diag(H,1);

n=length(D);
p(1)=1;
p(2)= -D(1);
esito=0;
q=0;
for i=2:n
    if p(i)==0
        return
    else
        p(i+1)= -D(i)*p(i)- C(i-1).^2*p(i-1);
        if p(i)*p(i-1)<0
            q=q+1;
        end

    end
end

if p(n+1)*p(n)<0 && q==n-1
    esito=1;
end
end

function [x] = Gradiente(A,b,Nmax,toll)
%Metodo del gradiente
n=length(b);
x0=zeros(n,1);
r=A*x0-b;
k=1;
e(1)=norm(r);

while k<Nmax && e(k)>toll
    a=(r'*r)/(r'*A*r);
    x =x0 - a*r;
    r= A*x-b;
    k=k+1;
    e(k)=norm(r);
    x0=x;

end

end

function [x0] = GradienteC(A,b,Nmax,toll)
%Metodo del gradiente coniugato
n=length(b);
x0=zeros(n,1);
r=A*x0-b;
d=-r;
k=1;
e(1)=norm(r);

while k<Nmax && e(k)>toll
    a=(r'*r)/(r'*A*r);
    x =x0 - a*r;
    r= A*x-b;
    b=d0'*A*d0/(d0'*A*d0);
    d= -r +b*d0;
    k=k+1;
    e(k)=norm(r);
    x0=x;
    d0=d;
end
end