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