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.