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 :

SOR

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:

SOR

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 :

Relaxation_MatVect.cxx

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 :

BiCg.cxx

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 :

BiCgcr.cxx

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 :

BiCgStab.cxx

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 :

BiCgStabl.cxx

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 :

Cg.cxx

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 :

Cgne.cxx

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 :

Cgs.cxx

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 :

CoCg.cxx

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 :

Gcr.cxx

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 :

Gmres.cxx

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 :

Lsqr.cxx

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 :

MinRes.cxx

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 :

QCgs.cxx

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 :

Qmr.cxx

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 :

QmrSym.cxx

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 :

Symmlq.cxx

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.

Location :

TfQmr.cxx

Email Contact