function [Q,R]=AddColumnQR(Q,R,k,w);
% ADDCOLUMNQR_update QR decomposition when a new column is added
%   [Q,R]=AddColumnQR(Q,R,k,w); finds the QR decomposition of
%   [A(:,1),A(:,2),... A(:,k-1),w,A(:,k),...A(:,n)] when [Q,R]=qr(A).
%   Uses Matlab PLANEROT

[m,n]=size(R);
R=[R(:,1:k-1),Q'*w,R(:,k:n)];
for j=m-1:-1:k
  G=planerot(R(j:j+1,k));
  R(j:j+1,k)=G*R(j:j+1,k);                    % annihilate new column
  if j<=n,                                    % transform also R
    R(j:j+1,j+1:n+1)=G*R(j:j+1,j+1:n+1); 
  end
  Q(:,j:j+1)=Q(:,j:j+1)*G';                   % update Q
end
