  function [y,z] = SolveSym(A,C,f,g)
% Solves [A C;C' 0][y;z] = [f;g]
% Compute L = [G 0 ; H K]..
G = chol(A,'lower');
H = (G\C)';
K = chol(H*H','lower');
% Solve L[u;v] = [f;g]
u = G\f;
v = K\(g - H*u);
% Solve U[y;z] = [u;v]
z = -K'\v;
y = G'\(u -H'*z);

[n,r] = size(C);
condBlockMatrix = cond([A C;C' zeros(r,r)])
condA = cond(A)
normH = norm(H*H')
% The accuracy of the backslash method squares with what we know about
% Gaussian elimination with pivoting: 
%               relerr approx = cond(matrix)*unit roundoff = 10^7 *10^-16
% For the SolveSym method, because A is ill conditioned H is huge: O(10^8)
% Indeed, H*H' will have O(1) errors since 10^16*10^-16 = 1
% This means that we will have KK' - HH' = O(1) and we see the same level
% of error in the SolveSym solution.

% 1 for noticing conditioning of A could be a problem
% 1 for observing that it causes H to be huge and for computation of K to
%   have large error

% -1 for constructing L and U to find y and z
% Up to -2 for computing [y;z] or intermediate quantities incorrectly

