%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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);