% demo.m
%
% Demonstrates the use of solve.m. 

close all

m = 2;
n_list = [2, 4, 8, 16, 32]'; % , 64, 128, 256, 512]';
fl_list = [];
er_list = [];

for i = 1 : length(n_list)
   % Define problem size.
   n = n_list(i);
   fprintf('n = %d\n', n)
   
   % Define a matrix A, A_ij = B_ij \kron C_ij
   fprintf('Generating B\n')
   B = get_spd(m*n);
   fprintf('Generating C\n')
   C = get_spd(m*n);
   % The matrix A is
   %   A = [ B11 \kron C11   (B21 \kron C21)' ]
   %       [ B21 \kron C21   B22 \kron C22    ]
   
   % Correct answer.
   x0 = 1 - 2 * rand(m*n*n, 1);
   
   % Compute RHS.
   b = afun(x0, B, C, [], [], m, n);
   
   % Solve using solve.m.
   flops(0);
   [x, it] = solve(B, C, b, m, n);
   fl = flops;
   err = norm(x - x0) / norm(x0);
   
   fprintf('\t%d\t%e\n', fl, err)
   fl_list = [fl_list; fl];
   er_list = [er_list; err];
end

figure(1)
loglog(n_list, fl_list, n_list, fl_list, 'bo')
title('Flops vs. n')

figure(2)
loglog(n_list, er_list, 'bo')
title('Error vs. n')
