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