Iterative Solvers

The solvers library comes with some iterative solver methods which target a linear system.

These solvers can be divided into two kinds of group. One is for symmetric linear system and another is for a non-symmetric system. Each solver uses largeMatrix and Vector as inputs and returns the result in a Vector.

The type of problem, real or complex, needs specified also as another input of a solver. In some cases, a solver can take advantage of a Preconditioner as an extra input to decrease the conditioner number of the system. A Preconditioner, which is also a LargeMatrix involved with matrix factorization, is only available in the type of SkylineStorage and CsStorage.

All solvers in the library are based on the class IterativeSolver, which provides some basic properties and methods. Specific iterative algorithms are implemented in the inherited class for two cases: with and without a preconditioner.

The IterativeSolver class

As the base class of all other solver classes, the IterativeSolver provides some basic properties and methods for all descendants.

class IterativeSolver
{
  public:
    //! Constructor by name
    IterativeSolver(const string_t& name, number_t vb = 0);
    //! Constructor by type
    IterativeSolver(IterativeSolverType ist, number_t vb = 0);
    //! Full constructor without type
    IterativeSolver(const string_t& name, number_t maxOfIt, real_t eps, number_t vb = 0, bool prec = true);
    //! Full constructor without name
    IterativeSolver(IterativeSolverType ist, number_t maxOfIt, real_t eps, number_t vb = 0, bool prec = true);
    //! Full constructor
    IterativeSolver(const string_t& name, IterativeSolverType ist, number_t maxOfIt, real_t eps, number_t vb = 0, bool prec = true);

    //! Destructor
    virtual ~IterativeSolver();

    //! Iterative solver name for documentation purposes
    string_t name() { return name_; }
    //! Reset Solver
    void resetSolver();

   protected:
   //! Define a default value for the maximum number of iterations
   number_t maximumOfIterations(const size_t nbRows);
   ...
}

A solver object can be constructed in several ways: With default values or with user-input values.

Information of solver can be printed out with some auxiliary functions.

void printHeader(const size_t) const;
void printHeader(const size_t, const String&) const;
void printHeader(const size_t, const Real) const;
void printHeader(const size_t, const Number) const;
void printHeader(const size_t, const Number, const String&) const;
void printOutput() const;
void printIteration() const;

In order to avoid invalid cast from complex to real and vice versa, some inline functions are provided.

inline void assign(Real& x, const Real& y) { x = y; }
inline void assign(Complex& x, const Real& y) { x = y; }
inline void assign(Complex& x, const Complex& y) { x = y; }
inline void assign(Real& x, const Complex& y) { x = y.real(); }

See also

• library= solvers,

• header= IterativeSolver.hpp,

• implementation= IterativeSolver.cpp,

• test= test_BicgSolver.hpp, test_BicgStabSolver.hpp, test_CgSolver.hpp, test_CgsSolver.hpp, test_GmresSolver.hpp, test_QmrSolver.hpp, test_SorSolver.hpp, test_SsorSolver.hpp

• header dep= Preconditioner.hpp, config.h, utils.h

The Preconditioner class

In some cases, the rate of convergence of an iterative solver can increase as a result of preconditioning. A preconditioner is basically a special matrix which can be factorized. There exist several types of matrix-factorization for Preconditioner: LU, LDLt, LDL*, Diag and SOR. Because of constraints of factorizing, only matrices of SkylineStorage and CsStorage can be used as a preconditioner.

In detail, a preconditioner is based on a template class Mat, on which all the preconditioned-solver are processed.

template<class Mat>
class Preconditioner
{
 public:
  Mat* precondMatrix_p; // pointer to precondition matrix
 private:
  //! Type of preconditioner
  PreconditionerType type_;
  //! Relaxation parameter for SSOR preconditioners
  Real omega_;
 public:
  Preconditioner()
   : precondMatrix_p(0), type_(_noPreconditioner), omega_(0.)  { }

   Preconditioner(PreconditionerType pt, Mat& A, const Real w = 1.)
    : precondMatrix_p(&A), type_(pt), omega_(w)  { }

   Preconditioner(Mat& A, PreconditionerType pt, const Real w = 1.)
     : precondMatrix_p(&A), type_(pt), omega_(w)  {  }

   ~Preconditioner() {};
   //! Return type of preconditionner
   static String name(PreconditionerType pt)
   {
     String s("");
     switch ( pt )
     {
       case _lu:       s = "LU"; bre
       ak;
       case _ldlt:     s = "LDLt"; break;
       case _ldlStar:  s = "LDL*"; break;
       case _ssor:     s = "SSOR"; break;
       case _diag:     s = "Diagonal"; break;
        //case myPreconditioner: s = "User Supplied"; break;
       default: break;
     }
     return s;
  }
 ...

Hint

A preconditioner can only templated with a specific kind of matrix that can be decomposed. There are only two types of storage, with which matrix can be factorized: SkylineStorage with LU, LDLt, LDL* factorization; CsStorage with SOR and DIAG factorization.

There are several constructors:

Preconditioner()
 : precondMatrix_p(0), type_(_noPreconditioner), omega_(0.)  { }

Preconditioner(PreconditionerType pt, Mat& A, const Real w = 1.)
 : precondMatrix_p(&A), type_(pt), omega_(w)  {  }

Preconditioner(Mat& A, PreconditionerType pt, const Real w = 1.)
 : precondMatrix_p(&A), type_(pt), omega_(w)  {  }

The PreconditionerType is an enumeration.

enum PreconditionerType {_noPreconditioner, _lu, _ldlt, _ldlStar, _ssor, _diag, _myPreconditioner};

See also

• library= solvers,

• header= Preconditioner.hpp,

• implementation= Preconditioner.cpp,

• test= test_LargeMatrixSkylineStorage.cpp, test_LargeMatrixCsStorage.cpp

• header dep= config.h, utils.h

Solvers

Each solver class in the solvers library implements a specific iterative method to solve a linear system: A*X=B in two cases: Real or Complex. The algorithm of all solvers except SOR and SSOR can work with and without a Preconditioner. Each solver takes largeMatrix \(matA\) and Vector \(vecB\) as inputs and a Vector \(vecX\) as output whose value is initialized with the input Vector \(vecX0\). A solver without a preconditioner is invoked with the overloading operator():

template<class Mat, class VecB, class VecX>
VecX operator()(Mat& matA, VecB& vecB, VecX& vecX0, ValueType solType)

Meanwhile, solver with a preconditioner \(pc\) can be used with the call

template<class Mat, class VecB, class VecX, class Prec>
VecX operator()(Mat& matA, VecB& vecB, VecX& vecX0, Prec& pc, ValueType solType)

User must specify the type of solver solType, which can be _real or _complex. The algorithm corresponds to without-preconditioner solver

template<typename K, class Mat, class VecB, class VecX>
void algorithm(Mat& matA, VecB& vecB, VecX& vecX, VecX& vecR)

and to a preconditioned one

template<typename K, class Mat, class VecB, class VecX, class Prec>
void algorithm(Mat& matA, VecB& vecB, VecX& vecX, VecX& vecR, Prec& pc)

The BicgSolver class

This class implements the Biconjugate Gradient method to solve a linear system. There are several basic constructors calling IterativeSolver constructors

class BicgSolver : public IterativeSolver
{
 public:
  //! Default Constructor
  BicgSolver() : IterativeSolver(_bicg) {}
  //! Full constructor
  BicgSolver(real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
  : IterativeSolver(_bicg, maxOfIt, eps, vb) {}

  //! contructors with key-value system
  BicgSolver(const Parameter& p1) : IterativeSolver(_bicg)
  {
    std::vector<Parameter> ps(1, p1);
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_bicg;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  BicgSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_bicg)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_bicg;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  BicgSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_bicg)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_bicg;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_BicgSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.hpp

The BicgStabSolver class

This class implements the Biconjugate Gradient Stabilized method to solve a linear system. There are several basic constructors calling IterativeSolver constructors

class BicgStabSolver : public IterativeSolver
{
  public:
   //! Default Constructor
   BicgStabSolver() : IterativeSolver(_bicgstab) {}
   //! Full Constructor
   BicgStabSolver(real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
   : IterativeSolver(_bicgstab, maxOfIt, eps, vb) {}

   //! contructors with key-value system
   BicgStabSolver(const Parameter& p1) : IterativeSolver(_bicgstab)
   {
    std::vector<Parameter> ps(1, p1);
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_bicgstab;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
   }
   BicgStabSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_bicgstab)
   {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_bicgstab;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
   }
   BicgStabSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_bicgstab)
   {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_bicgstab;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
   }
  ...
}

See also

• library= solvers,

• test= test_BicgStabSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.hpp

The CgSolver class

This class implements the Conjugate Gradient method to solve a linear system. There are several basic constructors calling IterativeSolver constructors

class CgSolver : public IterativeSolver
{
 public:
  //! Default constructor
  CgSolver(): IterativeSolver(_cg) {}
  //! Full constructor
  CgSolver(real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
  : IterativeSolver(_cg, maxOfIt, eps, vb) {}

  //! contructors with key-value system
  CgSolver(const Parameter& p1) : IterativeSolver(_cg)
  {
    std::vector<Parameter> ps(1, p1);
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_cg;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  CgSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_cg)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_cg;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  CgSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_cg)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_cg;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_CgSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.hpp

The CgsSolver class

This class implements the Conjugate Gradient Squared method to solve a linear system. There are several basic constructors calling IterativeSolver constructors

class CgsSolver : public IterativeSolver
{
 public:
  //! Default constructor
  CgsSolver()
    : IterativeSolver(_cgs) {}
  //! Full constructor
  CgsSolver(real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
    : IterativeSolver(_cgs, maxOfIt, eps, vb) {}

  //! contructors with key-value system
  CgsSolver(const Parameter& p1) : IterativeSolver(_cgs)
  {
    std::vector<Parameter> ps(1, p1);
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_cgs;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  CgsSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_cgs)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_cgs;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  CgsSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_cgs)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_cgs;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_CgsSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.H2-approximation

The GmresSolver class

This class implements the Generalized Minimal Residual method to solve a linear system. Different from other solvers, the GmresSolver has more constructors: user can specify the Krylov dimension in the constructor.

class GmresSolver : public IterativeSolver
{
 private:
  number_t krylovDim_;  //!< Dimension of Krylov space

 public:
  //! Constructor with Krylov dimension
  GmresSolver(number_t kd=defaultKrylovDimension)
    : IterativeSolver(_gmres, defaultMaxIterations, theDefaultConvergenceThreshold),
      krylovDim_(kd) {}
  //! Full constructor
  GmresSolver(number_t kd, real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
  : IterativeSolver(_gmres, maxOfIt, eps, vb), krylovDim_(kd) {}

  //! contructors with key-value system
  GmresSolver(const Parameter& p1) : IterativeSolver(_gmres)
  {
    std::vector<Parameter> ps(1, p1);
    real_t omega=1.;
    IterativeSolverType ist=_gmres;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim_, verboseLevel_, name_, ist);
  }
  GmresSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_gmres)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    real_t omega=1.;
    IterativeSolverType ist=_gmres;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim_, verboseLevel_, name_, ist);
  }
  GmresSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_gmres)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    real_t omega=1.;
    IterativeSolverType ist=_gmres;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim_, verboseLevel_, name_, ist);
  }
  GmresSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4) : IterativeSolver(_gmres)
  {
    std::vector<Parameter> ps(4);
    ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4;
    real_t omega=1.;
    IterativeSolverType ist=_gmres;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim_, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_GmresSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.H2-approximation

The QmrSolver class

This class implements the Quasi-Minimal Residual method to solve a linear system. There are several basic constructors calling IterativeSolver constructors

class QmrSolver : public IterativeSolver
{
 public:
  //! Constructor
  QmrSolver()
    : IterativeSolver(_qmr) {}
  //! Full Constructor
  QmrSolver(real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
  : IterativeSolver(_qmr, maxOfIt, eps, vb) {}

  //! contructors with key-value system
  QmrSolver(const Parameter& p1) : IterativeSolver(_qmr)
  {
    std::vector<Parameter> ps(1, p1);
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_qmr;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  QmrSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_qmr)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_qmr;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  QmrSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_qmr)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    real_t omega=1.; number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_qmr;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega, krylovDim, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_QmrSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.H2-approximation

The SorSolver class

This class implements the Successive Over Relaxation method to solve a linear system. For this solver, user can have more choices to specify the relaxation factor \(omega\) in the constructor. Different from other solvers above, SorSolver doesn’t work with a preconditioner.

class SorSolver : public IterativeSolver
{
 private:
  real_t omega_;  //!< relaxation factor of SOR method

 public:
  //! Constructor with omega
  SorSolver(real_t w=1.) : IterativeSolver(_sor), omega_(w) {}
  //! Full constructor
  SorSolver(real_t w, real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
  : IterativeSolver(_sor, maxOfIt, eps, vb), omega_(w) {}

  //! contructors with key-value system
  SorSolver(const Parameter& p1) : IterativeSolver(_sor)
  {
    std::vector<Parameter> ps(1, p1);
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_sor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  SorSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_sor)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_sor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  SorSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_sor)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_sor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  SorSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4) : IterativeSolver(_sor)
  {
    std::vector<Parameter> ps(4);
    ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4;
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_sor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_SorSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.H2-approximation

The SsorSolver class

This class implements the Symmetric Successive Over Relaxation method to solve a linear system. Like SorSolver, the SsorSolver only works without a preconditioner and user can specify the relaxation factor \(omega\) in the constructor.

class SsorSolver : public IterativeSolver
{
 private:
  real_t omega_;  //!< relaxation factor of SOR method

 public:
  //! Constructor with omega
  SsorSolver(real_t w=1.) : IterativeSolver(_ssor), omega_(w) {}
  //! Full constructor
  SsorSolver(real_t w, real_t eps, number_t maxOfIt = defaultMaxIterations, number_t vb = 0)
  : IterativeSolver(_ssor, maxOfIt, eps, vb), omega_(w) {}

  //! contructors with key-value system
  SsorSolver(const Parameter& p1) : IterativeSolver(_ssor)
  {
    std::vector<Parameter> ps(1, p1);
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_ssor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  SsorSolver(const Parameter& p1, const Parameter& p2) : IterativeSolver(_ssor)
  {
    std::vector<Parameter> ps(2);
    ps[0]=p1; ps[1]=p2;
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_ssor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  SsorSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3) : IterativeSolver(_ssor)
  {
    std::vector<Parameter> ps(3);
    ps[0]=p1; ps[1]=p2; ps[2]=p3;
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_ssor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  SsorSolver(const Parameter& p1, const Parameter& p2, const Parameter& p3, const Parameter& p4) : IterativeSolver(_ssor)
  {
    std::vector<Parameter> ps(4);
    ps[0]=p1; ps[1]=p2; ps[2]=p3; ps[3]=p4;
    number_t krylovDim=defaultKrylovDimension;
    IterativeSolverType ist=_ssor;
    buildSolverParams(ps, epsilon_, maxOfIterations_, omega_, krylovDim, verboseLevel_, name_, ist);
  }
  ...
}

See also

• library= solvers,

• test= test_SsorSolver.cpp

• header dep= IterativeSolver.hpp, Preconditioner.H2-approximation