%% Iterative Solution to Linear System
% Solve a square linear system using cgs with default settings, and then adjust the tolerance and number of iterations used in the solution process.
%
% Create a random sparse matrix A with 50% density. Also create a random vector for the right-hand side of \textrm{Ax}=b.

rng default
A = sprand(600,600,.5);
A = A'*A + speye(size(A));
b = rand(600,1);

%
% Solve \textrm{Ax}=b using cgs. The output display includes the value of the relative residual error \frac{\left\|\textrm{Ax}-b\right\|}{\left\|b\right\|}.

x = cgs(A,b);

% By default cgs uses 20 iterations and a tolerance of 1e-6, and the algorithm is unable to converge in those 20 iterations for this matrix. Since the residual is still large, it is a good indicator that more iterations (or a preconditioner matrix) are needed. You also can reduce the tolerance to make it easier for the algorithm to converge.
%
% Solve the system again using a tolerance of 1e-4 and 40 iterations.

x = cgs(A,b,1e-4,40);

% Even with a looser tolerance and more iterations cgs does not converge. When an iterative algorithm stalls in this manner it is a good indication that a preconditioner matrix is needed.
%
% Calculate the incomplete Cholesky factorization of A, and use the L factor as a preconditioner input to cgs.

L = ichol(A);
x = cgs(A,b,1e-4,40,L);

% Using a preconditioner improves the numerical properties of the problem enough that cgs is able to converge.
%
% Copyright 2019 The MathWorks, Inc.
%
%