% cdupdate.m
%
% Updates the centroid decomposition of a matrix H when appending a
% row or column.
%
% Jason R. Blevins <jrb11@duke.edu>
%
% $Id: cdupdate.m,v 1.11 2004/08/04 21:32:25 jrblevin Exp $
%
% Usage:
%
%   [B1,V1,Z1,count,time,Count,Time] = cdupdate1(B0,V0,Z0,a,b,k)
%
% Parameters:
%
%   B0,V0,Z0 - the centroid decomposition of H0 with sign vectors Z0
%   a,b      - column vectors which form the rank-1 update a+b'
%   k        - the desired truncation rank (or full rank if k==0)
%
% Returns:
%
%   B1,V1,Z1 - the updated CD, B1*V1', with sign vectors Z1
%   count    - the total number of sign vector search iterations
%   time     - the total time taken to find the CD
%   Count    - a vector containing
%   Time     - a vector containing

function [B1,V1,Z1,count,time,Count,Time] = cdupdate(B0,V0,Z0,a,b,k_new);

% Initialize Variables
epsilon = 10^(-13);

% Determine matrix dimensions
[m_b, k_b] = size(B0);
[l_v, k_v] = size(V0);
l_a = length(a);
l_b = length(b);

m = m_b;
l = l_v;
k_old = k_b;

if (k_new == 0)
  k_new = rank(B0);
end

% Sanity Checks:
if (k_b ~= k_v)
  error('cdupdate.m: B0 and V0 must have the same number of columns.');
elseif (l_a ~= m_b)
  error('cdupdate.m: length of a and the number of rows of B0 must match.');
elseif (l_b ~= l_v)
  error('cdupdate.m: length of b and the number of rows of V0 must match.');
%elseif (k_new > k_old)
%  error('cdupdate.m: requested rank larger than original CD rank plus one.');
end

% Begin the timer
t0 = clock;

% Calculate necessary components of a and b
n = V0'*b;
q = b - V0*n;
qnorm = norm(q);

% Determine if we need an additional basis vector Q, then construct S
if (qnorm < epsilon)
  Q = [];
  S = [B0] + a*[V0'*b]';
else
  Q = q/qnorm;
  S = [B0, zeros(m,1)] + a*[V0'*b; Q'*b]';
end

% Temporary change to test effectiveness of S
%S = B0*V0'+a*b';

% Calculate the CD of S
%[Bs,Vs,Zs,count_s,time_s,Count_s,Time_s] = centroidsign(S,Z0,k_new);
[Bs,Vs,Zs,count_s,time_s,Count_s,Time_s] = centroidsign(S,ones(size(Z0)),k_new);

% Return values
B1 = Bs;
% Temporary change to test effectiveness of S
V1 = [V0,Q]*Vs;
%V1 = Vs;
Z1 = Zs;


count = count_s;
time = etime(clock,t0);
Count = Count_s;
Time = Time_s;