function x = sherman_mw(L,U,b,Z,V,P)

% function x = sherman_mw(L,U,b,Z,V,P)
%
% This function uses the Sherman-Morrison-Woodbury formula
% to solve the linear system
%   (A - Z V^T) x = b
% where A has been factored as  
%    PA = LU
% where  P  is a permutation matrix and linear systems involving  
% L  and  U  are easy to solve.  ([L,U,P]=lu(A) is one way to
% compute these factors.  We could also use [L,U]=qr(A) and set 
% P = I.)
% Both A and A - Z V^T are assumed to be nonsingular.
% 
% The Sherman-Morrison-Woodbury formula is
%   inv(A - Z V^T) = inv(A) + inv(A) Z inv(I-V^T inv(A) Z) V^T inv(A),
% and we apply this to b without forming inverses.
%
% Dianne P. O'Leary  12/2005


[n,k]=size(Z);

y  = U \ (L \ (P*b));
Zh = U \ (L \ (P*Z));
t  = (eye(k) - V'*Zh) \ (V'*y);
x  = y + Zh*t;