% -------------------------------------------
% Primo esercizio del 13 febbario 2006
%--------------------------------------------

function lambda=metGivensSuccSturm(d,b,ind,tol)

n=length(d);

if (ind > n | ind < 1)
    disp('"ind" non corretto: 1<= ind <= length(d) ! ' );
    return
end

%-----------------------------------------------------
% Calcolo dell'intervallo [alfa,beta] da cui iniziare
% la ricerca di autovalori.
%-----------------------------------------------------

b1=[0 b 0];             %completo il vettore b
d1(1)=d(1)-abs(b1(1));  
d1(n)=d(n)-abs(b1(n-1));  
d2(1)=d(1)+abs(b1(1));  
d2(n)=d(n)+abs(b1(n-1));  
for i=2:n-1,
  d1(i)=d(i)-(abs(b1(i-1))+abs(b1(i)));  
  d2(i)=d(i)+(abs(b1(i-1))+abs(b1(i)));
end
alfa=min(d1);
beta=max(d2);

% ------------------------------------------------
% Uso la bisezione per isolare l' autovalore
% ind-esimo
%-------------------------------------------------

aa=alfa; bb=beta;

dist=abs(bb-aa);
ss=abs(aa)+abs(bb);
while (dist > tol*ss),
c=(aa+bb)/2;
[p,s]=succSturm(d,b,c);
if s>n-ind,
    bb=c;
else
    aa=c;
end
dist=abs(bb-aa);
ss=abs(aa)+abs(bb);
end
lambda=(aa+bb)/2;

%---------------------------------------------
%---------------------------------------------
function [p,cs]=succSturm(d,b,x)
%-----------------------------------------------
% Calcolo della successione di Sturm 'p' in x 
% a partire dai vettori d e b 
% e dei corrispondenti cambiamenti di segno 'cs'
%-----------------------------------------------
n=length(d);
p(1)=1; p(2)=d(1)-x;
for i=2:n,
    p(i+1)=(d(i)-x)*p(i)-b(i-1)^2*p(i-1);
end
cs=0; %contatore cambi di segno
s=0;  %contatore dei segni costanti
for i=2:length(p),
    if(p(i)*p(i-1)<=0), 
        cs=cs+1;
    end
    if(p(i)==0),
        s=s+1;
    end
end
    cs=cs-s;    
return

% -------------------------------------------
% Secondo esercizio del 13 febbario 2006
%--------------------------------------------


function [D,iter,phi1]=symJacobi(A,tol)
%----------------------------------------------------
% Data la matrice simmetrica A e una tolleranza tol,
% determina, mediante il metodo di Jacobi
% tutti i suoi autovalori che memorizza nella
% matrice diagonale D. Determina inoltre il numero
% di iterazioni, iter, e la quantit� phi che
% contiene la norma degli elementi extra-diagonali 
% e come converge il metodo.
%----------------------------------------------------
n=max(size(A)); D=A; phiD=norm(A,'fro');
%epsi=tol*phiD;
epsi=tol;
phiD=phiNorm(D);
iter=0;
[phi1]=phiD;
while (phiD > epsi)
    iter=iter+1;
    for p=1:n-1,
        for q=p+1:n
            [c,s]=calcoloCeS(D,p,q);
            D=calcoloProdottoGtDG1(D,c,s,p,q);
        end;
    end
    phiD=phiNorm(D);
    phi1=[phi1; phiD];
end
return

%---------------------------------
function [c,s]=calcoloCeS(A,p,q)
if abs(A(p,q))<1.e-15,
 c=1; s=0;
else
    eta=(A(q,q)-A(p,p))/(2*A(p,q));
    if (eta>=0),
        t=1/(eta+sqrt(1+eta^2));
    else
        t=-1/(-eta+sqrt(1+eta^2));
    end
    c=1/(sqrt(1+t^2));
    s=c*t;
end;
return

%------------------------------------
function phin=phiNorm(A)
phin=0;
n=max(size(A));
for i=1:n-1,
    for j=i+1:n,
        phin=phin+A(i,j)^2+A(j,i)^2;
    end
end
phin=sqrt(phin);
return

%--------------------------------------------

 function [D]=calcoloProdottoGtDG(D,c,s,p,q)
 %
 %  Prima eseguo Gt*D, quindi D*G.
 %
 n=max(size(D));
 
 %Costruisco G(p,q)
 G=eye(n);
 G(p,p)=c;
 G(p,q)=s; G(q,p)=-s;
 G(q,q)=c;
 
 D=(G'*D)*G;
 
 return
      
