function P2()
% Checks out BlockSolve

% Example 1
% A small well conditioned example that we can solve via explicit set-up...
clc
randn('seed',0)
n = 100; A = randn(n,n); A = A'*A + 10*n*eye(n,n);
q = 5; Y = randn(n,q); Z = randn(n,q);
N = (q+1)*n;
b = randn(N,1);
T = zeros(N,N);
for k=1:q+1
    rows = 1+(k-1)*n:k*n;
    T(rows,rows) = A;
    if k<=q
        T(rows,rows+n) = Y(:,k)*Z(:,k)';
        T(rows+n,rows) = Z(:,k)*Y(:,k)';
    end
end
x0 = T\b;
A = sparse(A);
xCVL = BlkSolveCVL(A,Y,Z,b);
xYOU = BlkSolve(A,Y,Z,b);
fprintf('Example1_ErrorCVL = %10.3e\n', norm(xCVL-x0));
fprintf('Example1_ErrorYOU = %10.3e\n\n\n', norm(xYOU-x0));

% Example 2
iProb = 39;
iProb = 4;
Problem = UFget(iProb);
A = 100*Problem.A;
[n,n] = size(A);

fprintf('Example 2:  UFget(%1d), n = %5d, nnz(A) = %10d\n\n',iProb,n,nnz(A))
disp('   q     tCVL    tYOU     diffx')
disp('------------------------------------')
for q = [25 50 100 200]
   Y = randn(n,q)/100;
   Z = randn(n,q)/100;
   b = rand(n*(q+1),1);
   tic,xCVL = BlkSolveCVL(A,Y,Z,b); tCVL = toc;
   tic,xYOU = BlkSolve(A,Y,Z,b); tYOU = toc;
   diffx = norm(xCVL-xYOU)/norm(xCVL);
   fprintf('%5d  %6.2f  %6.2f  %10.3e\n',q,tCVL,tYOU,diffx)
end



 function x = BlkSolveCVL(A,Y,Z,b)
% A is n-by-n, symmetric,positive definite, in sparse form
% Y and Z are n-by-q in full format
% b is column n(q+1)-by-1 in full format
% x is a column n(q+1)-by-1 that solves Tx = b where T is a (q+1)-by(q+1)
% Block tridiagonal matrix where each diagonal block is A, block (i,i+1)
% is Y(:,i)Z(:,i)', and block (i+1,i) is Z(:,i)Y(:,i)'

[n,q] = size(Y);
N = n*(q+1);

% Note that T = S + UMU' where U = [Ymat Zmat] with
%
%    S = kron(eye(q+1),A), 
%                       
%    M =  [ 0 I ; I 0 ]                   (q-by-q blocks)
%
%            Y(:,1)    0     0           
%    Ymat =     0    Y(:,2)   0           (q = 3 case, all blocks n-by-1)
%               0       0   Y(:,3)
%               0       0     0
%                  
%               0       0      0
%    Zmat =   Z(:,1)    0      0          (q = 3 case, all blocks n-by-1)
%               0     Z(:,2)   0
%               0       0    Z(:,3)

% If (S + UMU')x = b, then (I + WMU')x = z where SW = U and Sz = b.

% Using Sherman-Morrison-Woodbury and the fact that inv(M) = M we have
%
%         x = inv(I + WMU')z 
%           = (I - W*M*inv(I+U'*WM)U')z
%           = z - W*inv(M + U'*W)*U'*z

% We now evaluate the right hand side,
%
%           Compute Cholesky PAP' = GG'
%           Solve Sz = b where S = kron(eye(q+1,q+1),A)
%           Solve SW = U
%           z1 = U'z
%           z2 = (M + U'*W)\z1
%           z3 = W*z2;
%           x  =  z - z3
%
% exploiting structure when it makes sense. In particular,
% we never form S), U or W explictly. In terms of memory these
% matrices require nnz(A), nq, and nq storage locations to represent. 
% Not ((q+1)n)^2, (q+1)nq, (q+1)nq).

% Compute PAP' = GG' via ...

[G,flag,p] = chol(A,'lower','vector');  % Note: P = I(:,p)
    
%G = chol(A)';

% Solve Sz = b...
%z = reshape(G'\(G\reshape(b,n,q+1)),N,1);
b = reshape(b,n,q+1);
b(p,:) = G'\(G\b(p,:));
z = reshape(b,N,1);

% Compute z1 = U'z = [Ymat Zmat]'z = [Ymat'*z ; Zmat'*z]

z1 = zeros(2*q,1);
for k=1:q
    rows    = 1+(k-1)*n:k*n;
    z1(k)   = Y(:,k)'*z(rows);
    z1(k+q) = Z(:,k)'*z(rows+n);
end

% Next we need W, but we do not form it explicitly.
% Solve A*Ytilde = Y and A*Ztilde = Z ...
Y_tilde = zeros(n,q);
Z_tilde = zeros(n,q);
Y_tilde(p,:) = G'\(G\Y(p,:));
Z_tilde(p,:) = G'\(G\Z(p,:));

% Note that W = [Ymat_tilde Zmat_tilde] where the recipes for Ymat_tilde
% and Zmat_tilde are the same as the recipes for Ymat and Zmat with Y and
% Z replaced by Ytilde and Ztilde.

% Compute the 2q-by-2q matrix 
%                       Ymat'*Ymat_tilde             eye(q,q)+ Ymat'*Zmat_tilde
%    C = M + U'*W   = 
%                       eye(q,q) + Zmat_tilde'*Ymat       Zmat'*Zmat_tilde

% The diagonal blocks are diagonal.
C11 = zeros(q,q); C22 = zeros(q,q);
for k=1:q
    C11(k,k) = Y(:,k)'*Y_tilde(:,k);
    C22(k,k) = Z(:,k)'*Z_tilde(:,k);
end
% The (1,2) block is lower bidiagonal with ones on the diagonal...
C12 = eye(q,q);
for k=1:q-1
    C12(k+1,k) = Y(:,k+1)'*Z_tilde(:,k);
end
% C is symmetric...
C = [C11 C12; C12' C22];

% Solve C*z2 = z1. This is a small system so it is OK just to use \
z2 = C\z1;

% Compute z3 = W*z2 = [Ymat_tilde Zmat_tilde]*z2 
%                   = Ymat_tilde*z2(1:q) + Zmat_tilde*z2(q+1:2*q)
z3 = zeros(N,1);
for k=1:q
    rows = 1+(k-1)*n:k*n;
    z3(rows)   = z3(rows)   + Y_tilde(:,k)*z2(k);
    z3(rows+n) = z3(rows+n) + Z_tilde(:,k)*z2(k+q);
end
x = z - z3;

% -2 if your code is correct but very slow. This usually happened
% because you mis-manaqed
% the use of Y and Z
