%% Iterative Solution to Linear System
% Solve a square linear system using bicgstabl 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 b for the right-hand side of \textrm{Ax}=b.

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

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

x = bicgstabl(A,b);

% By default bicgstabl 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 use a larger tolerance to make it easier for the algorithm to converge.
%
% Solve the system again using a tolerance of 1e-4 and 100 iterations.

x = bicgstabl(A,b,1e-4,100);

% Even with a looser tolerance and more iterations, the residual error does not improve much. 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 bicgstabl.

L = ichol(A);
x = bicgstabl(A,b,1e-4,100,L');

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