function H=dobleqr(H)
%
%   H=dobleqr(H)
%
%     Copyright: Eduardo Casas (Octubre, 1999)
%
%   Esta función realiza una iteración implícita del método QR con
%   doble shift de Wilkinson aplicada a una matriz Hessenber superior H,
%   obteniendo de nuevo una matriz Hessenberg superior semejante a la
%   primera.
%
n=size(H,1);
if n == 2, return,end
n1=n-1;
t=H(n1,n1)+H(n,n);
s=H(n1,n1)*H(n,n)-H(n1,n)*H(n,n1);
u=[(H(1,1)-t)*H(1,1)+s+H(1,2)*H(2,1);H(2,1)*(H(1,1)+H(2,2)-t);H(2,1)*H(3,2)];
s=norm(u);
if u(1) < 0; s=-s;end
u(1)=u(1)+s;
t=2/(u'*u);
p=eye(3,3)-t*u*u';
H(1:3,:)=p*H(1:3,:);
if n == 3
   H=H*p;
else
   H(1:4,1:3)=H(1:4,1:3)*p;
end
for j=1:n-3
   j1=j+1;
   j3=j+3;
   u=H(j1:j3,j);
   s=norm(u);
   if u(1) < 0; s=-s;end
   u(1)=u(1)+s;
   t=2/(u'*u);
   p=eye(3,3)-t*u*u';
   H(j1:j3,j)=[-s;0;0];
   H(j1:j3,j1:n)=p*H(j1:j3,j1:n);
   jn=min([j3+1;n]);
   H(1:jn,j1:j3)=H(1:jn,j1:j3)*p;
end
n2=n-2;
u=H(n1:n,n2);
s=norm(u);
if u(1) < 0;s=-s;end
u(1)=u(1)+s;
t=2/(u'*u);
p=eye(2,2)-t*u*u';
H(n1:n,n2)=[-s;0];
H(n1:n,n1:n)=p*H(n1:n,n1:n);
H(1:n,n1:n)=H(1:n,n1:n)*p;