Biconjugate gradients stabilized (l) method
x = bicgstabl(A,b)
x = bicgstabl(afun,b)
x = bicgstabl(A,b,tol)
x = bicgstabl(A,b,tol,maxit)
x = bicgstabl(A,b,tol,maxit,M)
x
= bicgstabl(A,b,tol,maxit,M1,M2)
x = bicgstabl(A,b,tol,maxit,M1,M2,x0)
[x,flag] = bicgstabl(A,b,...)
[x,flag,relres] = bicgstabl(A,b,...)
[x,flag,relres,iter] = bicgstabl(A,b,...)
[x,flag,relres,iter,resvec] = bicgstabl(A,b,...)
x = bicgstabl(A,b)
attempts to solve the system of linear equations
A*x=b
for x
. The
n
byn
coefficient matrix A
must be
square and the righthand side column vector b
must have length
n
.
x = bicgstabl(afun,b)
accepts a function handle afun
instead of the matrix A
. afun(x)
accepts a vector input
x
and returns the matrixvector product A*x
. In all of the
following syntaxes, you can replace A
by afun
.
x = bicgstabl(A,b,tol)
specifies the tolerance of the method. If
tol
is [] then bicgstabl
uses the default, 1e6.
x = bicgstabl(A,b,tol,maxit)
specifies the maximum number of
iterations. If maxit
is [] then bicgstabl
uses the
default, min(N,20)
.
x = bicgstabl(A,b,tol,maxit,M)
and
x
= bicgstabl(A,b,tol,maxit,M1,M2)
use preconditioner M
or M=M1*M2
and effectively solve the system A*inv(M)*x = b
for x. If M
is [] then a preconditioner is not applied. M
may be a function handle returning M\x
.
x = bicgstabl(A,b,tol,maxit,M1,M2,x0)
specifies the initial guess. If
x0
is [] then bicgstabl
uses the default, an all zero
vector.
[x,flag] = bicgstabl(A,b,...)
also returns a
convergence flag
:
Flag  Convergence 

 bicgstabl converged to the desired tolerance
tol within maxit iterations. 
 bicgstabl iterated maxit times but did not
converge. 
 Preconditioner M was illconditioned. 


 One of the scalar quantities calculated during 
[x,flag,relres] = bicgstabl(A,b,...)
also returns
the relative residual norm(bA*x)/norm(b)
. If flag
is
0
, relres <= tol
.
[x,flag,relres,iter] = bicgstabl(A,b,...)
also
returns the iteration number at which x
was computed, where 0 <=
iter <= maxit
. iter
can be k/4
where
k
is some integer, indicating convergence at a given quarter
iteration.
[x,flag,relres,iter,resvec] = bicgstabl(A,b,...)
also returns a vector of the residual norms at each quarter iteration, including
norm(bA*x0)
.
You can pass inputs directly to bicgstabl
:
n = 21; A = gallery('wilk',n); b = sum(A,2); tol = 1e12; maxit = 15; M = diag([10:1:1 1 1:10]); x = bicgstabl(A,b,tol,maxit,M);
You can also use a matrixvector product function:
function y = afun(x,n) y = [0; x(1:n1)] + [((n1)/2:1:0)'; (1:(n1)/2)'].*x+[x(2:n); 0];
and a preconditioner backsolve function:
function y = mfun(r,n) y = r ./ [((n1)/2:1:1)'; 1; (1:(n1)/2)'];
as inputs to bicgstabl
:
x1 = bicgstabl(@(x)afun(x,n),b,tol,maxit,@(x)mfun(x,n));
This example demonstrates the use of a preconditioner.
Load west0479
, a real 479by479 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 = 1e12; maxit = 20;
Use bicgstabl
to find a solution at the requested tolerance and number of iterations.
[x0,fl0,rr0,it0,rv0] = bicgstabl(A,b,tol,maxit);
fl0
is 1 because bicgstabl
does not converge to the requested tolerance 1e12
within the requested 20 iterations. In fact, the behavior of bicgstabl
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 bicgstabl
.
semilogy(0:0.25: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',1e5));
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',1e6)); [x1,fl1,rr1,it1,rv1] = bicgstabl(A,b,tol,maxit,L,U);
fl1
is 0 because bicgstabl
drives the relative residual to 1.0257e015
(the value of rr1
). The relative residual is less than the prescribed tolerance of 1e12
at the sixth iteration (the value of it1
) when preconditioned by the incomplete LU factorization with a drop tolerance of 1e6
. The output rv1(1)
is norm(b)
, and the output rv1(9)
is norm(bA*x2)
since bicgstabl
uses quarter iterations.
You can follow the progress of bicgstabl
by plotting the relative residuals at each iteration starting from the initial estimate (iterate number 0).
semilogy(0:0.25:it1,rv1/norm(b),'o'); h = gca; h.XTick = 0:0.25:it1; xlabel('Iteration number'); ylabel('Relative residual');