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

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

% Copyright 2015 The MathWorks, Inc.


%%
% Load |A = 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 |qmr| to find a solution at the requested tolerance and number of
% iterations.
[x0,fl0,rr0,it0,rv0] = qmr(A,b,tol,maxit);

%%
% |fl0| is 1 because |qmr| does
% not converge to the requested tolerance |1e-12| within
% the requested 20 iterations. The seventeenth iterate is the best approximate
% solution and is the one returned as indicated by |it0 = 17|. MATLAB stores the residual
% history in |rv0|.

%%
% Plot the behavior of |qmr|.
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 the preconditioner with |ilu|, since the matrix |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] = qmr(A,b,tol,maxit,L,U);

%%
% |fl1| is 0 because |qmr| drives
% the relative residual to |4.1410e-014| (the value
% of |rr1|). The relative residual is less than the
% prescribed tolerance of |1e-12| at the sixth 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(7)| is |norm(b-A*x2)|.

%%
% You can follow the progress of |qmr| 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');