From this page you can get a Matlab® implementation of the BiCGstab(ell)
algorithm.
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 nonHermitian, .... The algorithm
is effective especially in case A is sparse and of large size.
The BiCGstab(ell) algorithm has been introduced in
File:  cgstab.m  
Requires:  Matlab Version 5.1.  
Function:  [X,HISTORY] = CGSTAB(A,B,X0,TR0,OPTIONS,L,U);  
Description:  
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 NbyP,
where A is NbyN.


CGSTAB will start iterating from
an initial guess which by default is an array of size NbyP
of
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 1e8. The default maximum
number of iterations is 300. No preconditioning is used. CGSTAB
uses
BiCGstab(ell)
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 Mfile 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) where 

X0  An NbyP array of doubles, an initial guess.  
TR0  An Nby1 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
= ZEROS(N,P).
If TR0 is not specified then TR0 = (BA*X0)*RAND(P,1). The OPTIONS structure specifies certain parameters in the algorithm. 

FIELD NAME  PARAMETER  DEFAULT  
OPTIONS.Tol  Tol is the convergence tolerance.  1e8  
OPTIONS.MaxIt  Maximum number of iterations.  300  
OPTIONS.Disp  Shows size of intermediate residuals ('yes'/'no').  'no'  
OPTIONS.ell  ell for BiCGstab(ell).  4  
OPTIONS.
TypeAcc 
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*MUB*MU) < Tol*NORM(B*MU) for all MU of size Pby1. 
0  
OPTIONS.
Simultane 
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 righthand side B ) ('no').  'yes'  
OPTIONS.
TypePrecond 
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 (=twosided 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'  
OPTIONS.
Omega 
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.

0.97  
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 NbyN matrix M (then M is the preconditioner), or as NbyN matrices L and U (then L*U is the preconditioner), or as NbyN matrices L, U , and P (then P\(L*U) is the preconditioner), or as one or two strings containing the name of Mfiles ('M', or 'L' and 'U') which apply a linear operator to a given NbyP array. 

CGSTAB (without input arguments) lists the options and the defaults.  