Iterative solvers
Syntax
The syntax of all iterative solvers is the same
int Gmres(const Matrix& A, Vector& x, const Vector& b,
Preconditioner& M, Iteration& iter);
The first argument A is the matrix to solve. The type of this argument is template,
therefore you can provide your own class defining the matrix. The second and third argument
are vectors, x contains the initial guess on input, the solution on output,
b contains the right-hand-side. Again, the type of x and b is template,
so you can have your own class to define the vectors. The fourth template argument is the
preconditioner. For the moment, there is an implementation of the identity preconditioner
(no preconditioning) and a SOR (Successive Over Relaxation)
preconditioner. The last argument is a Seldon structure defining the parameters of the iteration.
Basic use
We provide an example of iterative resolution using seldon structures for the matrices and the vectors.
// first we construct a sparse matrix int n = 1000; // size of linear system // we assume that you know how to fill arrays values, ptr, ind Matrix<double, General, RowSparse< A(n, n, nnz, values, ptr, ind); // then we declare vectors Vector<double> b(n), x(n); // you fill right-hand side and initial guess b.Fill(); x.Fill(0); // initialization of iteration parameters int nb_max_iter = 1000; // relative stopping criterion double tolerance = 1e-6; Iteration<double> iter(nb_max_iter, tolerance); // identity preconditioner -> no preconditioning Preconditioner_Base precond; // then you can call an iterative solver, Cg for example Cg(A, x, b, precond, iter); // if you are using Gmres, you can set the restart parameter // by default, this parameter is equal to 10 iter.SetRestart(5); Gmres(A, x, b, precond, iter);
Advanced use
Now, we will show an example where we construct a new class for a matrix that we don't store. If you want to test iterative solvers with that kind of strategy, you can take a look at the file test/program/iterative_test.cpp. The main thing to do is to overload the functions Mlt and MltAdd that are called by iterative solvers to perform the matrix-vector product.
// Class for a new type of matrix.
template<class T>
class BlackBoxMatrix
{
protected :
// internal data used to represent the matrix
int n;
T beta;
Seldon::Vector<T> lambda;
public :
// basic constructor
BlackBoxMatrix(int n_, const T& beta_)
{
beta = beta_; n = n_;
lambda.Reallocate(n);
for (int i = 0; i < n; i++)
lambda(i) = i+1;
}
// this method is used by iterative solvers
// this should return the number of rows of the matrix
int GetM() const
{
return n;
}
};
// We need to define MltAdd(alpha, A, X, beta, Y) : Y <- beta Y + alpha A X
template<class T, class Vector>
void MltAdd(const T& alpha, const BlackBoxMatrix<T>& A,
const Vector& X, const T& beta, Vector& Y)
{
if (beta == T(0))
Y.Fill(T(0));
else
Mlt(beta, Y);
int n = A.GetM();
Vector C(X);
// C = S^{-1} C
for (int i = (n-2); i >= 0; i--)
C(i) -= A.beta*C(i+1);
// C = B C
for (int i = 0; i < n; i++)
C(i) *= A.lambda(i);
// C = S C
for (int i = 0; i < (n-1); i++)
C(i) += A.beta*C(i+1);
// Y = Y + alpha C
Add(alpha, C, Y);
}
// and the transpose...
template<class T, class Vector>
void MltAdd(const T& alpha, const class_SeldonTrans& Trans,
const BlackBoxMatrix<T>& A,
const Vector& X, const T& beta, Vector& Y)
{
if (beta == T(0))
Y.Fill(T(0));
else
Mlt(beta, Y);
int n = A.GetM();
// Transpose of S B S^{-1} is S^{-t} B S^t
Vector C(X);
// Y = S^t Y
for (int i = (n-1); i >= 1; i--)
C(i) += A.beta*C(i-1);
// Y = B Y
for (int i = 0; i < n; i++)
C(i) *= A.lambda(i);
// Y = S^{-t} Y
for (int i = 1; i < n; i++)
C(i) -= A.beta*C(i-1);
// Y = Y + alpha C
Add(alpha, C, Y);
}
template<class T, class Vector>
void Mlt(const BlackBoxMatrix<T>& A, const Vector& X, Vector& Y)
{
Y.Zero();
MltAdd(T(1), A, X, T(0), Y);
}
template<class T, class Vector>
void Mlt(const class_SeldonTrans& Trans, const BlackBoxMatrix<T>& A,
const Vector& X, Vector& Y)
{
Y.Zero();
MltAdd(T(1), Trans, A, X, T(0), Y);
}
// class for preconditioning
template<class T>
class MyPreconditioner : public BlackBoxMatrix<T>
{
public :
MyPreconditioner(int n, T beta) : BlackBoxMatrix<T>(n, beta) { }
// solving M r = z
template<class Matrix, class Vector>
void Solve(const Matrix& A, const Vector& r, Vector& z)
{
Mlt(*this, r, z);
}
// solving transpose(M) r = z
template<class Matrix, class Vector>
void TransSolve(const Matrix& A, const Vector& r, Vector& z)
{
Mlt(SeldonTrans, *this, r, z);
}
};
int main()
{
// now it is very classical like the previous example
int n = 20; double beta = 0.5;
BlackBoxMatrix<double>; A(n, beta);
Vector<double> b(n), x(n);
// you fill right-hand side and initial guess
b.Fill();
x.Fill(0);
// initialization of iteration parameters
int nb_max_iter = 1000;
double tolerance = 1e-6;
Iteration<double> iter(nb_max_iter, tolerance);
// your own preconditioner
MyPreconditioner<double> precond(n, 1.2);
// then you can call an iterative solver, Qmr for example
Qmr(A, x, b, precond, iter);
return 0;
}
If you want to use your own class of vector, there are other functions to define : Add, Norm2, DotProd, DotProdConj , Copy, Copy and the method Zero.
Methods of Preconditioner_Base:
| Solve | Applies the preconditioner |
| TransSolve | Applies the transpose of the preconditioner |
Methods of SorPreconditioner:
| InitSymmetricPreconditioning | Symmetric SOR will be used |
| InitUnsymmetricPreconditioning | SOR will be used |
| SetParameterRelaxation | changes the relaxation parameter |
| SetNumberIterations | changes the number of iterations |
| Solve | Applies the preconditioner |
| TransSolve | Applies the transpose of the preconditioner |
Methods of Iteration:
| Constructors | |
| GetRestart | returns restart parameter |
| SetRestart | changes restart parameter |
| GetTolerance | returns stopping criterion |
| SetTolerance | changes stopping criterion |
| SetMaxIterationNumber | changes maximum number of iterations |
| GetIterationNumber | returns iteration number |
| SetIterationNumber | changes iteration number |
| ShowMessages | displays residual each 100 iterations |
| HideMessages | displays nothing |
| ShowFullHistory | displays residual each iteration |
| Init | provides right hand side |
| First | returns true for the first iteration |
| IsInitGuess_Null | returns true if the initial guess is zero |
| SetInitGuess | informs if the initial guess is zero or no |
| Finished | returns true if the stopping criteria are satisfied |
| Fail | informs that the iterative solver failed |
| ErrorCode | returns error code |
Functions :
| SOR | Performs SOR iterations |
| BiCg | BIConjugate Gradient |
| BiCgcr | BIConjugate Gradient Conjugate Residual |
| BiCgStab | BIConjugate Gradient STABilized |
| BiCgStabl | BIConjugate Gradient STABilized (L) |
| Cg | Conjugate Gradient |
| Cgne | Conjugate Gradient Normal Equation |
| Cgs | Conjugate Gradient Squared |
| CoCg | Conjugate Orthogonal Conjugate Gradient |
| Gcr | Generalized Conjugate Residual |
| Gmres | Generalized Minimum RESidual |
| Lsqr | Least SQuaRes |
| MinRes | Minimum RESidual |
| QCgs | Quasi Conjugate Gradient Squared |
| Qmr | Quasi Minimum Residual |
| QmrSym | Quasi Minimum Residual SYMmetric |
| Symmlq | SYMMetric Least sQuares |
| TfQmr | Transpose Free Quasi Minimum Residual |
Solve, TransSolve for Preconditioner_Base
Syntax :
void Solve(const Matrix&, const Vector&, Vector&); void TransSolve(const Matrix&, const Vector&, Vector&);
These methods should be overloaded if you want to define your own preconditioner since they define the application of the preconditioner and its transpose to a vector.
Example :
// constructor of a matrix int n = 20; Matrix<double> A(n, n); A.Fill(); // declaration of a preconditioner Preconditioner_Base M; // vectors Vector<double> r(n), z(n); r.Fill(); // now we apply preconditioning, i.e. solving M z = r // for Preconditioner_Base, it will set z = r (identity preconditioner) M.Solve(A, r, z); // we can also apply the transpose of preconditioner // i.e. solving transpose(M) z = r M.TransSolve(A, r, z);
Location :
Class Preconditioner_Base
Iterative.hxx
Iterative.cxx
Solve, TransSolve for SorPreconditioner
Syntax :
void Solve(const Matrix&, const Vector&, Vector&); void TransSolve(const Matrix&, const Vector&, Vector&);
These methods define the application of the preconditioner and its transpose to a vector. The used preconditioner is SOR (Successive Over Relaxation). It can be used in its symmetric version (SSOR), or the unsymmetric version. In this last case, the application of the preconditioner consists of a forward sweep while the transpose consists of a backward sweep. If you are using Seldon structures of sparse matrices, the function SOR is available. If you are using other structures, it is necessary to overload the function SOR.
Example :
// constructor of a matrix int n = 20; Matrix<double> A(n, n); A.Fill(); // declaration of a preconditioner SorPreconditioner<double> M; // by default, relaxation parameter omega = 1 // number of iterations = 1 // you can change that M.SetParameterRelaxation(1.5); M.SetNumberIterations(2); // if you prefer to use symmetric version M.InitSymmetricPreconditioning(); // vectors Vector<double> r(n), z(n); r.Fill(); // now we apply preconditioning, i.e. solving M z = r // for Preconditioner_Base, it will set z = r (identity preconditioner) M.Solve(A, r, z); // we can also apply the transpose of preconditioner // i.e. solving transpose(M) z = r M.TransSolve(A, r, z);
Related topics :
Location :
Class SorPreconditioner
Precond_Ssor.cxx
InitSymmetricPreconditioning, InitUnSymmetricPreconditioning for SorPreconditioner
Syntax :
void InitSymmetricPreconditioning(); void InitUnsymmetricPreconditioning();
InitSymmetricPreconditioning sets SSOR as preconditioning. The symmetric
SOR consists of a forward sweep followed by a backward sweep. InitSymmetricPreconditioning
sets SOR, it consists of a forward sweep for the preconditioner, and a backward sweep
for the transpose of the preconditioner.
Example :
// declaration of a preconditioner SorPreconditioner<double> M; // by default, symmetric preconditioning M.InitUnsymmetricPreconditioning();
Related topics:
Location :
Class SorPreconditioner
Precond_Ssor.cxx
SetParameterRelaxation for SorPreconditioner
Syntax :
void SetParameterRelaxation(const T& omega);
This methods changes the relaxation parameter.
Location :
Class SorPreconditioner
Precond_Ssor.cxx
SetNumberIterations for SorPreconditioner
Syntax :
void SetNumberIterations(int);
This methods changes the number of SOR iterations.
Location :
Class SorPreconditioner
Precond_Ssor.cxx
Constructors for Iteration
Syntax :
Iteration(); Iteration(int, const T&);
The second constructor specifies the maximum number of iterations and the stopping criterion
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
GetRestart, SetRestart for Iteration
Syntax :
int GetRestart(); void SetRestart(int);
These methods give access to the restart parameter, which is used by some iterative solvers, e.g. BiCgSTAB(l), Gmres and Gcr. The default value is equal to 10.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
GetTolerance, SetTolerance for Iteration
Syntax :
T GetTolerance(); void SetTolerance(T);
These methods give access to the stopping criterion.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
SetMaxIterationNumber for Iteration
Syntax :
void SetMaxIterationNumber(int);
This method gives access to the maximum iteration number.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
GetIterationNumber, SetIterationNumber for Iteration
Syntax :
int GetIterationNumber(); void SetIterationNumber(int);
This method gives access to the iteration number.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
ShowMessages, HideMessages, ShowFullHistory for Iteration
Syntax :
void ShowMessages(); void HideMessages(); void ShowFullHistory();
If ShowMessages is called, it will display residual each 100 iterations
during the resolution. If HideMessages is called, there is no display at all,
whereas ShowFullHistory displays residual each iteration.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
Init for Iteration
Syntax :
void Init(const Vector&);
This method is used by iterative solvers to initialize the computation of residuals. Since relative residual is computed, we need to know the norm of the first residual.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
First for Iteration
Syntax :
bool First();
This method returns true for the first iteration.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
SetInitGuess, IsInitGuess_Null for Iteration
Syntax :
bool IsInitGuess_Null(); void SetInitGuess(bool);
SetInitGuess allows you to inform the iterative solver
that your initial guess is null. This can spare one matrix-vector product, which is
good if the expected number of iterations is small.
Example :
int n = 10; Matrix<double> A(n, n); Vector<double> x(n), b(n); A.Fill(); b.Fill(); Iteration<double> iter(100, 1e-6); // you inform that initial guess is null iter.SetInitGuess(true); // if the initial guess is null // and x is non-null, the solver enforces x to be 0 Preconditioner_Base M; Gmres(A, x, b, M, iter);
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
Finished for Iteration
Syntax :
bool Finished(const Vector&);
This method is used by iterative solvers to know if the stopping criterion is reached. This method also displays the residual if required.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
Fail for Iteration
Syntax :
bool Fail(int, const string&);
This method is used by iterative solvers when a breakdown occurs, often due to a division by 0.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
ErrorCode for Iteration
Syntax :
int ErrorCode();
This method returns the error code after the resolution. If the resolution has been successful, it should return 0.
Location :
Class Iteration
Iterative.hxx
Iterative.cxx
SOR
Syntax :
void SOR(const Matrix&, Vector&, const Vector&, const T&, int); void SOR(const Matrix&, Vector&, const Vector&, const T&, int, int);
This method tries to solve A x = b with n iterations
of SOR.
Example :
int n = 1000; Matrix<double, General, RowSparse> A; // you initialize A as you want (SetData for example) // then vectors Vector<double> x(n), b(n); x.Zero(); b.Fill(); // we want to solve A x = b // 2 Sor iterations (forward sweep) with omega = 0.5 double omega = 0.5; int nb_iterations = 2; SOR(A, x, b, omega, nb_iterations); // you can ask for symmetric SOR: forward sweeps followed // by backward sweeps SOR(A, x, b, omega, nb_iterations, 0); // if you need backward sweep SOR(A, x, b, omega, nb_iterations, 3);
Location :
BiCg
Syntax :
int BiCg(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICG algorithm.
This algorithm can solve complex general linear systems and calls matrix-vector
products with the transpose matrix.
Location :
BiCgcr
Syntax :
int BiCgcr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGCR algorithm.
This algorithm can solve symmetric complex linear systems.
Location :
BiCgStab
Syntax :
int BiCgStab(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
BiCgStabl
Syntax :
int BiCgStabl(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB(l) algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
Cg
Syntax :
int Cg(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using CG algorithm.
This algorithm can solve real symmetric or hermitian linear systems.
Location :
Cgne
Syntax :
int Cgne(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using CGNE algorithm.
This algorithm can solve complex general linear systems and calls matrix
vector products with the transpose matrix.
Location :
Cgs
Syntax :
int Cgs(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using CGS algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
Cocg
Syntax :
int Cocg(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm.
This algorithm can solve complex symmetric linear systems.
Location :
Gcr
Syntax :
int Gcr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
Gmres
Syntax :
int Gmres(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using restarted GMRES algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
Lsqr
Syntax :
int Lsqr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using LSQR algorithm.
This algorithm can solve complex general linear systems and calls matrix
vector products with the transpose matrix.
Location :
MinRes
Syntax :
int MinRes(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using BICGSTAB algorithm.
This algorithm can solve complex symmetric linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
QCgs
Syntax :
int QCgs(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using QCGS algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
Location :
Qmr
Syntax :
int Qmr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using QMR algorithm.
This algorithm can solve complex general linear systems and calls matrix
vector products with the transpose matrix.
Location :
QmrSym
Syntax :
int QmrSym(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using symmetric QMR algorithm.
This algorithm can solve complex symmetric linear systems.
Location :
Symmlq
Syntax :
int Symmlq(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using SYMMLQ algorithm.
This algorithm can solve complex symmetric linear systems.
Location :
TfQmr
Syntax :
int TfQmr(const Matrix&, Vector&, const Vector&,
Preconditioner&, Iteration&);
This method tries to solve A x = b by using TFQMR algorithm.
This algorithm can solve complex general linear systems and doesn't call matrix
vector products with the transpose matrix.
