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');