The conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is symmetric and positive-definite. More...

#include <STK_CG.h>

Public Types
typedef ColVector::Type	Type

Public Member Functions
	CG ()
	Default Constructor.

	CG (MultFunctor const &mult, ColVector const &b, InitFunctor *const &p_init=0, Type eps=Arithmetic< Type >::epsilon())
	Constructor.

	CG (CG const &cg)
	Copy constructor.

	~CG ()
	destructor

CG *	clone () const
	clone pattern

ColVector const &	x () const

Real const &	x (int const &i) const

int const &	iter () const

int const &	nbStart () const

ColVector const &	r () const

void	setEps (Type const &eps)
	Set the tolerance.

void	setB (ColVector const &b)
	Set the constant vector.

void	setMultFunctor (MultFunctor const &mult)
	Set functor computing Ax.

void	setInitFunctor (InitFunctor *const &p_init)
	Set functor computing x at initialization.

int	run ()

String const &	error () const
	get the last error message.

Protected Member Functions
int	cg ()

Protected Attributes
String	msg_error_
	String with the last error message.

Private Attributes
ColVector	x_
	solution of the system

ColVector	r_
	residuals of the system

Type	eps_
	tolerance

int	iter_
	number of iterations

int	nbStart_
	number of restart_

MultFunctor const *	p_mult_
	pointer on the functor performing Ax

InitFunctor const *	p_init_
	pointer on the functor initializing x

ColVector const *	p_b_
	constant pointer on the second member of the system

Detailed Description

template<class MultFunctor, class ColVector, class InitFunctor = DefaultFunctor<ColVector>>
class STK::CG< MultFunctor, ColVector, InitFunctor >

The conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is symmetric and positive-definite.

The conjugate gradient method is an iterative method, so it can be applied to sparse systems that are too large to be handled by direct methods such as the Cholesky decomposition.

The algorithm is detailed below for solving $\mathbf{A}\mathbf{x} = \mathbf{b}$ where A is a real, symmetric, positive-definite matrix. The input vector x can be an approximate initial solution or 0.

$\begin{array}{l} \mathbf{r}_0 := \mathbf{b} - \mathbf{A x}_0 \\ \mathbf{p}_0 := \mathbf{r}_0 \\ k := 0 \\ \hbox{repeat} \\ \qquad \alpha_k := \frac{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k}{\mathbf{p}_k^\mathrm{T} \mathbf{A p}_k} \, \\ \qquad \mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k \mathbf{p}_k \\ \qquad \mathbf{r}_{k+1} := \mathbf{r}_k - \alpha_k \mathbf{A p}_k \\ \qquad \hbox{if } r_{k+1} \hbox{ is sufficiently small then exit loop} \\ \qquad \beta_k := \frac{\mathbf{r}_{k+1}^\mathrm{T} \mathbf{r}_{k+1}}{\mathbf{r}_k^\mathrm{T} \mathbf{r}_k} \\ \qquad \mathbf{p}_{k+1} := \mathbf{r}_{k+1} + \beta_k \mathbf{p}_k \\ \qquad k := k + 1 \\ \hbox{end repeat} \end{array}$

The result is $x_{k+1}$

Template Parameters

MultFunctor	A functor computing the result of Ax.
InitFunctor	A functor computing the initial value x.
ColVector	The type of the containers for the vectors.

Definition at line 86 of file STK_CG.h.

Member Typedef Documentation

◆ Type

template<class MultFunctor , class ColVector , class InitFunctor = DefaultFunctor<ColVector>>

typedef ColVector::Type STK::CG< MultFunctor, ColVector, InitFunctor >::Type

Definition at line 89 of file STK_CG.h.

Constructor & Destructor Documentation

◆ CG() [1/3]

template<class MultFunctor , class ColVector , class InitFunctor = DefaultFunctor<ColVector>>

STK::CG< MultFunctor, ColVector, InitFunctor >::CG ( )

inline

Default Constructor.

Definition at line 91 of file STK_CG.h.

91: x_(), r_(), eps_(0.), iter_(0), nbStart_(0), p_mult_(0), p_init_(0), p_b_(0) {}

STK::CG::nbStart_

int nbStart_

number of restart_

Definition STK_CG.h:213

STK::CG::r_

ColVector r_

residuals of the system

Definition STK_CG.h:207

STK::CG::eps_

Type eps_

tolerance

Definition STK_CG.h:209

STK::CG::p_mult_

MultFunctor const * p_mult_

pointer on the functor performing Ax

Definition STK_CG.h:215

STK::CG::p_b_

ColVector const * p_b_

constant pointer on the second member of the system

Definition STK_CG.h:219

STK::CG::p_init_

InitFunctor const * p_init_

pointer on the functor initializing x

Definition STK_CG.h:217

STK::CG::iter_

int iter_

number of iterations

Definition STK_CG.h:211

STK::CG::x_

ColVector x_

solution of the system

Definition STK_CG.h:205

Referenced by STK::CG< MultFunctor, ColVector, InitFunctor >::clone().

◆ CG() [2/3]

template<class MultFunctor , class ColVector , class InitFunctor = DefaultFunctor<ColVector>>

STK::CG< MultFunctor, ColVector, InitFunctor >::CG	(	MultFunctor const &	mult,
		ColVector const &	b,
		InitFunctor *const &	p_init = `0`,
		Type	eps = `Arithmetic<Type>::epsilon()`
	)

inline

Constructor.

Parameters

mult	functor which compute Ax with A a matrix and x a vector
b	from Ax=b
p_init	functor which initialize x
eps	tolerance

Definition at line 98 of file STK_CG.h.

      : x_(), r_()
      , eps_(eps), iter_(0), nbStart_(0)
      , p_mult_(&mult)
      , p_init_(p_init)
      , p_b_(&b)
    {}

◆ CG() [3/3]

template<class MultFunctor , class ColVector , class InitFunctor = DefaultFunctor<ColVector>>

STK::CG< MultFunctor, ColVector, InitFunctor >::CG ( CG< MultFunctor, ColVector, InitFunctor > const & cg )

inline

Copy constructor.

The constructor to copy.

Parameters

cg	the conjugate gradient to copy

Definition at line 112 of file STK_CG.h.

      : x_(cg.x_)
      , r_(cg.r_)
      , eps_(cg.eps_), iter_(cg.iter_), nbStart_(cg.nbStart_)
      , p_mult_(cg.p_mult_)
      , p_init_(cg.p_init_)
      , p_b_(cg.p_b_)
    {}

◆ ~CG()

template<class MultFunctor , class ColVector , class InitFunctor = DefaultFunctor<ColVector>>

STK::CG< MultFunctor, ColVector, InitFunctor >::~CG ( )

inline

destructor

Definition at line 121 of file STK_CG.h.

121{}

Member Function Documentation

◆ cg()

template<class MultFunctor , class ColVector , class InitFunctor = DefaultFunctor<ColVector>>

int STK::CG< MultFunctor, ColVector, InitFunctor >::cg ( )

inlineprotected

Returns: the number of iterations

Definition at line 156 of file STK_CG.h.

    {
      iter_ = 0;
      int nbStart = 0;
      ColVector xOld, z, p_;
 
      // initialization
      if(!p_init_) {x_ = *p_b_;}
      else { x_ = (*p_init_)();}
      while(nbStart<3)
      {
        int step = 0; //number of step
        Real bnorm2 = p_b_->norm2(), alpha, beta; //
        if (bnorm2 == 0.) bnorm2 = 1.;
        //compute the residuals
        r_ = *p_b_ - (*p_mult_)(x_);
        if (r_.norm2()/bnorm2 <eps_) { break;}
        //initialization of the conjugate direction
        p_= r_;
        while(1)
        {
          Real rnorm2 = r_.norm2();
          //compute z=A p
          z.move((*p_mult_)(p_));
          //compute alpha
          alpha = rnorm2/p_.dot(z);
          //update x_
          xOld.exchange(x_);
          x_ = xOld + alpha * p_;
          iter_++;
          //update residuals
          r_ = r_ - (alpha * z);
          //compute beta
          beta = 1/rnorm2;
          if ((rnorm2=r_.norm2())/bnorm2 <eps_) { nbStart = 2; break;}
          beta *= rnorm2;
          //update p_
          p_ = (p_ * beta) + r_;
          step++;
          if( step > 50 ) { break;}
        }
        nbStart++;
      }
      // return an error
      return iter_;
    }