function TestSemiSep()
% Checks out SemiSep
clc
randn('seed',0)

% A small example to check for correctness...
n = 100;
d = randn(n,1); c = randn(n,1); u = randn(n,1); v = randn(n,1); x = randn(n,1);
y0 = SemiSep(d,c,u,v,x,4);
y  = SemiSepCVL(d,c,u,v,x,4);
Example1Error = norm(y - y0)

% A large example to check for efficiency...
n = 100000;
d = randn(n,1); c = randn(n,1); u = randn(n,1); v = randn(n,1); x = randn(n,1);
disp(' ')
disp('  nmin      tCVL    tYou')
disp('----------------------------')
for nmin = [4 8 16 32 64 128]
   tic
   y = SemiSepCVL(d,c,u,v,x,nmin);
   tCVL = toc;
   tic
   y = SemiSep(d,c,u,v,x,nmin);
   tYou = toc;
   fprintf('%5d     %6.3f   %6.3f\n',nmin,tCVL,tYou)
end

function y = SemiSepCVL(d,c,u,v,x,nmin)
n = length(x);
if n<=nmin
    % Explicitly form A...
    A = zeros(n,n);
    for i=1:n
      for j=1:n
        if i==j
            A(i,j) = d(i)+u(i)*v(i);
        elseif j==i+1
            A(i,j) = c(i+1)+u(i)*v(j);
        elseif i==j+1
            A(i,j) = c(i) + v(i)*u(j);
        elseif j>i+1
            A(i,j) = u(i)*v(j);
        else
            A(i,j) = v(i)*u(j);
        end
      end
    end
    y = A*x;
else
    % Split the problem in half...
    n1= floor(n/2); n2 = n-n1;  
    T = 1:n1; B = n1+1:n;
    % yT      A(T,T) A(T,B)   x(T)
    %     =                 *
    % yT      A(B,T) A(B,B)   x(B)
    zT = SemiSepCVL(d(T),c(T),u(T),v(T),x(T),nmin);  % A(T,T)*x(T)
    yT = zT + u(T)*(v(B)'*x(B)); yT(n1) = yT(n1)+ c(n1+1)*x(n1+1); 
    zB = SemiSepCVL(d(B),c(B),u(B),v(B),x(B),nmin);  % A(B,B)*x(B)
    yB = zB + v(B)*(u(T)'*x(T)); yB(1)  = yB(1) + c(n1+1)*x(n1);
    y = [yT;yB];
end
   
    
    
