function [e,q]=QRStep(l,n,e,q)
% QRSTEP compute QR iteration step on a bidiagonal matrix 
%   [e,q]=QRStep(l,n,e,q) computes one QR iteration step on the
%   bidiagonal matrix for the trailing pricipal submatrix (l:n,l:n)

an=e(n)^2+q(n)^2;                       % compute Wilkinson shift
an1=e(n-1)^2+q(n-1)^2;                  % with code from ImplicitQR.m
bn1=q(n-1)*e(n);
delta =(an-an1)/2;
if delta==0,
  sigma=an-abs(bn1);
else
  sigma=an-bn1^2/(delta+sign(delta)*sqrt(delta^2+bn1^2));
end
[c s]=GivensRotation(q(l)^2-sigma,q(l)*e(l+1));     
h=q(l)*c-e(l+1)*s;                      % first transformation from the
e(l+1)=q(l)*s+e(l+1)*c; q(l)=h;         % right generates bulge
bulge=-q(l+1)*s; q(l+1)=q(l+1)*c; 
for i=l:n-2                             % Chasing the bulge
  [c s]=GivensRotation(q(i),bulge);     % from left
  q(i)=q(i)*c-bulge*s; h=e(i+1)*c-q(i+1)*s; 
  q(i+1)=e(i+1)*s+q(i+1)*c; e(i+1)=h; 
  bulge=-e(i+2)*s; e(i+2)=e(i+2)*c;
  [c s]=GivensRotation(e(i+1),bulge);   % from right
  e(i+1)=e(i+1)*c-bulge*s; 
  h=q(i+1)*c-e(i+2)*s; e(i+2)=q(i+1)*s+e(i+2)*c; 
  q(i+1)=h; 
  bulge=-q(i+2)*s; q(i+2)=q(i+2)*c;
end
[c s]=GivensRotation(q(n-1),bulge);     % last from left
q(n-1)=q(n-1)*c-bulge*s; h=e(n)*c-q(n)*s; 
q(n)=e(n)*s+q(n)*c; e(n)=h; 
