% demo_fullrank.m
%
% Weighted LS with full rank C and D.

fprintf('Index\tfl\tit\tres\n')

%idx_list = [4, 6, 8, 12, 16, 24, 32, 48, 64]';
idx_list = 48;
fl_list = [];
fl0_list = [];

for k = 1 : length(idx_list)

   idx = idx_list(k);

   mC = ceil(1.5 * idx);
   nC = idx;
   mD = mC;
   nD = nC;

   m = mC * mD;
   n = nC * nD;

   C = 5 - 10 * rand(mC, nC);
   D = 5 - 10 * rand(mD, nD);
   w = 1 + 9 * rand(m, 1);
   b = 1 - 2 * rand(m, 1);

   % Compute exact answer.
   flops(0);
   M = kron(C,D);
   for i = 1 : m
     M(i,:) = M(i,:) / sqrt(w(i));
   end
   [Q, T, S, r] = cod(M, 0, n);
   b2 = b ./ sqrt(w);
   b3 = Q' * b2;
   b4 = [ T \ b3(1:r); zeros(n-r,1) ];
   x0 = S * b4;
   fl0 = flops;
   res0 = norm(b2 - M * x0) / norm(b);
   clear M b2 b3 b4

   % Use CG to solve the system.
   flops(0);
   [x, it] = wls(C, D, w, b);
   fl = flops;
   res = norm((b - KPmult(C, D, x)) ./ sqrt(w)) / norm(b);
   err = norm(x0 - x) / norm(b);

   % Results.
   fprintf('%2d & %9d & %2d & %5.2f & %12d & %5.2f & %0.2e \\\\\n', idx, fl, it, res, fl0, res0, err)
   fl_list = [fl_list; fl];
   fl0_list = [fl0_list; fl0];

end % for k
