www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/UsingIcholAsAPreconditionerExample.m
%% Using ichol as a Preconditioner % This example shows how to use an incomplete Cholesky factorization as a % preconditioner to improve convergence. %% % Create a symmetric positive definite matrix, |A|. % Copyright 2015 The MathWorks, Inc. N = 100; A = delsq(numgrid('S',N)); %% % Create an incomplete Cholesky factorization as a preconditioner for % |pcg|. Use a constant vector as the right hand side. As a baseline, % execute |pcg| without a preconditioner. b = ones(size(A,1),1); tol = 1e-6; maxit = 100; [x0,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit); %% % Note that |fl0 = 1| indicating that |pcg| did not drive the relative residual % to the requested tolerance in the maximum allowed iterations. Try the % no-fill incomplete Cholesky factorization as a preconditioner. L1 = ichol(A); [x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L1,L1'); %% % |fl1 = 0|, indicating that |pcg| converged to the requested tolerance and % did so in 59 iterations (the value of |it1|). Since this matrix is a % discretized Laplacian, however, using modified incomplete Cholesky can % create a better preconditioner. A modified incomplete Cholesky % factorization constructs an approximate factorization that preserves the % action of the operator on the constant vector. That is, % |norm(A*e-L*(L'*e))| will be approximately zero for |e = % ones(size(A,2),1)| even though |norm(A-L*L','fro')/norm(A,'fro')| is not % close to zero. It is not necessary to specify type for this syntax since % |nofill| is the default, but it is good practice. opts.type = 'nofill'; opts.michol = 'on'; L2 = ichol(A,opts); e = ones(size(A,2),1); norm(A*e-L2*(L2'*e)) [x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L2,L2'); %% % |pcg| converges (|fl2 = 0|) but in only 38 iterations. Plotting all three % convergence histories shows the convergence. semilogy(0:maxit,rv0./norm(b),'b.'); hold on semilogy(0:it1,rv1./norm(b),'r.'); semilogy(0:it2,rv2./norm(b),'k.'); legend('No Preconditioner','IC(0)','MIC(0)'); %% % The plot shows that the modified incomplete Cholesky preconditioner % creates a much faster convergence. % % You can also try incomplete Cholesky factorizations with threshold % dropping. The following plot shows convergence of |pcg| with % preconditioners constructed with various drop tolerances. L3 = ichol(A, struct('type','ict','droptol',1e-1)); [x3,fl3,rr3,it3,rv3] = pcg(A,b,tol,maxit,L3,L3'); L4 = ichol(A, struct('type','ict','droptol',1e-2)); [x4,fl4,rr4,it4,rv4] = pcg(A,b,tol,maxit,L4,L4'); L5 = ichol(A, struct('type','ict','droptol',1e-3)); [x5,fl5,rr5,it5,rv5] = pcg(A,b,tol,maxit,L5,L5'); figure semilogy(0:maxit,rv0./norm(b),'b-','linewidth',2); hold on semilogy(0:it3,rv3./norm(b),'b-.','linewidth',2); semilogy(0:it4,rv4./norm(b),'b--','linewidth',2); semilogy(0:it5,rv5./norm(b),'b:','linewidth',2); legend('No Preconditioner','ICT(1e-1)','ICT(1e-2)', ... 'ICT(1e-3)','Location','SouthEast'); %% % Note the incomplete Cholesky preconditioner constructed with drop % tolerance |1e-2| is denoted as |ICT(1e-2)|. % % As with the zero-fill incomplete Cholesky, the threshold dropping % factorization can benefit from modification (i.e. |opts.michol = 'on'|) % since the matrix arises from an elliptic partial differential equation. % As with |MIC(0)|, the modified threshold based dropping incomplete % Cholesky will preserve the action of the preconditioner on constant % vectors, that is |norm(A*e-L*(L'*e))| will be approximately zero.