% given an m-by-n matrix A and a nonzero n-vector v with v(1) = 1, the
% following algorithm returns AP where P = I - 2*v*v'/(v'*v)
% (this col.house(A,v) in Golub and Van Loan, page 197)

function A = righthouse(A,v)

beta = -2/(v'*v);
w = beta*A*v;
A = A + w*v';