function [x, it] = solve(B, C, b, m, n, tol, maxit, guess)
% [x, it] = solve(B, C, b, m, n, tol, maxit, guess)
%
% Solves the mn^2-by-mn^2 linear system
%   A x = b
% where
%   A_ij = B_ij \kron C_ij
% for i,j = 1 : m. The matrix A is SPD, and the system is solved
% using Preconditioned Conjugate Gradients.
%
% The built-in Matlab function pcg.m is used to implement the PCG
% algorithm. The functions afun.m, m1fun.m, and m2fun.m must be
% available.
%
% Input:
%   B, C  : mn-by-mn SPD matrices
%   b     : mn^2 vector
%   m, n  : integers >= 1
% Optional inputs:
%   tol   : convergence tolerance for PCG. Default is 1e-10.
%   maxit : maximum interations for PCG. Default is 10 m n^2.
%   guess : mn^2 vector. Guess for solution to linear system. Default is 0.
% Output:
%   x     : mn^2 vector. Approximate solution to linear system.
%   it    : integer >= 0. Number of iterations in PCG.

% Check number of input arguments.
if nargin < 5
   error('solve error: at least 5 arguments required')
end

% Check sizes.
mn = m * n;
mn2 = mn * n;

[mB, nB] = size(B);
if (mB ~= mn) | (nB ~= mn)
   error('solve error: B must be mn-by-mn')
end

[mC, nC] = size(B);
if (mC ~= mn) | (nC ~= mn)
   error('solve error: C must be mn-by-mn')
end

[mb, nb] = size(b);
if (mb ~= mn2) | (nb ~= 1)
   error('solve error: b must be length mn^2 column vector')
end

% Compensate for missing arguments.

if nargin < 6
   tol = 1e-10;
end

if nargin < 7
   maxit = 10 * m * n * n;
end

if nargin < 8
   guess = [];
end

% Allocate G and H, both mn-by-n.
G = zeros(mn, n);
H = zeros(mn, n);

% Cholesky factor the diagonal blocks.
for i = 1 : m
   ii = (i-1)*n + (1:n);
   G(ii, :) = chol(B(ii,ii))';
   H(ii, :) = chol(C(ii,ii))';
end

% Solve using Preconditioned Conjugate Gradients.
[x, junk1, junk2, it_pcg] = pcg('afun', b, tol, maxit, ...
   'm1fun', 'm2fun', guess, B, C, G, H, m, n);

if nargout > 1
   it = it_pcg;
end
