%es 1 19.09.2022
clear all
close all
clc

A=[0 0 1 0 0; 0 10 2 1 1; 1 2 6 0 0; 0 1 0 4 0; 0 1 0 0 2];
k=0;
n=length(A);

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

if k~= n^2/2 -n/2
    disp('Matrice non simmetrica');
else
    [H]=HessTri(A);
    [r]=Pstrum(H);
    
    fprintf ('Ho %d autovalori positivi\n', r);
    Nmax=100;
    toll=sqrt(eps);
    l0=4;
    [l] = MetodoPot(H,l0,Nmax,toll);

end

function [H]=HessTri(A)
n=length(A);

Q=eye(n);

for i=1:n-1
    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 [r]=Pstrum(H)

D=diag(H);
C=diag(H,1);
n=length(D);
p(1)=1;
p(2)= -D(1);
K=[p(1)];

k=1;

for i=2:n
    p(i+1)= -D(i)*p(i) - C(i-1)^2*p(i-1);

    if p(i)~=0
        K=[p(i) K];
    end
  
end

if p(n+1)~=0
     K=[p(n+1) K];
end

k=length(K);
r=0;

for j=2:k
    if K(j)*K(j-1)<0
        r=r+1;
    end
end

end

function [h] = MetodoPot(H,p,Nmax,toll)
%metodo delle potenze con shift
y=ones(length(H),1);
T=H-p*eye(length(H));
k=0;
l=1;
e=1;
[L,U]=lu(T);

while e>toll && k<Nmax
    k=k+1;
    l0=l;
    [w]=Solve(L,U,y)
    [~,m]=max(w);
    l=y(m)/w(m);
    y= w/norm(w,inf)
    e = abs(l-l0)/abs(l);
end
h= l+p;
end

function [L,U]=lu(T)
%fattorizzazione LU doolittle
n=length(T);

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

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

function [w]=Solve(L,U,y)
%Risoluzione sistema triang inf e triang sup;
n=length(y);
t=eye(n,1);
t(1) = y(1)/L(1,1);

for i=2:n
    t(i)= (y(i)- L(i, 1:i-1)*t(1:i-1))/L(i,i);
end
w=eye(n,1);
w(n) = t(n)/U(n,n);

for j=n-1:-1:1
    w(j)= (t(j)- U(j, j+1:n)*w(j+1:n))/U(j,j);
end
end