function [H,P]=qrt(H,P)
%
%  [H,P]=qrt(H,P)
%
%     Copyright: Eduardo Casas (Octubre, 1999)
% 
%  Esta función realiza una iteracción del algoritmo QR aplicado a
%  una matriz simétrica y tridiagonal H, obteniendo a su vez otra
%  matriz simétrica y tridiagonal semejante a la primera. Además si
%  el número de argumentos de salida es 2 y solamente se suministra
%  H como argumento de entrada entonces la salida P es la matriz
%  utilizada en la iteración: P'*H*P es la matriz de salida.
%  Si el número de argumentos de entrada y salida es 2, entonces la
%  matriz de salida P es el producto (por la derecha) de la de entrada
%  y la matriz utilizada en la iteración QR.
%
n=size(H,1);
if nargout == 2 & nargin == 1
   P=eye(n,n);
end
c=[];s=[];n1=n-1;
for j=1:n1
   j1=j+1;
   if H(j1,j) == 0
      c=[c;1];
      s=[s;0];
   else
      nx=norm(H(j:j1,j));
      c1=H(j,j)/nx;
      s1=-H(j1,j)/nx;
      c=[c;c1];
      s=[s;s1];
      H(j,j)=nx;
      H(j1,j)=0;
      if j == n1
         j2=n;
      else
         j2=j1+1;
      end
      H(j:j1,j1:j2)=[c1 -s1;s1 c1]*H(j:j1,j1:j2);
   end
end
for j=1:n1
   if c(j) ~= 1 | s(j) ~= 0 
      j1=j+1;
      if j > 1
         j0=j-1;
         H(j0,j)=H(j,j0);
         H(j0,j1)=0;
      end
      H(j,j)=c(j)*H(j,j)-s(j)*H(j,j1);
      H(j1,j)=-s(j)*H(j1,j1);
      H(j,j1)=H(j1,j);
      H(j1,j1)=c(j)*H(j1,j1);
      if nargout == 2
         P(:,j:j1)=P(:,j:j1)*[c(j) s(j);-s(j) c(j)];
      end
   end
end