% %============================================================================%
% % Duke University                                                            %
% % K. P. Trofatter                                                            %
% % kpt2@duke.edu                                                              %
% %============================================================================%

%% Environment ================================================================%
clear();
clc();

%% Parameters =================================================================%
% linear system
m = 10000;
n = 10000;
vr = 1.0;

% zlsq
gamma = 1.0;

% iterations
tol = 1.0e-3; % residual norm
imax = 100; % max iterations

% gmres
rst = [];

%% Run ========================================================================%
% roll linear system
A = vr / sqrt(2) * (randn(m, n) + 1.0j * randn(m, n));
b = randn(m, 1);
x0 = zeros(n, 1);
%x0 = randn(n, 1);

% build zlsq equation
[A_lsq, b_lsq] = zlsq(@(x) A * x, @(x) A' * x, b, gamma);
bn = norm(b_lsq);

% jacobi preconditioner
P = diag(A' * A) + gamma;
P = 1.0 ./ P;
Pif = @(x) P .* x;

% solve zlsq equation
% zgd_quad
fprintf('zgd_quad(): ');
tic();
[x_zgd, rn_zgd] = zgd_quad(A_lsq, b_lsq, x0, tol, imax);
toc();
% zcg
fprintf('zcg(): ');
tic();
[x_zcg, rn_zcg] = zcg(A_lsq, b_lsq, x0, [], tol, imax);
toc();
% pcg (matlab's version of zcg)
fprintf('pcg(): ');
tic();
[x_pcg, ~, ~, ~, rn_pcg] = pcg(A_lsq, b_lsq, tol / bn, imax, [], [], x0);
toc();
% zgmres
fprintf('zgmres(): ');
tic();
[x_zgm, rn_zgm] = zgmres(A_lsq, b_lsq, x0, tol, imax);
toc();
% gmres (matlab's version of zgmres)
fprintf('gmres(): ');
tic();
[x_gmr, ~, ~, ~, rn_gmr] = gmres(A_lsq, b_lsq, rst, tol / bn, imax, [], [], x0);
toc();

%% Draw =======================================================================%
% figure
fh = figure(1);
clf(fh);
% axis
ah = axes('Parent', fh);
% draw
hold(ah, 'on');
gh(1) = semilogy(rn_zgd, 'Color', 'b');
gh(2) = semilogy(rn_pcg, 'Color', 'k');
gh(3) = semilogy(rn_zcg, 'Color', 'r');
gh(4) = semilogy(rn_gmr, 'Color', 'k');
gh(5) = semilogy(rn_zgm, 'Color', 'g');
hold(ah, 'off');
% annotate
grid(ah, 'on');
legend(ah, gh, {'zgd\_quad()', 'pcg()', 'zcg()', 'gmres()', 'zgmres()'});
xlabel(ah, 'iteration');
ylabel(ah, 'residual norm');


%==============================================================================%
%                                                                              %
%                                                                              %
%                                                                              %
%==============================================================================%