Matrix-free version of the GaussSeidel iterative solver. More...

#include <ProductGaussSeidel.hpp>

Inheritance diagram for bogus::ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt >:

Public Types
enum	{ has_trivial_diagonal = IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value }

typedef GaussSeidelBase < ProductGaussSeidel, BlockMatrixType >	Base

typedef Base::GlobalProblemTraits	GlobalProblemTraits

typedef GlobalProblemTraits::Scalar	Scalar

typedef block_solvers_impl::DiagonalMatrixWrapper < DiagonalType,!!has_trivial_diagonal >	DiagWrapper

typedef block_solvers_impl::DMtStorage < BlockMatrixType, DiagWrapper, PrecomputeDMt >	DMtStorage

typedef LocalProblemTraits < Base::BlockTraits::RowsPerBlock, Scalar >	BlockProblemTraits

typedef BlockMatrixTraits < BlockMatrixType >	BlockTraits

typedef Signal< unsigned, Scalar >	CallBackType

Public Member Functions
	ProductGaussSeidel ()
	Default constructor – you will have to call setMatrix() before using the solve() function.

	ProductGaussSeidel (const BlockObjectBase< BlockMatrixType > &matrix)
	Constructor with the main system matrix M.

	ProductGaussSeidel (const BlockObjectBase< BlockMatrixType > &matrix, const DiagonalType &diagonal)
	Constructor with both the main matrix M and the diagonal D.

ProductGaussSeidel &	setMatrix (const BlockObjectBase< BlockMatrixType > &matrix)
	Sets the system matrix (M) and initializes internal structures.

ProductGaussSeidel &	setDiagonal (const DiagonalType &diagonal)
	Sets the system diagonal (D) and initializes internal structures.

template<typename NSLaw , typename RhsT , typename ResT >
Scalar	solve (const NSLaw &law, const RhsT &b, ResT &x, bool tryZeroAsWell=true) const
	Finds an approximate solution for a constrained linear problem.

template<typename NSLaw , typename RhsT , typename ResT , typename LSDerived , typename HDerived >
Scalar	solveWithLinearConstraints (const NSLaw &law, const BlockObjectBase< LSDerived > &Cinv, const BlockObjectBase< HDerived > &H, const Scalar alpha, const RhsT &b, const RhsT &c, ResT &x, bool tryZeroAsWell=true, unsigned solveEvery=1) const

template<typename NSLaw , typename RhsT , typename ResT , typename WDerived >
Scalar	solveWithLinearConstraints (const NSLaw &law, const BlockObjectBase< WDerived > &W, const RhsT &b, ResT &x, bool tryZeroAsWell=true, unsigned solveEvery=1) const

void	setMaxThreads (unsigned maxThreads=0)
	Sets the maximum number of threads that the solver can use. More...

void	setAutoRegularization (Scalar maxRegul)
	Sets the auto-regularization (a.k.a. proximal point) coefficient. More...

void	setEvalEvery (unsigned evalEvery)
	Sets the number of iterations that should be performed between successive evaluations of the global error function. More...

void	setSkipTol (Scalar skipTol)
	Sets the minimum iteration step size under which local problems are temporarily frozen.

void	setSkipIters (unsigned skipIters)
	Sets the number of iterations for temporarily freezing local problems.

void	useInfinityNorm (bool useInfNorm)
	Sets whether the solver will use the infinity norm instead of the l1 one to compute the global residual from the local ones.

bool	usesInfinityNorm () const

Scalar	eval (const NSLaw &law, const ResT &y, const RhsT &x) const
	Eval the current global residual as a function of the local ones. More...

void	projectOnConstraints (const NSLaw &projector, VectorT &x) const
	Projects the variable `x` on the constraints defined by `projector`.

void	dualityCOV (const NSLaw &law, const RhsT &b, ResT &x) const
	Compute associated change of variable (see NSLaw)

Scalar	solve (const NSLaw &law, const RhsT &b, ResT &x) const

void	setMaxIters (unsigned maxIters)
	For iterative solvers: sets the maximum number of iterations.

unsigned	maxIters () const

void	setTol (Scalar tol)
	For iterative solvers: sets the solver tolerance.

Scalar	tol () const

CallBackType &	callback ()
	Callback hook; will be triggered every N iterations, depending on the solver. More...

const CallBackType &	callback () const

const BlockObjectBase < BlockMatrixType > &	matrix () const

Protected Types
typedef Base::Index	Index

typedef Base::BlockProblemTraits::Matrix	DiagonalBlockType

Protected Member Functions
void	updateLocalMatrices ()

template<typename NSLaw , typename VecT , typename ResT >
void	innerLoop (bool parallelize, const NSLaw &law, const VecT &b, std::vector< unsigned char > &skip, Scalar &ndxRef, VecT &Mx, ResT &x) const

const BlockMatrixBase < BlockMatrixType > &	explicitMatrix () const

const IterableBlockObject < BlockMatrixType > &	iterableMatrix () const

void	processLocalMatrices ()

Scalar	evalAndKeepBest (const NSLaw &law, const ResT &x, const typename GlobalProblemTraits::DynVector &y, typename GlobalProblemTraits::DynVector &x_best, Scalar &err_best) const

bool	tryZero (const NSLaw &law, const RhsT &b, ResT &x, typename GlobalProblemTraits::DynVector &x_best, Scalar &err_best) const

void	updateScalings ()

Protected Attributes
DiagWrapper	m_diagonal

DMtStorage	m_DMt

ResizableSequenceContainer < DiagonalBlockType >::Type	m_localMatrices

GlobalProblemTraits::DynVector	m_regularization

unsigned	m_maxThreads
	See setMaxThreads(). Defaults to 0 .

unsigned	m_evalEvery
	See setEvalEvery(). Defaults to 25.

Scalar	m_skipTol
	See setSkipTol(). Defaults to 1.e-6.

unsigned	m_skipIters
	See setSkipIters() Defaults to 10.

Scalar	m_autoRegularization

GlobalProblemTraits::DynVector	m_scaling

bool	m_useInfinityNorm
	See useInfinityNorm(). Defaults to false.

const BlockObjectBase < BlockMatrixType > *	m_matrix
	Pointer to the matrix of the system.

unsigned	m_maxIters
	See setMaxIters()

Scalar	m_tol
	See setTol()

CallBackType	m_callback

Detailed Description

template<typename BlockMatrixType, typename DiagonalType = typename BlockMatrixType::Scalar, bool PrecomputeDMt = !( IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value )>
class bogus::ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt >

Matrix-free version of the GaussSeidel iterative solver.

Assumes that the system matrix is defined as the product (M D M'), with D a block-diagonal matrix whose block sizes coincide with those of the columns of M.

Template Parameters

BlockMatrixType	The type of the main solver matrix M (the constructor's first argument).
DiagonalType	The type of the diagonal matrix D. The default type (M's Scalar type) means using a constant times the identity matrix.
PrecompupeDMt	Whether to precompute and store the (D M') part of the product (M D M'). Defaults to true if DiagonalType is not the default value.

Warning: Parallelization is supported, but dangerous. If in doubt, use setMaxThreads(1)

Note: D must be block-diagonal. Works best when M is row-major and its columns are quite sparse

See Also: GaussSeidel

Member Function Documentation

template<typename BlockMatrixType >

CallBackType& bogus::BlockSolverBase< BlockMatrixType >::callback ( )

inherited

Callback hook; will be triggered every N iterations, depending on the solver.

Useful to monitor the convergence of the solver. Can be connected to a function that takes an unsigned and a Scalar as parameters. The first argument will be the current iteration number, and the second the current residual.

See Also: Signal< unsigned, Scalar >

Scalar bogus::ConstrainedSolverBase< ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt > , BlockMatrixType >::eval	(	const NSLaw &	law,
		const ResT &	y,
		const RhsT &	x
	)		const

inherited

Eval the current global residual as a function of the local ones.

y should be such that y = m_matrix * x + rhs

Returns

the current residual err defined as follow :

if m_useInfinityNorm is true, then err $:= \max\limits_{1 \leq i \leq n } law.eval(i,x_i,y_i)$
else err := $\frac 1 {n+1} \sum\limits_{1 \leq i \leq n } law.eval(i,x_i,y_i)$

void bogus::GaussSeidelBase< ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt > , BlockMatrixType >::setAutoRegularization ( Scalar maxRegul )

inherited

Sets the auto-regularization (a.k.a. proximal point) coefficient.

The regularization works by slightly altering the local problems, so at each iteration we try to solve

$\left\{ \begin{array}{rcl} y^{k+1} &=& ( M + \alpha I ) x^{k+1} - \alpha x^k + b^{k} \\ &s.t.& law (x^{k+1},y^{k+1}) \end{array} \right.$

where $\alpha$ is the regularization coefficient.

Note that as $| x^{k+1} - x^{k} | \rightarrow 0$ when the algorithm converges, we are still trying to find a solution of the same global problem.

For under-determined problems, regularization might helps preventing x reaching problematically high values. Setting $\alpha$ to a too big value will however degrade the convergence of the global algorithm.

Parameters

maxRegul If greater than zero, then positive terms will be added to the diagonal of the local matrices so that their minimal eigenvalue become greater than maxRegul.

void bogus::GaussSeidelBase< ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt > , BlockMatrixType >::setEvalEvery ( unsigned evalEvery )

inherited

Sets the number of iterations that should be performed between successive evaluations of the global error function.

( Those evaluations require a full matrix/vector product, and are therfore quite costly )

void bogus::GaussSeidelBase< ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt > , BlockMatrixType >::setMaxThreads ( unsigned maxThreads = 0 )

inherited

Sets the maximum number of threads that the solver can use.

If \p maxThreads is zero, then it will use the current OpenMP setting.

Warning: If multi-threading is enabled without coloring, the result will not be deterministic, as it will depends on the order in which threads solve contacts.

On the other hand, the algorithm will run much faster.

template<typename BlockMatrixType, typename DiagonalType = typename BlockMatrixType::Scalar, bool PrecomputeDMt = !( IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value )>

template<typename NSLaw , typename RhsT , typename ResT , typename LSDerived , typename HDerived >

Scalar bogus::ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt >::solveWithLinearConstraints	(	const NSLaw &	law,
		const BlockObjectBase< LSDerived > &	Cinv,
		const BlockObjectBase< HDerived > &	H,
		const Scalar	alpha,
		const RhsT &	b,
		const RhsT &	c,
		ResT &	x,
		bool	tryZeroAsWell = `true`,
		unsigned	solveEvery = `1`
	)		const

Solves

$\left\{ \begin{array}{rclll} y &=& MDM^T x &+ p &+ b \\ 0 &=& H^T x &+ \frac{1}{\alpha}C p &+ c \\ &s.t.& law (x,y) \end{array} \right.$

where Cinv is such that Cinv * x = $C^{-1} x$

Warning: Requires: m_evalEvery multiple of solveEvery ;

template<typename BlockMatrixType, typename DiagonalType = typename BlockMatrixType::Scalar, bool PrecomputeDMt = !( IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value )>

template<typename NSLaw , typename RhsT , typename ResT , typename WDerived >

Scalar bogus::ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt >::solveWithLinearConstraints	(	const NSLaw &	law,
		const BlockObjectBase< WDerived > &	W,
		const RhsT &	b,
		ResT &	x,
		bool	tryZeroAsWell = `true`,
		unsigned	solveEvery = `1`
	)		const

Solves

$\left\{ \begin{array}{rcl} y &=& MDM^T x + W x + b \\ &s.t.& law (x,y) \end{array} \right.$

with W arbitrary linear operator ( matrix or expression )

Warning: Requires: m_evalEvery multiple of solveEvery ;

Member Data Documentation

Scalar bogus::GaussSeidelBase< ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt > , BlockMatrixType >::m_autoRegularization

protectedinherited

See Also: setAutoRegularization(). Defaults to 0.

The documentation for this class was generated from the following file:

ProductGaussSeidel.hpp

Public Types

Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Detailed Description

template<typename BlockMatrixType, typename DiagonalType = typename BlockMatrixType::Scalar, bool PrecomputeDMt = !( IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value )> class bogus::ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt >

Member Function Documentation

Member Data Documentation

template<typename BlockMatrixType, typename DiagonalType = typename BlockMatrixType::Scalar, bool PrecomputeDMt = !( IsSame<DiagonalType, typename BlockMatrixType::Scalar>::Value )>
class bogus::ProductGaussSeidel< BlockMatrixType, DiagonalType, PrecomputeDMt >