function [T,V]=Tridiagonalize(A)
% TRIDIAGONALIZE transforms a symmetric matrix to tridiagonal form 
%   [T,V]=Tridiagonalize(A) transforms the symmetric matrix A into
%   tridiagonal form T using rotations such that T=V'*A*V

n=length(A); V=eye(n); 
for p=2:n-1
  for q=p+1:n 
    [c,s]=GivensRotation(A(p,p-1),A(q,p-1));
    if s~=0
      A(p,p-1)=A(p,p-1)*c-A(q,p-1)*s;         % one element on row p
      z=(A(p,p)-A(q,q))*s+2*A(q,p)*c;         % corner elements
      A(p,p)=A(p,p)-z*s; A(q,q)=A(q,q)+z*s;
      A(q,p)=-A(q,p)+z*c;
      h=A(p+1:q-1,p)*c-A(q,p+1:q-1)'*s;       % colum p and row q
      A(q,p+1:q-1)=A(p+1:q-1,p)*s+A(q,p+1:q-1)'*c;
      A(p+1:q-1,p)=h;
      h=A(q+1:n,p)*c-A(q+1:n,q)*s;            % rest of colums p and q
      A(q+1:n,q)=A(q+1:n,p)*s+A(q+1:n,q)*c;
      A(q+1:n,p)=h;
      h=V(:,p)*c-V(:,q)*s;                    % update V=V*U
      V(:,q)=V(:,p)*s+V(:,q)*c;
      V(:,p)=h;
    end
  end
end
c=diag(A,-1);
T=diag(diag(A))+diag(c,1)+diag(c,1)';
