%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% problem4.m
% Solve Ax = b using GMRES with different preconditioners.  
% 
% Sungwoo Park
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;
close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate the problem A*x = b
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load west0479;
A = west0479;
n = size(A,1);
x_true = ones(n,1);
b = A*x_true;

restart = 100;
tol = 1e-2;
fprintf('#of restarting = %d\ttolerance = %.0e\n\n',restart,tol);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (0) GMRES without preconditioner
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic;
[x,flag,relres,iter]= gmres(A,b,restart,tol);
t0 = toc;
fprintf('[0] GMRES without preconditioner\n');
fprintf('(time for GMRES)\t(outer iter)\t(inner iter)\t(relative residual)\n');
fprintf('%.4f(sec)\t\t\t%d\t\t\t\t%d\t\t\t\t%.3e\n\n',t0,iter(1),iter(2),relres);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (1)luinc changing the drop tolerance.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('[1] GMRES with luinc changing the drop tolerance.\n');
fprintf('(drop_tol)\t(t_luinc)\t(t_GMRES)\t(t_total)\t(itr_out)\t(itr_in)\t(rel res)\t(nz of L)\t(nz of U)\n');
drop_tol = 1e-4;

for i=1:10,
    drop_tol = drop_tol*1e-2;
    tic;
    [L1,U1] = luinc(A,drop_tol);
    t1 =toc;
    tic;
    [x,flag,relres,iter] = gmres(A,b,restart,tol,[],L1,U1);
    t2 = toc;
    fprintf('%.0e\t\t%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\t%d\t\t%d\n'...
        ,drop_tol,t1,t2,t1+t2,iter(1),iter(2),relres,nnz(L1),nnz(U1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (2)luinc using the modified version to preserve row sums.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n[2] GMRES with luinc using the modified version to preserve row sums.\n');
fprintf('(drop_tol)\t(t_luinc)\t(t_GMRES)\t(t_total)\t(itr_out)\t(itr_in)\t(rel res)\t(nz of L)\t(nz of U)\n');
options.milu = 0;
options.droptol = 1e-6;

tic;
[L1,U1] = luinc(A,options);
t1 =toc;
tic;
[x,flag,relres,iter] = gmres(A,b,restart,tol,[],L1,U1);
t2 = toc;
fprintf('%.0e\t\t%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\t%d\t\t%d\n'...
    ,options.droptol,t1,t2,t1+t2,iter(1),iter(2),relres,nnz(L1),nnz(U1));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (3)luinc using a matrix reordering before factorization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n[3] GMRES with luinc changing the drop tolerance.\n');
fprintf('(drop_tol)\t(t_luinc)\t(t_GMRES)\t(t_total)\t(itr_out)\t(itr_in)\t(rel res)\t(nz of L)\t(nz of U)\n');
drop_tol = 1e-6;

% reordering
p = colamd(A);
A_reorder = A(:,p);

tic;
[L1,U1] = luinc(A_reorder,drop_tol);
t1 =toc;
tic;
[x,flag,relres,iter] = gmres(A_reorder,b,restart,tol,[],L1,U1);
t2 = toc;
fprintf('%.0e\t\t%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\t%d\t\t%d\n'...
    ,drop_tol,t1,t2,t1+t2,iter(1),iter(2),relres,nnz(L1),nnz(U1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (4)luinc using left, right, and two-sided preconditioning.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n[4] GMRES with luinc using LEFT, RIGHT, and TWO-sided.\n');
drop_tol = 1e-6;
tic;
[L,U] = luinc(A,drop_tol);
t1 =toc;
fprintf('(type)\t(drop_tol)\t(t_luinc)\t(t_GMRES)\t(t_total)\t(itr_out)\t(itr_in)\t(rel res)\t(nz of L)\t(nz of U)\n');


% a)left preconditioning
tic;
[x,flag,relres,iter] = run_gmres(A,b,restart,tol,L,U,[],'left');
t2 = toc;
fprintf('LEFT  \t%.0e\t\t%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\t%d\t\t%d\n'...
    ,drop_tol,t1,t2,t1+t2,iter(1),iter(2),relres,nnz(L),nnz(U));

% b)right preconditioning
tic;
[x,flag,relres,iter] = run_gmres(A,b,restart,tol,L,U,[],'right');
t2 = toc;
fprintf('RIGHT\t%.0e\t\t%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\t%d\t\t%d\n'...
    ,drop_tol,t1,t2,t1+t2,iter(1),iter(2),relres,nnz(L),nnz(U));

% c)two-sided preconditioning
tic;
[x,flag,relres,iter] = run_gmres(A,b,restart,tol,L,U,[],'two');
t2 = toc;
fprintf('TWO   \t%.0e\t\t%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\t%d\t\t%d\n'...
    ,drop_tol,t1,t2,t1+t2,iter(1),iter(2),relres,nnz(L),nnz(U));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (5) approximate inverse preconditioning.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n[5] GMRES with Approximate Inverse Preconditioning.\n');
W = [];
tic;
for i = 1:n
    e = spalloc(n,1,1);
    e(i)=1;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Check the sparsity pattern of A's i-th row
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    col_nz = find(A(i,:));
    nz = nnz(A(i,:));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Construct A_nz cosisting of all columns of A corresponding to nonzero
    % entries of A's i-th row. And solve LS problem of min||A_nz*w-e||
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    A_nz = A(:,col_nz);
    [C,R] = qr(A_nz,e);
    w = spalloc(n,1,nz);
    w_ls = R\C;
    w(col_nz) = w_ls;
    W = [W w];
end
t1 = toc;
fprintf('non-zero element of W = %d\n',nnz(W));
fprintf('non-zero element of A = %d\n',nnz(A));
fprintf('||I-AW||=%f\n',norm(A*W-eye(n,n),'fro'));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gmres with right preconditioning using Approximate Inverse W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tic;
[x,flag,relres,iter] = run_gmres(A,b,restart,tol,[],[],W,'right_aip');
t2 = toc;
fprintf('(t_AIP)  \t(t_GMRES)\t(t_total)\t(itr_out)\t(itr_in)\t(rel res)\n');
fprintf('%.4f(s)\t%.4f(s)\t%.4f(s)\t%d\t\t\t%d\t\t\t%.3e\n\n'...
     ,t1,t2,t1+t2,iter(1),iter(2),relres);