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-A*X0)*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(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. 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 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   <G.L.G.Sleijpen@uu.nl>