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