From this page you can get a Matlab® implementation of the BiCGstab(ell)
The BiCGstab(ell) algorithm can be used for computing the solution X of a linear system A*X=B, where A is a square n by n matrix and B is an n by p array (with p small).The matrix can be real or complex, Hermitian or non-Hermitian, .... The algorithm is effective especially in case A is sparse and of large size.
The BiCGstab(ell) algorithm has been introduced in
|Requires:||Matlab Version 5.1.|
|Function:||[X,HISTORY] = CGSTAB(A,B,X0,TR0,OPTIONS,L,U);|
|X = CGSTAB(A,B) attempts to solve
the system of linear equations A*X = B for X.
The coefficient matrix A must be square and the
right hand side. B must by an N-by-P,
where A is N-by-N.
|CGSTAB will start iterating from
an initial guess which by default is an array of size N-by-P
all zeros. Iterates are produced until the method either converges, fails,
or has computed the maximum number of iterations. Convergence is achieved
when an iterate X has relative residual NORM(B(:,I)-A*X(:,I))/NORM(B(:,I))
less than or equal to the tolerance of the method for all I=1:P.
The default tolerance is 1e-8. The default maximum
number of iterations is 300. No preconditioning is used. CGSTAB
with ell= 4 as default.
|[X,HIST] = CGSTAB(A,B) returns
also the convergence history (the norms of the subsequential residuals).
|... = CGSTAB('Afun',B).
The first input argument is either a square matrix (which can be full or sparse, symmetric or non symmetric, real or complex), or a string containing the name of an M-file which applies a linear operator to a given array of size SIZE(B). In the latter case, N = SIZE(B,1).
|The remaining input arguments are optional and can be given
in practically any order:
... = CGSTAB(...,X0,TR0,OPTIONS,M)
|X0||An N-by-P array of doubles, an initial guess.|
|TR0||An N-by-1 array of doubles, the shadow residual.|
|OPTIONS||A structure containing additional parameters.|
|M||String(s) or array(s) defining the preconditioner.|
|If X0 is not specified then X0
If TR0 is not specified then TR0 = (B-A*X0)*RAND(P,1).
The OPTIONS structure specifies certain parameters in the algorithm.
|OPTIONS.Tol||Tol is the convergence tolerance.||1e-8|
|OPTIONS.MaxIt||Maximum number of iterations.||300|
|OPTIONS.Disp||Shows size of intermediate residuals ('yes'/'no').||'no'|
|OPTIONS.ell||ell for BiCGstab(ell).||4|
|specifies in what sense the systems of equations are to be solved
0: each equation with a relative accuracy Tol:
NORM(A*X(:,I)-B(:,I)) < Tol*NORM(B(:,I)) for all I=1:P;
1: each equation with an accuracy of Tol:
NORM(A*X(:,I)-B(:,I)) < Tol for all I=1:P;
2: each equation in the space with a relative accuracy Tol:
NORM(A*X*MU-B*MU) < Tol*NORM(B*MU) for all MU of size P-by-1.
|De equations can be solved simultaneously (using a banded version of the BiCGstab(ell) algorithm) ('yes' ) or sequentially (applying BiCGstab(ell ) on to the column vectors of the right-hand side B ) ('no').||'yes'|
|Specifies in what way the preconditioner is to be applied. The preconditioning will always be explicit, but it can be left ('left' ), right ('right'), central (=two-sided if two matrices have been specified) ('central',), or via the Eisenstat trick (if A is matrix and the preconditioner is a diagonal matrix) ('Eisenstat').||'left'|
|If A is a matrix (not a string), and no preconditioner
has been specified, but a type of preconditioning has been specified then
RILU(Omega) will be used.
|If M is not specified then M = I (no preconditioning).
A preconditoner can be specified in the argument list:
... = CGSTAB(...,M,...)
... = CGSTAB(...,L,U,...)
... = CGSTAB(...,L,U,P,...)
... = CGSTAB(...,'M',...)
... = CGSTAB(...,'L','U'...)
as an N-by-N matrix M (then M is the preconditioner), or as N-by-N matrices L and U (then L*U is the preconditioner), or as N-by-N matrices L, U , and P (then P\(L*U) is the preconditioner), or as one or two strings containing the name of M-files ('M', or 'L' and 'U') which apply a linear operator to a given N-by-P array.
|CGSTAB (without input arguments) lists the options and the defaults.|