function mit18086_ilu(n) %MIT18086_ILU % Sets up a 2d Poisson problem, and computes the LU % as well as ILU decomposition for different tolerances % (including ILU 0 and tol=0, which is normal LU). % Run times are clocked and errors are measured. % 03/2008 by Benjamin Seibold % Feel free to modify for teaching and learning. %------------------------------------------------------------------- tols = {'0' 1e-1 1e-3 1e-5 0}; % set of ILU tolerances n_runs = 5; % number of runs for cpu-time %------------------------------------------------------------------- if nargin<1, n = 16; end A = delsq(numgrid('S',n+2)); % 2d Poisson matrix for square t = 1:n_runs; for k = 1:n_runs, tic, [L,U] = lu(A); t(k) = toc; end t_lu = median(t); disp('========================================') disp(sprintf('cond(A) = %0.2f',condest(A))) disp(sprintf('runtime LU = %0.3f ms',t_lu*1e3)) disp(sprintf('norm(A-L*U) = %0.2e',normest(A-L*U))) disp('========================================') clf subplot(length(tols)+1,3,1), spy(A), title('2d Poisson matrix A') subplot(length(tols)+1,3,2), spy(L,'b') hold on, spy(U,'r'), hold off, title('LU factorization') drawnow for itol = 1:length(tols) for k = 1:n_runs tic, [Lapp,Uapp] = luinc(A,tols{itol}); t(k) = toc; end t_ilu = median(t); Aapp = Lapp*Uapp; for k = 1:n_runs, tic, Iapp = Aapp\A; t(k) = toc; end t_iapp = median(t); if tols{itol}=='0', name = 'ILU 0'; else name = sprintf('%0.1e',tols{itol}); end disp(['tolerance = ',name]) disp('----------------------------------------') disp(sprintf('cond(inv(Aapp)*A) = %0.2f',condest(Iapp))) disp(sprintf('runtime ILU = %0.3f ms',t_ilu*1e3)) disp(sprintf('runtime inv(Aapp)*A = %0.3f ms',t_iapp*1e3)) disp(sprintf('norm(A-Aapp) = %0.2e',normest(A-Aapp))) disp(sprintf('norm(A-Aapp.*spones(A)) = %0.2e',... normest(A-Aapp.*spones(A)))) disp('========================================') subplot(length(tols)+1,3,itol*3+1), spy(Lapp,'b'), hold off text(n^2*.9,n^2*.2,name,'HorizontalAlignment','right') title('L_{app}') subplot(length(tols)+1,3,itol*3+2), spy(Uapp,'r'), hold off title('U_{app}') subplot(length(tols)+1,3,itol*3+3), spy(Aapp,'b'), hold off title('L_{app}*U_{app}') drawnow end