Multigrid matrix inversion

\( A x = b\) More...

Collaboration diagram for Multigrid matrix inversion:

Classes
struct	dg::NestedGrids< Geometry, Matrix, Container >
	Hold nested grids and provide dg fast interpolation and projection matrices. More...
struct	dg::MultigridCG2d< Geometry, Matrix, Container >
	Solve. More...

Functions
template<class MatrixType0 , class ContainerType0 , class ContainerType1 , class MatrixType1 , class NestedGrids >
void	dg::nested_iterations (std::vector< MatrixType0 > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< MatrixType1 > &inverse_ops, NestedGrids &nested_grids)
	Full approximation nested iterations.
template<class NestedGrids , class MatrixType0 , class MatrixType1 , class MatrixType2 >
void	dg::multigrid_cycle (std::vector< MatrixType0 > &ops, std::vector< MatrixType1 > &inverse_ops_down, std::vector< MatrixType2 > &inverse_ops_up, NestedGrids &nested_grids, unsigned gamma, unsigned p)
	EXPERIMENTAL Full approximation multigrid cycle.
template<class MatrixType0 , class MatrixType1 , class MatrixType2 , class NestedGrids , class ContainerType0 , class ContainerType1 >
void	dg::full_multigrid (std::vector< MatrixType0 > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< MatrixType1 > &inverse_ops_down, std::vector< MatrixType2 > &inverse_ops_up, NestedGrids &nested_grids, unsigned gamma, unsigned mu)
	EXPERIMENTAL One Full multigrid cycle.
template<class NestedGrids , class MatrixType0 , class MatrixType1 , class MatrixType2 , class ContainerType0 , class ContainerType1 , class ContainerType2 >
void	dg::fmg_solve (std::vector< MatrixType0 > &ops, ContainerType0 &x, const ContainerType1 &b, std::vector< MatrixType1 > &inverse_ops_down, std::vector< MatrixType2 > &inverse_ops_up, NestedGrids &nested_grids, const ContainerType2 &weights, double eps, unsigned gamma)
	EXPERIMENTAL Full multigrid cycles.

Detailed Description

\( A x = b\)

Function Documentation

fmg_solve()

template<class NestedGrids , class MatrixType0 , class MatrixType1 , class MatrixType2 , class ContainerType0 , class ContainerType1 , class ContainerType2 >

void dg::fmg_solve	(	std::vector< MatrixType0 > &	ops,
		ContainerType0 &	x,
		const ContainerType1 &	b,
		std::vector< MatrixType1 > &	inverse_ops_down,
		std::vector< MatrixType2 > &	inverse_ops_up,
		NestedGrids &	nested_grids,
		const ContainerType2 &	weights,
		double	eps,
		unsigned	gamma )

EXPERIMENTAL Full multigrid cycles.

• Compute residual with given initial guess.
• If error larger than tolerance, do a full multigrid cycle with Chebeyshev iterations as smoother
• repeat

  Parameters

  ops	Index 0 is the MatrixType on the original grid, 1 on the half grid, 2 on the quarter grid, ...
  x	(read/write) contains initial guess on input and the solution on output
  b	The right hand side
  inverse_ops_down	a vector of inverse, smoothing operators (usually lambda functions combining operators and solvers) of size stages-1
  inverse_ops_up	a vector of inverse, smoothing operators (usually lambda functions combining operators and solvers) of size stages
  nested_grids	provides projection and interapolation operations and workspace
  weights	Defines the error norm
  eps	relative and absolute error tolerance
  gamma	The shape of the multigrid cycle: typically 1 (V-cycle) or 2 (W-cycle)

  Attention

  This method is rather unreliable, it only converges if the parameters are chosen correctly ( there need to be enough smooting steps for instance, and a large jump factor in the Elliptic class also seems to help) and otherwise just iterates to infinity. This behaviour is probably related to the use of the Chebyshev solver as a smoother

full_multigrid()

template<class MatrixType0 , class MatrixType1 , class MatrixType2 , class NestedGrids , class ContainerType0 , class ContainerType1 >

void dg::full_multigrid	(	std::vector< MatrixType0 > &	ops,
		ContainerType0 &	x,
		const ContainerType1 &	b,
		std::vector< MatrixType1 > &	inverse_ops_down,
		std::vector< MatrixType2 > &	inverse_ops_up,
		NestedGrids &	nested_grids,
		unsigned	gamma,
		unsigned	mu )

EXPERIMENTAL One Full multigrid cycle.

Parameters

ops	Index 0 is the f on the original grid, 1 on the half grid, 2 on the quarter grid, ...
x	(read/write) contains initial guess on input and the solution on output
b	The right hand side
inverse_ops_down	a vector of inverse, smoothing operators (usually lambda functions combining operators and solvers) of size stages-1
inverse_ops_up	a vector of inverse, smoothing operators (usually lambda functions combining operators and solvers) of size stages
nested_grids	provides projection and interapolation operations and workspace
gamma	The shape of the multigrid cycle: typically 1 (V-cycle) or 2 (W-cycle)
mu	The repetition of the multigrid cycle (1 is typically ok)

Attention

This method is rather unreliable, it only converges if the parameters are chosen correctly ( there need to be enough smooting steps for instance, and a large jump factor in the Elliptic class also seems to help) and otherwise just iterates to infinity. This behaviour is probably related to the use of the Chebyshev solver as a smoother

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ... If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also

See The dg dispatch system for a detailed explanation of our type dispatch system

multigrid_cycle()

template<class NestedGrids , class MatrixType0 , class MatrixType1 , class MatrixType2 >

void dg::multigrid_cycle	(	std::vector< MatrixType0 > &	ops,
		std::vector< MatrixType1 > &	inverse_ops_down,
		std::vector< MatrixType2 > &	inverse_ops_up,
		NestedGrids &	nested_grids,
		unsigned	gamma,
		unsigned	p )

EXPERIMENTAL Full approximation multigrid cycle.

See also

https://www.osti.gov/servlets/purl/15002749

Compute \( x_0^h \leftarrow C_{\nu_1}^{\nu_2}(x_0^h, b^h)\) for given \( b^h\) and initial guess \( x_0^h\):
Note that we need to save the initial guess \( w_0^h := x_0^h\) in a multigrid cycle because we need it to compute the error for the next upper level.

• Smooth \( f(x^h) = b^h\) with initial guess \(x_0^{h}\), overwrite \( x_0^h\)
• Compute \( r_0^h = b^h - f(x_0^h)\)
• Project \( r_0^{2h} = P r_0^h\) and \( x_0^{2h} = Px_0^h\)
• Compute \( b^{2h} = f(x_0^{2h}) + r_0^{2h}\)
• \( w_0^{2h} = x_0^{2h}\)
• recursively call itself \(\gamma\) times or solve \( f(x^{2h}) = b^{2h}\) with initial guess \( x_0^{2h}\), overwrite \( x_0^{2h}\)
  • ...
• \( \delta^{2h} = x_0^{2h} - w_0^{2h}\) (Here we need the saved initial guess \( w_0^{2h}\))
• \( x_0^h = x_0^h + I \delta^{2h}\)
• Smooth \( f(x^h) = b^h\) with initial guess \( x_0^h\), overwrite \( x_0^h\)

This algorithm forms the core of multigrid algorithms.

Parameters

ops	Index 0 is the Operator on the original grid, 1 on the half grid, 2 on the quarter grid, ...
inverse_ops_down	a vector of inverse, smoothing operators (usually lambda functions combining operators and solvers) of size stages-1
inverse_ops_up	a vector of inverse, smoothing operators (usually lambda functions combining operators and solvers) of size stages
nested_grids	provides projection and interapolation operations and workspace
gamma	The shape of the multigrid cycle: typically 1 (V-cycle) or 2 (W-cycle)
p	The current stage h

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ... If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also

See The dg dispatch system for a detailed explanation of our type dispatch system

nested_iterations()

template<class MatrixType0 , class ContainerType0 , class ContainerType1 , class MatrixType1 , class NestedGrids >

void dg::nested_iterations	(	std::vector< MatrixType0 > &	ops,
		ContainerType0 &	x,
		const ContainerType1 &	b,
		std::vector< MatrixType1 > &	inverse_ops,
		NestedGrids &	nested_grids )

Full approximation nested iterations.

Solve \( f(x_0^{h}) = b^{h}\) for given \( b^h\) and initial guess \( x_0^h\):

• Compute \( r_0^h = b^h - f(x_0^h)\)
• Project \( r_0^{2h} = P r_0^h\) and \( x_0^{2h} = Px_0^h\)
• Compute \( b^{2h} = f(x_0^{2h}) + r_0^{2h}\)
• Cycle: \( f(x^{2h}) = b^{2h}\) with initial guess \( x_0^{2h}\):
  • (residuum \( r_0^{2h}\) remains unchanged)
  • \( r_0^{4h} = Pr_0^{2h}\) and \( x_0^{4h} = Px_0^{2h}\)
  • Compute \( b^{4h} = f(x_0^{4h}) + r_0^{4h}\)
  • Cycle (on coarsest grid solve): \( f(x^{4h}) = b^{4h}\) with initial guess \(x_0^{4h}\):
    • ...
  • \( \delta^{4h} = x^{4h} - x_0^{4h}\)
  • \( \tilde x_0^{2h} = x_0^{2h} + I \delta^{4h}\)
  • Solve \( f(x^{2h}) = b^{2h}\) with initial guess \( \tilde x_0^{2h}\)
• \( \delta^{2h} = x^{2h} - x_0^{2h}\)
• \( \tilde x_0^h = x_0^h + I \delta^{2h}\)
• Solve \( f(x^h) = b^h\) with initial guess \( \tilde x_0^h\)

This algorithm is equivalent to a multigrid V-cycle with zero down-grid smoothing and infinite (i.e. solving) upgrid smoothing.

Parameters

ops	Operators f on the various grids, i.e. dg::apply( ops[0], x, b) computes b = f(x). Index 0 is on the original grid, 1 on the half grid, 2 on the quarter grid, ...
x	(read/write) contains initial guess on input and the solution on output (if the initial guess is good enough the solve may return immediately)
b	The right hand side
inverse_ops	a vector of inverse operators f^{-1}, i.e. dg::apply( inverse_ops[0], b, x) computes \( x = f^{-1}(b)\) (usually lambda functions combining operators and solvers). On call x contains the initial guess and should contain the solution on return.
nested_grids	provides projection and interapolation operations and workspace

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ... If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also

See The dg dispatch system for a detailed explanation of our type dispatch system