function varargout = impdsqr(H)
%IMPDSQR One step of the Implicit Double Shift QR Iteration
%
%   [Q,H] = impdsqr(H) returns modified upper Hessenberg matrix H,
%   using one step of implicit double shift QR iteration.
%   Q is a product of Householder reflections such that QHQ' = H2.
%   H2 is the result of a double shift using the eigenvalues of
%   the lower right-hand 2 x 2 matrix H(n-1:n,n-1:n).

[m,n] = size(H);
if m~=n
    disp('matrix H  is not square')  ;
    return;
end;

if nargout == 2
    Q = eye(n,n);
end

% trace of lower 2 x 2 matrix
trce = H(n-1,n-1) + H(n,n);
%det of lower 2 x 2 matrix
determ = H(n-1,n-1)*H(n,n) - H(n,n-1)*H(n-1,n);
% first nonzero entries of column 1 of S
x = H(1,1)^2 + H(1,2)*H(2,1) - H(1,1)*trce + determ;
y = H(1,1)*H(2,1) + H(2,2)*H(2,1) - H(2,1)*trce;
z = H(3,2)*H(2,1);

for i = 0:n-3
   [u,beta] = houseparms([x y z]');
   lengu = length(u);
   Hu_bar = eye(lengu, lengu) - beta* u * u';
   H(i+1:i+3,1:n) = Hu_bar*H(i+1:i+3,1:n);
   H(1:n,i+1:i+3) = H(1:n,i+1:i+3)*Hu_bar;
   if nargout == 2
      P = eye(n);
      P(i+1:i+3,i+1:i+3) = Hu_bar;
      Q = P*Q;
   end

   x = H(i+2,i+1);
   y = H(i+3,i+1);
   if i < n-3
      z = H(i+4,i+1);
    end;
end

[u,beta] = houseparms([x y ]');
lengu = length(u);
Hu_bar = eye(lengu, lengu) - beta * u * u';
H(n-1:n,1:n) = Hu_bar*H(n-1:n,1:n);
H(1:n,n-1:n) = H(1:n,n-1:n)*Hu_bar;
if nargout == 2
   P = eye(n);
   P(n-1:n,n-1:n) = Hu_bar;
   Q = P*Q;
end

if nargout == 2
   varargout{1} = Q';
   varargout{2} = H;
else
    varargout{1} = H;
end
