function [Q,R] = UpdateQR(Qc,Rc,u,v)
% Updating the QR factorization of a matrix that has be changed by a
% rank-1 matrix.
% Qc is nxn orthogonal and Rc is nxn upper triangular.
% u and v are column n-vectors.
% Q is orthogonal and R is upper triangular so QR = Qc*Rc + uv'
[n,n] = size(Qc);
% Generate and apply Givens rotations G{n-1},...,G{2},G{1} so that
% G{1}'*...*G{n-1}'u is zero except in the first component.
% Rc and Qc must be updated as follows..
% Rc = G{1}'*...*G{n-1}'*Rc
% Qc = Qc*G{n-1}*...G{1}
u = Qc'*u;
for k=n-1:-1:1
[c,s] = Givens(u(k),u(k+1));
Rc(k:k+1,k:n) = [c s;-s c]'*Rc(k:k+1,k:n);
u(k:k+1) = [c s;-s c]'*u(k:k+1);
Qc(:,k:k+1) =Qc(:,k:k+1)*[c s;-s c];
end
% Rc is now upper hessenberg. The same can be said of Rc + u*v'
Rc(1,:) = Rc(1,:) + u(1)*v';
% Now generate more rotations G{1},...,G{n-1} so that
% G{n-1}'*...*G{1}'*R is upper triangular. Update Qc along the way.
for k=1:n-1
[c,s] = Givens(Rc(k,k),Rc(k+1,k));
Rc(k:k+1,k:n) = [c s;-s c]'*Rc(k:k+1,k:n);
Qc(:,k:k+1) = Qc(:,k:k+1)*[c s;-s c];
end
R = Rc;
Q = Qc;

function [c,s] = Givens(a,b)
% function [c,s] = Givens(a,b)
% Givens rotation computation
% Determines cosine-sine pair (c,s) so that [c s;-s c]'*[a;b] = [r;0]
if b==0
c = 1; s = 0;
else
if abs(b)>abs(a)
tau = -a/b; s = 1/sqrt(1+tau^2); c = s*tau;
else
tau = -b/a; c = 1/sqrt(1+tau^2); s = c*tau;
end
end