%esercizio 1 02/12/22
clear all
close all
clc

a0=0;
a1=2;

[control0]=sylvestermio(a0);
[control1]=sylvestermio(a1);

if control0==1 && control1==1
    disp('La matrice non risulta definita positiva per nessuno dei due valori!')
elseif control0==0 && control1==1
    sprintf('La matrice risulta definita positiva per a=%d!',a1)
    a=a1;
elseif control0==1 && control1==0
    sprintf('La matrice risulta definita positiva per a=%d!',a0)
    a=a0;
else
    disp('La matrice risulta definita positiva per tutti e due valori!')
end

A=[2+a^2 5 6 3 1; 5 14 6 1 1; 6 6 12 7 1;3 1 7 5 1;1 1 1 1 2];
b=[1 0 0 0 0]';

[S]=CHOLESKY(A,b);

%A questo punto ho il sistema S*S'*x=b
%Passo a S*y=b e S'*x=y

%Prima trovo y e utilizzo la risuìoluzione di un sistema triangolare
%inferiore
[y]=SOLVETRIL(S,b);
[x]=SOLVETRIU(S',y);

function [c] = sylvestermio(a)
%funzione per la creazione delle sottomatrici, si può evitare e definirla
%nello script

A=[2+a^2 5 6 3 1; 5 14 6 1 1; 6 6 12 7 1;3 1 7 5 1;1 1 1 1 2];
n=length(A);


d=A(1,1);
for i=1:n
    if d>0
        Ai = A(1:i, 1:i);
        d=LUDET(Ai);
    else
        c=0;
        return
    end
end
c=1;
end

function [d] = LUDET(Ai)
%decomposizione LU con versione doolittle
n=length(Ai);

L=eye(n);
U(1,:)=Ai(1,:);

for k=1:n
    U(k,k:n)=Ai(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
    L(k+1:n,k)=(Ai(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
end
d=1;
for i=1:n
    d=d*U(i,i);
end
end


function [S] = CHOLESKY(A,b)
%Applicazione di CHoleshy
n=length(A);


for i=2:n
    for j=1:i-1
        if A(i,j)==A(j,i)
        else
            disp('Matrice non simmetrica');
        end
    end
end
%inizio fattorizzazione di cholesky
S=zeros(n);
S(1,1)=sqrt(A(1,1));
for i=2:n
    for j=1:i-1
        S(i,j)=(A(i,j) - S(i, 1:j-1)*S(j, 1:j-1)')/S(j,j);
        S(i,i)=sqrt(A(i,i)-S(i,1:i-1)*S(i,1:i-1)');
    end
end
end

function [y]=SOLVETRIL(S,b)

n=length(S);

y=zeros(n,1);
y(1)=b(1)/S(1,1);
for i=1:n
    y(i)= (b(i) - S(i,1:i-1)*y(1:i-1))/S(i,i);
end
end

function [x]=SOLVETRIU(H,y)

n=length(H);
x=zeros(n,1);
x(n)=y(n)/H(n,n);

for i=n-1:-1:1
    x(i)=(y(i) - H(i,i+1:n)*x(i+1:n))/H(i,i) ;
end
end