function [Qs,Rs]=UpdateQR(Q,R,u,v)
% UPDATEQR Rank-1 update of the QR-Decomposition
%   [Qs,Rs]=UpdateQR(Q,R,u,v); If As=A+u v' and [Q,R]=qr(A) then we
%   compute Qs and Rs such that As=Qs Rs. Uses Matlab's PLANEROT.

[m,n]=size(R); Qs=Q; Rs=R; w=Q'*u;
for k=m:-1:n+2                        % annihilate w(n+2:m)
  G=planerot(w(k-1:k));
  w(k-1:k)=G*w(k-1:k);
  Qs(:,k-1:k)=Qs(:,k-1:k)*G';  
end
for k=n+1:-1:2                        % annihilate w(2:n+1)
  G=planerot(w(k-1:k));
  w(k-1:k)=G*w(k-1:k);         
  Rs(k-1:k,k-1:n)=G*Rs(k-1:k,k-1:n);
  Qs(:,k-1:k)=Qs(:,k-1:k)*G';       
end
Rs(1,:)=Rs(1,:)+w(1)*v';              % Add rank-1 change to first row 
for k=1:n                             % reduce Hessenberg matrix to Rs
  G=planerot(Rs(k:k+1,k));   
  Rs(k:k+1,k:n)=G*Rs(k:k+1,k:n);    
  Qs(:,k:k+1)=Qs(:,k:k+1)*G';       
end
