  function [x1,x2,its] = SparseSolve(A1,A2,C,b1,b2,tol,itMax)
% Solves the system
%        A1   C   x1      b1
%               *      =
%        C'  A2   x2      b2
%
% where 
%       A1 is n1xn1, symmetric positive definite in sparse format
%       A2 is n2xn2, symmetric positive definite in sparse format
%       C  is n1xn2  in sparse formate
%       b1 is n1x1
%       b2 is n2x1
% Uses block Gauss-Seidel with a maximum of itMax iterations.
% Termination criteria: norm(Ax-b)<=tol*norm(b)
% its houses the number of required iterations.

% Precompute the sparse Cholesky factors of A1 and A2...
[G1,p1,S1] = chol(A1,'lower'); spy(S1), spy(G1)
[G2,p2,S2] = chol(A2,'lower'); spy(S2), spy(G2)
% For an initial guess, set x1 = 0 and x2 = 0...
[n1,n2] = size(C);
x1 = zeros(n1,1);
x2 = zeros(n2,1);

% Other initializations...
k = 1;
beta = tol*norm([b1;b2]);
% The Block Gauss-Seidel iteration...
while k<itMax && norm([A1*x1+C*x2-b1;C'*x1+A2*x2-b2])>=beta
    x1 = S1*(G1'\(G1\(S1'*(b1-C*x2))));
    x2 = S2*(G2'\(G2\(S2'*(b2-C'*x1))));
    k = k+1;
end
its = k-1;

% -3 for Cholesky factorization inside loop
% -1 for computing norm(b) repeatedly
% up to -2 for not taking advantage of triangular solves or S^{-1} = S'
% -1 for incorrect stopping criterion
% -1 for not using new x1

% Some students cleverly observed that norm(A*[x1_new; x2_new] - b) equals
%      norm(C*(x2_new - x2_old))!
