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 non-Hermitian, .... The algorithm is effective especially in case A is sparse and of large size.

The BiCGstab(ell) algorithm has been introduced in

G.L.G. Sleijpen and D.R. Fokkema,
BiCGstab(ell) for Linear Equations involving Unsymmetric Matrices with Complex Spectrum ,
ETNA, 1 (1993), pp. 11-32.
The BiCGstab(ell) here uses enhancements as explained in
G.L.G. Sleijpen and H.A. van der Vorst ,
Maintaining convergence properties of BiCGstab methods in finite precision arithmetic,
Numerical Algorithms, 10 (1995), pp. 203-223.

The Matlab® code here consists a number of function M-files that are bundled in one file cgstab.m. In this form,  cgstab.m requires Matlab version 5.1 or higher. We also provide a few very simple "test" files that can help to see how the cgstab.m can be used. cgstab.m.tar.gz contains the M-files cgstab.m plus the test files in tarred and zipped form (apply gunzip cgstab.m.tar.gz and tar xvf cgstab.m.tar to unpack). cgstab.m can also be obtained separately: see below, where you can also find a description of their use.

File: cgstab.m
Requires: Matlab Version 5.1.
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 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 1e-8. 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 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 = ZEROS(N,P)
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. 

  © Gerard L. G. Sleijpen   <>
  Last modified: Tuesday, 22-Jun-2010 13:19:49 CEST