www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/UsingCgsWithAPreconditionerExample.m

    %% Using cgs with a Preconditioner.
% This example demonstrates the use of a preconditioner.

% Copyright 2015 The MathWorks, Inc.


%%
% Load |west0479|, a real 479-by-479 nonsymmetric sparse matrix.
load west0479;
A = west0479;

%%
% Define |b| so that the true solution is a vector of all ones.
b = full(sum(A,2));

%%
% Set the tolerance and maximum number of iterations.
tol = 1e-12;
maxit = 20;

%%
% Use |cgs| to find a solution at the requested tolerance and number of
% iterations.
[x0,fl0,rr0,it0,rv0] = cgs(A,b,tol,maxit);

%%
% |fl0| is 1 because |cgs| does not converge to the requested tolerance
% |1e-12| within the requested 20 iterations. In fact, the behavior of
% |cgs| is so poor that the initial guess (|x0 = zeros(size(A,2),1)| is the
% best solution and is returned as indicated by |it0 = 0|. MATLAB stores
% the residual history in |rv0|.

%%
% Plot the behavior of |cgs|.
semilogy(0:maxit,rv0/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');

%%
% The plot shows that the solution does not converge. You can use a preconditioner to improve the outcome.

%%
% Create a preconditioner with |ilu|, since |A| is nonsymmetric.
%   
%   [L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));
%
%  Error using ilu
%  There is a pivot equal to zero.  Consider decreasing the 
%  drop tolerance or consider using the 'udiag' option.

%%
% MATLAB cannot construct the incomplete LU as it would result in a singular factor, which is useless as a preconditioner.

%%
% You can try again with a reduced drop tolerance, as indicated by the
% error message.
[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = cgs(A,b,tol,maxit,L,U);

%%
% |fl1| is 0 because |cgs| drives
% the relative residual to |4.3851e-014| (the value
% of |rr1|). The relative residual is less than the
% prescribed tolerance of |1e-12| at the third iteration
% (the value of |it1|) when preconditioned by the incomplete
% LU factorization with a drop tolerance of |1e-6|.
% The output |rv1(1)| is |norm(b)| and
% the output |rv1(14)| is |norm(b-A*x2)|.

%%
% You can follow the progress of |cgs| by plotting
% the relative residuals at each iteration starting from the initial
% estimate (iterate number 0).
semilogy(0:it1,rv1/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');