function P1Grade()
% Checks out SymSol
clc
% A very small example to check for correctness..
n = 5;
L = tril(randn(n,n),-1)+eye(n,n);
B = randn(n,n); B = B+B';
X0 = reshape((kron(L,L)+eye(n^2,n^2))\B(:),n,n);
disp('True Solution = ')
disp(X0)
X = SymSolve(L,B);
disp('SymSolve(L,B) = ')
disp(X)
relErr = norm(X-X0,'fro')/norm(X0,'fro');
fprintf('||X0 - X||/||X0|| = %10.3e\n',relErr)

% A large example for benchmarking..
n = 500;
L = tril(randn(n,n),-1)+eye(n,n);
B = randn(n,n); B = B+B';
tic
X = SymSolve(L,B);
tYOU = toc;
tic
X = SymSolveCVL(L,B);
tCVL = toc;
tic
X = SymSolveCVLSym(L,B);
tCVLSym = toc;
fprintf('\n\ntYOU/tCVL = %5.2f\n',tYOU/tCVL)
fprintf('\n\ntCVLSym/tCVL = %5.2f\n',tCVLSym/tCVL)

function X = SymSolveCVL(L,B)
% L is nxn and unit lower triangular and B is nxn and symmetric.
% X solves the systen LXL' + X = B
%
[n,n] = size(B);
X = zeros(n,n);
% Consider the first column in the equation SXT + X = B:
% SX(:,1) + X(:,1) = B(:,1)
% i.e., (S+I)X(:,1) = B(:,1)
Ltilde = L + eye(n,n);
X(:,1) = Ltilde\B(:,1);
% Now assume k>1 and that X(:,1:k-1) is known.
% Let Y = XL' and note that the kth column of Y is given by
% Y(:,k) = XL(k,:)' = X(:,1:k-1)L(k,1:k-1)') + X(:,k)
% Compare the k-th columns in B = LXL' + X = LY + X:
% B(:,k) = LY(:,k) + X(:,k)
% = L*(X(:,1:k-1)L(k,1:k-1)') + X(:,k)) + X(:,k)
% = (L+I)*X(:,k) + L*(X(:,1:k-1)L(k,1:k-1)')
% Solve this for X(:,k):
% (L + I)X(:,k) = B(:,k) - L*(X(:,1:k-1)L(k,1:k-1)')
for k=2:n
v = X(:,1:k-1)*L(k,1:k-1)';
b = B(:,k) - L*v;
X(:,k) = Ltilde\b;
% Recomputing L + I each time through this loop slows things down.
end

function X = SymSolveCVLSym(L,B)
% L is nxn and unit lower triangular and B is nxn and symmetric.
% X solves the systen LXL' + X = B
%
[n,n] = size(B);
X = zeros(n,n);
% Consider the first column in the equation SXT + X = B:
% SX(:,1) + X(:,1) = B(:,1)
% i.e., (S+I)X(:,1) = B(:,1)
Ltilde = L + eye(n,n);
X(:,1) = Ltilde\B(:,1);
% Now assume k>1 and that X(:,1:k-1) is known.
% Let Y = XL' and note that the kth column of Y is given by
% Y(:,k) = XL(k,:)' = X(:,1:k-1)L(k,1:k-1)') + X(:,k)
% Compare the k-th columns in B = LXL' + X = LY + X:
% B(:,k) = LY(:,k) + X(:,k)
% = L*(X(:,1:k-1)L(k,1:k-1)') + X(:,k)) + X(:,k)
% = (L+I)*X(:,k) + L*(X(:,1:k-1)L(k,1:k-1)')
% Solve this for X(:,k):
% (L + I)X(:,k) = B(:,k) - L*(X(:,1:k-1)L(k,1:k-1)')
for k=2:n
v = X(:,1:k-1)*L(k,1:k-1)';
% f we know X(:,1:k-1) then by symmetry we know X(1:k-1,k):
% X(1:k-1,k) = X(k,1:k-1)'
%b = B(:,k);
%for j=1:n
% b(j:n) = b(j:n) - v(j)*L(j:n,j);
%end
X(1:k-1,k) = X(k,1:k-1)';
X(k:n,k) = Ltilde(k:n,k:n)\(B(k:n,k)-L(k:n,:)*v - Ltilde(k:n,1:k-1)*X(k,1:k-1)');

end

% -1 if X(:,1:k-1)*T(1:k-1,k) is computed as a sum of vectors.

% -1 for computing L+I n times, an O(n^3) overhead
% -2 for B(:,k) - L*X(:,1:k-1)*L(k,1:k-1)' instead of
% B(:,k) - L*(X(:,1:k-1)*L(k,1:k-1)')