% cholinc_pcg_test.m % This MATLAB program demonstrates and evaluates % the use of an incomplete Cholesky decomposition % preconditioner in the conjugate gradient method. % The code demonstrates the use of the method for % a 3-D diffusion matrix on a grid of NxNxN points. % K. Beers. MIT ChE. 10/1/03. v 6/14/05. function iflag_main = cholinc_pcg_test(N,use_mod); iflag_main = 0; % First, set 3-D diffusion matrix Ntot = N^3; N_sq = N^2; % set matrix elements A = spalloc(Ntot,Ntot,7*Ntot); for ix=1:N for iy=1:N for iz=1:N n = get_label_3D(ix,iy,iz,N,N); % master label A(n,n) = 6; % diagonal element if(iz > 1) % U m = get_label_3D(ix,iy,iz-1,N,N); A(n,m) = -1; end if(iz < N) % L m = get_label_3D(ix,iy,iz+1,N,N); A(n,m) = -1; end if(ix > 1) % W m = get_label_3D(ix-1,iy,iz,N,N); A(n,m) = -1; end if(ix < N) % E m = get_label_3D(ix+1,iy,iz,N,N); A(n,m) = -1; end if(iy > 1) % S m = get_label_3D(ix,iy-1,iz,N,N); A(n,m) = -1; end if(iy < N) % N m = get_label_3D(ix,iy+1,iz,N,N); A(n,m) = -1; end end end end disp('Finished setting A matrix'); nnz_A = nnz(A); disp(['nnz(A) = ', int2str(nnz_A)]); % solve system using pcg with no preconditioner h_grid = (N+1)^(-2); b = h_grid*ones(size(A,1),1); % set options for pcg solver tol = 1e-6; maxit = 2*Ntot; [x,iflag,relres,iter_pcg_no_precond] = pcg(A,b,tol,maxit); % Now, use incomplete Cholesky with varying % values of the drop tolerance. drop_tol_logmin = -6; drop_tol_logmax = 0; drop_tol_vect = logspace(drop_tol_logmin,drop_tol_logmax,10); % store values of CPU time relative to use without % preconditioner and numbers of non-zero elements % in incomplete Cholesky factor iter_pcg_cholinc = zeros(size(drop_tol_vect)); nnzR_cholinc = zeros(size(drop_tol_vect)); % set parameters for cholinc if(use_mod) ICHOL_OPTS.michol = 1; end for k=1:length(drop_tol_vect) ICHOL_OPTS.droptol = drop_tol_vect(k); % generate incomplete Cholesky preconditioner R = cholinc(A,ICHOL_OPTS); nnzR_cholinc(k) = nnz(R); [x,iflag,relres,iter_cholinc(k)] = pcg(A,b,tol,maxit,R',R); end % Next, plot performance figure; subplot(2,1,1); semilogx(drop_tol_vect,iter_cholinc); axis([ min(drop_tol_vect), max(drop_tol_vect), ... round(0.9*min(iter_cholinc)), round(1.1*max(iter_cholinc))]); xlabel('drop tolerance in cholinc'); ylabel('# of pcg iterations'); title(['Effect of cholinc preconditioner on pcg, ', ... ' 3-D diffusion matrix, N = ', int2str(N)]); gtext(['# CG iterations with no precond = ', int2str(iter_pcg_no_precond)]); subplot(2,1,2); loglog(drop_tol_vect,nnzR_cholinc,'-.'); % add upper bound and lower bound hold on; loglog(drop_tol_vect,nnz_A*ones(size(drop_tol_vect)),':'); drop_tol_logmid = (drop_tol_logmin+drop_tol_logmax)/2; text(10^drop_tol_logmid,2*nnz_A,'nnz(A)'); nnz_Rup = 2*(N^5); loglog(drop_tol_vect,nnz_Rup*ones(size(drop_tol_vect)),':'); text(10^drop_tol_logmid,2*nnz_Rup,... 'upper bound on nnz(R), R=chol(A)'); axis([ min(drop_tol_vect), max(drop_tol_vect), ... 0.1*nnz_A, 10*nnz_Rup]); xlabel('drop tolerance in cholinc'); ylabel('nnz(R), R from cholinc'); % if use modified incomplete Cholesky, add note if(use_mod) gtext('modified cholinc'); end % make spy plots of A and of R for various values of the % drop tolerance figure; % spy plot of A subplot(2,2,1); spy(A); % spy plot of R with droptol = 1e-4; ICHOL_OPTS.droptol = 1e-4; R = cholinc(A,ICHOL_OPTS); subplot(2,2,2); spy(R); gtext('droptol = 10^{-4}'); % spy plot of R with droptol = 1e-3; ICHOL_OPTS.droptol = 1e-3; R = cholinc(A,ICHOL_OPTS); subplot(2,2,3); spy(R); gtext('droptol = 10^{-3}'); % spy plot of R with droptol = 1e-2; ICHOL_OPTS.droptol = 1e-2; R = cholinc(A,ICHOL_OPTS); subplot(2,2,4); spy(R); gtext('droptol = 10^{-2}'); iflag_main = 1; return; % ---------------------------------------- % This routine computes the label of a grid % point in the 3-D system. function label = get_label_3D(ix,iy,iz,Nx,Ny); label = (iz-1)*Nx*Ny + (ix-1)*Ny + iy; return;