Gerard L.G. Sleijpen

BiCGstab(ell)

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.
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 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) where
	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-AX0)RAND(P,1). The OPTIONS structure specifies certain parameters in the algorithm.
	FIELD NAME			PARAMETER	DEFAULT
	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
	OPTIONS. TypeAcc			specifies in what sense the systems of equations are to be solved 0: each equation with a relative accuracy Tol: NORM(AX(:,I)-B(:,I)) < TolNORM(B(:,I)) for all I=1:P; 1: each equation with an accuracy of Tol: NORM(AX(:,I)-B(:,I)) < Tol for all I=1:P; 2: each equation in the space with a relative accuracy Tol: NORM(AXMU-BMU) < TolNORM(BMU) for all MU of size P-by-1.	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 right-hand 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 (=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'
	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 N-by-N matrix M (then M is the preconditioner), or as N-by-N matrices L and U (then LU is the preconditioner), or as N-by-N matrices L, U , and P (then P\(LU) 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.