function [aa,bb]=OrthogonalQDStepShifted(a,b,s)
% ORTHOGONALQDSTEPSHIFTED computes one orthogonal QD step with shift
%   [aa,bb]=OrthogonalQDStepShifted(a,b,s) computes the LR-step
%   R*R^T=Rt^T*Rt+s*I with R = diag(a)+diag(b,1) and shift s>0.  The
%   result is Rt=diag(aa)+diag(bb)

n=length(a); 
sig=sqrt(s);
for k=1:n-1
  [si,co,y]=GeneralizedGivensRotation([a(k);0], sig);
  a(k)=y(1);
  [co,si]=GivensRotation(a(k),b(k));
  a(k)=co*a(k)-si*b(k);
  b(k)=-si*a(k+1);
  a(k+1)=co*a(k+1);
end
[si,co,y]=GeneralizedGivensRotation([a(n);0],sig);
a(n)=y(1);
aa=a; bb=b;
