Algebraic multigrid solver for SpMatrix / AlgVector systems. More...

#include <AMReX_AlgMG.H>

Public Types
using	TS = std::conditional_t< std::is_same_v< T, double >, float, T >

using	BottomSolver = AlgMGBottomSolver

using	KrylovSolver = AlgMGKrylovSolver

using	InterpType = AlgMGInterpType

using	Smoother = AlgMGSmoother

Public Member Functions
	AlgMG ()=default

	AlgMG (SpMatrix< T > &a_mat)

void	define (SpMatrix< T > &a_mat)
	Define the solver with a square matrix. Must be called before solve.

void	solve (AlgVector< T > &a_sol, AlgVector< T > const &a_rhs)
	Solve A x = b. `a_sol` holds the initial guess on entry.

void	setVerbose (int v)

void	setPrintIndentation (std::string s)
	Prefix of the lines the solver prints, e.g. to nest them in MLMG output.

void	setMaxIter (int n)

void	setFixedIter (int n)
	If positive, run exactly this many V-cycles regardless of tolerance.

void	setRelTol (T t)

void	setAbsTol (T t)
	Converged when the residual 2-norm is below max(reltol*bnorm, abstol).

void	setSmoother (Smoother sm)
	Takes effect at the next solve; the hierarchy is kept.

void	setRelaxWeight (T w)
	Relaxation weight of the smoother; negative restores the default. Next solve.

void	setChebyshevDegree (int d)
	Chebyshev polynomial degree. Next solve.

void	setChebyshevRatio (T r)
	Chebyshev eigenvalue ratio lambda_max / lambda_min. Next solve.

void	setPreSmooth (int nu)
	Number of pre-smoothing sweeps; negative restores the smoother's default.

void	setPostSmooth (int nu)
	Number of post-smoothing sweeps; negative restores the smoother's default.

void	setBottomSmooth (int nu)

void	setStrongThreshold (T theta)

void	setMaxLevels (int n)

void	setMaxCoarseSize (Long n)
	Stop coarsening when a level has at most this many rows. Next setup.

void	setAggressiveNumLevels (int n)
	Number of levels with aggressive coarsening, starting at the finest. Next setup.

void	setAggressiveDirectInterp (bool b)

void	setSingular (bool b)

void	setInterpType (InterpType it)
	Takes effect at the next setup.

void	setPMaxElmts (int n)

void	setTruncFactor (T f)
	Drop interpolation weights below this fraction of the row maximum. Next setup.

void	setKrylovSolver (KrylovSolver a)
	Krylov solver with the V-cycle as preconditioner (see KrylovSolver). Default None.

void	precond (AlgVector< T > &x, AlgVector< T > const &b)

void	setBottomSolver (BottomSolver bs)

void	setBottomTol (T t)

void	setBottomMaxIter (int n)

void	setBottomVerbose (int v)

void	setThrowException (bool b)

int	getNumIters () const

int	getNumSetups () const
	Number of setups run so far (for checks that the setup is reused).

T	getResidualNorm () const

int	numLevels () const
	Number of levels; 0 before setup.

std::size_t	getSetupPeakBytes () const
	Peak bytes in use in The_Arena during setup (0 if not available).

Long	numGlobalRows (int lev) const
	Global number of rows on a level. Valid after setup.

SpMatrix< T > const &	getMatrix (int lev) const

SpMatrix< T > const &	getInterp (int lev) const

void	setup ()

bool	coarsen_level (int lev, SpMatrix< T > &A, SpMatrix< T > &P, AlgPartition &cpart, bool &aggressive, Long &ncoarse, bool &truncated, Long &nnz_untruncated)

SpMatrix< TS >	create_soc (SpMatrix< T > const &A, T scale) const

void	pmis (SpMatrix< TS > &S, SpMatrix< TS > &ST, Gpu::DeviceVector< int > &state, bool isolated_as_coarse=false) const

void	project_out (AlgVector< T > &x)

SpMatrix< T >	interpolation (SpMatrix< T > &A, SpMatrix< TS > &S, SpMatrix< TS > const &ST, Gpu::DeviceVector< int > const &state, Gpu::DeviceVector< Long > const &cidx, AlgPartition const &cpart, T scale) const
	Interpolation of the selected type from the C-points in `state`.

Static Public Member Functions
static T	strength_scale (SpMatrix< T > const &A)
	Largest \|a_ii\|, the unit of the strength matrix entries.

static AlgPartition	coarse_numbering (Gpu::DeviceVector< int > const &state, Gpu::DeviceVector< Long > &cidx, int which=PointType::coarse)
	Local index among the rows of type `which` and their partition.

static SpMatrix< T >	build_interp (SpMatrix< T > &A, SpMatrix< TS > const &S, Gpu::DeviceVector< int > const &state, Gpu::DeviceVector< Long > const &cidx, AlgPartition const &cpart)
	Classical direct interpolation.

static SpMatrix< T >	truncate_interp (SpMatrix< T > &P, AlgPartition const &cpart, int max_elmts, T trunc_factor)

static void	normalize_rows (SpMatrix< T > &P)
	Scales the rows of P to sum to one; aborts on an empty row.

static SpMatrix< TS >	second_pass_strength (SpMatrix< TS > &S, Gpu::DeviceVector< int > const &state, Gpu::DeviceVector< Long > const &cidx, AlgPartition const &cpart)

static SpMatrix< T >	select_coarse_rows (SpMatrix< T > const &P, Gpu::DeviceVector< int > const &state, Gpu::DeviceVector< Long > const &cidx, AlgPartition const &cpart)
	Rows of P at the C-points, in C-point numbering.

static SpMatrix< T >	build_interp_mm (SpMatrix< T > &A, SpMatrix< TS > &S, SpMatrix< TS > const &ST, Gpu::DeviceVector< int > const &state, Gpu::DeviceVector< Long > const &cidx, AlgPartition const &cpart, T scale, bool plus_i, int const row_state=nullptr, int max_elmts=0, T trunc_factor=T(0), Long nnz_untruncated=nullptr)

Static Public Attributes
static constexpr int	max_p_elmts = 32

static constexpr int	default_p_max_elmts = 4

Detailed Description

template<typename T>
class amrex::AlgMG< T >

Algebraic multigrid solver for SpMatrix / AlgVector systems.

The setup selects coarse points with PMIS [2] on the classical strength of connection [1], builds the interpolation P (extended+i [3] in the matrix-matrix form of [4] by default, or direct [1]) with at most p_max_elmts entries per row [3], and forms A_c = P^T A P until the coarsest level has at most max_coarse_size rows. The first aggressive_num_levels levels use aggressive coarsening [5]. The solve runs V-cycles with Chebyshev [7], l1-Jacobi [6], weighted Jacobi or, on CPUs, l1 hybrid Gauss-Seidel [6] smoothing, on their own or as the preconditioner of a Krylov solver (setKrylovSolver).

Singular matrices with the constant null vector (setSingular) are handled by removing the mean of the right-hand side and of the coarsest-level residuals.

The matrix must be square with rows sorted by column index (see SpMatrix::sortCSR), and it must stay alive while the solver is in use. With R = P^T, the coarsening assumes a symmetric or nearly symmetric matrix. Parameters are set per object with the setters below; the Linear Solvers chapter of the documentation describes their effect.

References (full citations in the documentation): [1] Ruge & Stüben, Algebraic multigrid, SIAM 1987. doi:10.1137/1.9781611971057.ch4 [2] De Sterck, Yang & Heys, SIAM J. Matrix Anal. Appl. 27 (2006). doi:10.1137/040615729 [3] De Sterck, Falgout, Nolting & Yang, Numer. Linear Algebra Appl. 15 (2008). doi:10.1002/nla.559 [4] Li, Sjögreen & Yang, SIAM J. Sci. Comput. 43 (2021). doi:10.1137/20M134931X [5] Yang, Numer. Linear Algebra Appl. 17 (2010). doi:10.1002/nla.689 [6] Baker, Falgout, Kolev & Yang, SIAM J. Sci. Comput. 33 (2011). doi:10.1137/100798806 [7] Adams, Brezina, Hu & Tuminaro, J. Comput. Phys. 188 (2003). doi:10.1016/S0021-9991(03)00194-3

Member Typedef Documentation

◆ BottomSolver

template<typename T >

using amrex::AlgMG< T >::BottomSolver = AlgMGBottomSolver

◆ InterpType

template<typename T >

using amrex::AlgMG< T >::InterpType = AlgMGInterpType

◆ KrylovSolver

template<typename T >

using amrex::AlgMG< T >::KrylovSolver = AlgMGKrylovSolver

◆ Smoother

template<typename T >

using amrex::AlgMG< T >::Smoother = AlgMGSmoother

◆ TS

template<typename T >

using amrex::AlgMG< T >::TS = std::conditional_t<std::is_same_v<T,double>, float, T>

Value type of the strength matrices: float when T is double. Their entries are the matrix entries divided by strength_scale.

Constructor & Destructor Documentation

◆ AlgMG() [1/2]

template<typename T >

amrex::AlgMG< T >::AlgMG ( )

default

◆ AlgMG() [2/2]

template<typename T >

amrex::AlgMG< T >::AlgMG ( SpMatrix< T > & a_mat )

explicit

Member Function Documentation

◆ build_interp()

template<typename T >

SpMatrix< T > amrex::AlgMG< T >::build_interp	(	SpMatrix< T > &	A,
		SpMatrix< TS > const &	S,
		Gpu::DeviceVector< int > const &	state,
		Gpu::DeviceVector< Long > const &	cidx,
		AlgPartition const &	cpart
	)

static

Classical direct interpolation.

◆ build_interp_mm()

template<typename T >

SpMatrix< T > amrex::AlgMG< T >::build_interp_mm	(	SpMatrix< T > &	A,
		SpMatrix< TS > &	S,
		SpMatrix< TS > const &	ST,
		Gpu::DeviceVector< int > const &	state,
		Gpu::DeviceVector< Long > const &	cidx,
		AlgPartition const &	cpart,
		T	scale,
		bool	plus_i,
		int const *	row_state = `nullptr`,
		int	max_elmts = `0`,
		T	trunc_factor = `T(0)`,
		Long *	nnz_untruncated = `nullptr`
	)

static

Extended (plus_i = false) or extended+i interpolation in matrix-matrix form. ST is the transpose of S, scale its unit. With row_state, only its C-points get a row of P; the other rows are empty. In CPU builds, P is also truncated as in truncate_interp if max_elmts or trunc_factor is positive, and nnz_untruncated (if not null) gets the number of weights before.

◆ coarse_numbering()

template<typename T >

AlgPartition amrex::AlgMG< T >::coarse_numbering	(	Gpu::DeviceVector< int > const &	state,
		Gpu::DeviceVector< Long > &	cidx,
		int	which = `PointType::coarse`
	)

static

Local index among the rows of type which and their partition.

◆ coarsen_level()

template<typename T >

bool amrex::AlgMG< T >::coarsen_level	(	int	lev,
		SpMatrix< T > &	A,
		SpMatrix< T > &	P,
		AlgPartition &	cpart,
		bool &	aggressive,
		Long &	ncoarse,
		bool &	truncated,
		Long &	nnz_untruncated
	)

Selects the C-points of level lev and builds its interpolation P (rows partitioned by cpart). Returns false when coarsening stops. The strength matrices are freed on return, before the triple product. If truncated, P is already truncated and nnz_untruncated holds its number of weights before.

◆ create_soc()

template<typename T >

SpMatrix< typename AlgMG< T >::TS > amrex::AlgMG< T >::create_soc	(	SpMatrix< T > const &	A,
		T	scale
	)		const

Strong part of A: S_ij = a_ij / scale if i strongly depends on j, else absent. scale is strength_scale(A).

◆ define()

template<typename T >

void amrex::AlgMG< T >::define ( SpMatrix< T > & a_mat )

Define the solver with a square matrix. Must be called before solve.

◆ getInterp()

template<typename T >

SpMatrix< T > const & amrex::AlgMG< T >::getInterp ( int lev ) const

inline

◆ getMatrix()

template<typename T >

SpMatrix< T > const & amrex::AlgMG< T >::getMatrix ( int lev ) const

inline

◆ getNumIters()

template<typename T >

int amrex::AlgMG< T >::getNumIters ( ) const

inline

◆ getNumSetups()

template<typename T >

int amrex::AlgMG< T >::getNumSetups ( ) const

inline

Number of setups run so far (for checks that the setup is reused).

◆ getResidualNorm()

template<typename T >

T amrex::AlgMG< T >::getResidualNorm ( ) const

inline

◆ getSetupPeakBytes()

template<typename T >

std::size_t amrex::AlgMG< T >::getSetupPeakBytes ( ) const

inline

Peak bytes in use in The_Arena during setup (0 if not available).

◆ interpolation()

template<typename T >

SpMatrix< T > amrex::AlgMG< T >::interpolation	(	SpMatrix< T > &	A,
		SpMatrix< TS > &	S,
		SpMatrix< TS > const &	ST,
		Gpu::DeviceVector< int > const &	state,
		Gpu::DeviceVector< Long > const &	cidx,
		AlgPartition const &	cpart,
		T	scale
	)		const

Interpolation of the selected type from the C-points in state.

◆ normalize_rows()

template<typename T >

void amrex::AlgMG< T >::normalize_rows ( SpMatrix< T > & P )

static

Scales the rows of P to sum to one; aborts on an empty row.

◆ numGlobalRows()

template<typename T >

Long amrex::AlgMG< T >::numGlobalRows ( int lev ) const

inline

Global number of rows on a level. Valid after setup.

◆ numLevels()

template<typename T >

int amrex::AlgMG< T >::numLevels ( ) const

inline

Number of levels; 0 before setup.

◆ pmis()

template<typename T >

void amrex::AlgMG< T >::pmis	(	SpMatrix< TS > &	S,
		SpMatrix< TS > &	ST,
		Gpu::DeviceVector< int > &	state,
		bool	isolated_as_coarse = `false`
	)		const

PMIS C/F selection. state gets PointType values. Points without any strong connection become F, or C if isolated_as_coarse.

◆ precond()

template<typename T >

void amrex::AlgMG< T >::precond	(	AlgVector< T > &	x,
		AlgVector< T > const &	b
	)

Apply one V-cycle to b with zero initial guess: x = M^{-1} b. This is the preconditioner used by the Krylov solver; usable stand-alone after the first solve or setup.

◆ project_out()

template<typename T >

void amrex::AlgMG< T >::project_out ( AlgVector< T > & x )

inline

Singular systems: remove the mean of x (the constant null vector is the same on every level since interpolation preserves constants).

◆ second_pass_strength()

template<typename T >

SpMatrix< typename AlgMG< T >::TS > amrex::AlgMG< T >::second_pass_strength	(	SpMatrix< TS > &	S,
		Gpu::DeviceVector< int > const &	state,
		Gpu::DeviceVector< Long > const &	cidx,
		AlgPartition const &	cpart
	)

static

Strength graph between C-points connected through a strong path of length two, in C-point numbering (rows partitioned by cpart).

◆ select_coarse_rows()

template<typename T >

SpMatrix< T > amrex::AlgMG< T >::select_coarse_rows	(	SpMatrix< T > const &	P,
		Gpu::DeviceVector< int > const &	state,
		Gpu::DeviceVector< Long > const &	cidx,
		AlgPartition const &	cpart
	)

static

Rows of P at the C-points, in C-point numbering.

◆ setAbsTol()

template<typename T >

void amrex::AlgMG< T >::setAbsTol ( T t )

inline

Converged when the residual 2-norm is below max(reltol*bnorm, abstol).

◆ setAggressiveDirectInterp()

template<typename T >

void amrex::AlgMG< T >::setAggressiveDirectInterp ( bool b )

inline

Use classical direct interpolation for the first stage of the two-stage aggressive interpolation (cheaper setup, more cycles).

◆ setAggressiveNumLevels()

template<typename T >

void amrex::AlgMG< T >::setAggressiveNumLevels ( int n )

inline

Number of levels with aggressive coarsening, starting at the finest. Next setup.

◆ setBottomMaxIter()

template<typename T >

void amrex::AlgMG< T >::setBottomMaxIter ( int n )

inline

◆ setBottomSmooth()

template<typename T >

void amrex::AlgMG< T >::setBottomSmooth ( int nu )

inline

Number of smoother sweeps on the coarsest level with the Jacobi bottom solver; negative restores the default.

◆ setBottomSolver()

template<typename T >

void amrex::AlgMG< T >::setBottomSolver ( BottomSolver bs )

inline

Coarsest-level solver (see BottomSolver). Default direct, which uses smoother sweeps instead beyond 1024 rows.

◆ setBottomTol()

template<typename T >

void amrex::AlgMG< T >::setBottomTol ( T t )

inline

◆ setBottomVerbose()

template<typename T >

void amrex::AlgMG< T >::setBottomVerbose ( int v )

inline

◆ setChebyshevDegree()

template<typename T >

void amrex::AlgMG< T >::setChebyshevDegree ( int d )

inline

Chebyshev polynomial degree. Next solve.

◆ setChebyshevRatio()

template<typename T >

void amrex::AlgMG< T >::setChebyshevRatio ( T r )

inline

Chebyshev eigenvalue ratio lambda_max / lambda_min. Next solve.

◆ setFixedIter()

template<typename T >

void amrex::AlgMG< T >::setFixedIter ( int n )

inline

If positive, run exactly this many V-cycles regardless of tolerance.

◆ setInterpType()

template<typename T >

void amrex::AlgMG< T >::setInterpType ( InterpType it )

inline

Takes effect at the next setup.

◆ setKrylovSolver()

template<typename T >

void amrex::AlgMG< T >::setKrylovSolver ( KrylovSolver a )

inline

Krylov solver with the V-cycle as preconditioner (see KrylovSolver). Default None.

◆ setMaxCoarseSize()

template<typename T >

void amrex::AlgMG< T >::setMaxCoarseSize ( Long n )

inline

Stop coarsening when a level has at most this many rows. Next setup.

◆ setMaxIter()

template<typename T >

void amrex::AlgMG< T >::setMaxIter ( int n )

inline

◆ setMaxLevels()

template<typename T >

void amrex::AlgMG< T >::setMaxLevels ( int n )

inline

Maximum number of levels. Two means a two-level method. Takes effect at the next setup, i.e., the next solve after define.

◆ setPMaxElmts()

template<typename T >

void amrex::AlgMG< T >::setPMaxElmts ( int n )

inline

Maximum number of entries per row of P (0: no limit; at most max_p_elmts; negative restores the default). Next setup.

◆ setPostSmooth()

template<typename T >

void amrex::AlgMG< T >::setPostSmooth ( int nu )

inline

Number of post-smoothing sweeps; negative restores the smoother's default.

◆ setPreSmooth()

template<typename T >

void amrex::AlgMG< T >::setPreSmooth ( int nu )

inline

Number of pre-smoothing sweeps; negative restores the smoother's default.

◆ setPrintIndentation()

template<typename T >

void amrex::AlgMG< T >::setPrintIndentation ( std::string s )

inline

Prefix of the lines the solver prints, e.g. to nest them in MLMG output.

◆ setRelaxWeight()

template<typename T >

void amrex::AlgMG< T >::setRelaxWeight ( T w )

inline

Relaxation weight of the smoother; negative restores the default. Next solve.

◆ setRelTol()

template<typename T >

void amrex::AlgMG< T >::setRelTol ( T t )

inline

◆ setSingular()

template<typename T >

void amrex::AlgMG< T >::setSingular ( bool b )

inline

Singular matrix with the constant vector in its null space (e.g. Poisson with periodic or Neumann boundaries). The mean is removed from the right-hand side and the coarsest residuals, and the solution is returned with zero mean. Next setup.

◆ setSmoother()

template<typename T >

void amrex::AlgMG< T >::setSmoother ( Smoother sm )

inline

Takes effect at the next solve; the hierarchy is kept.

◆ setStrongThreshold()

template<typename T >

void amrex::AlgMG< T >::setStrongThreshold ( T theta )

inline

Strength threshold, typically 0.1 to 0.5 (0.25): connections at least this fraction of the strongest one in the row are strong. A row whose sum exceeds 0.9 times its diagonal has no strong connection. Takes effect at the next setup, i.e., the next solve after define.

◆ setThrowException()

template<typename T >

void amrex::AlgMG< T >::setThrowException ( bool b )

inline

◆ setTruncFactor()

template<typename T >

void amrex::AlgMG< T >::setTruncFactor ( T f )

inline

Drop interpolation weights below this fraction of the row maximum. Next setup.

◆ setup()

template<typename T >

void amrex::AlgMG< T >::setup ( )

◆ setVerbose()

template<typename T >

void amrex::AlgMG< T >::setVerbose ( int v )

inline

◆ solve()

template<typename T >

void amrex::AlgMG< T >::solve	(	AlgVector< T > &	a_sol,
		AlgVector< T > const &	a_rhs
	)

Solve A x = b. a_sol holds the initial guess on entry.

◆ strength_scale()

template<typename T >

T amrex::AlgMG< T >::strength_scale ( SpMatrix< T > const & A )

static

Largest |a_ii|, the unit of the strength matrix entries.

◆ truncate_interp()

template<typename T >

SpMatrix< T > amrex::AlgMG< T >::truncate_interp	(	SpMatrix< T > &	P,
		AlgPartition const &	cpart,
		int	max_elmts,
		T	trunc_factor
	)

static

Keep the max_elmts largest weights of each row of P (0: all) and drop weights below trunc_factor times the row maximum. Kept weights are rescaled to preserve the row's positive and negative sums.

Member Data Documentation

◆ default_p_max_elmts

template<typename T >

constexpr int amrex::AlgMG< T >::default_p_max_elmts = 4

staticconstexpr

◆ max_p_elmts

template<typename T >

constexpr int amrex::AlgMG< T >::max_p_elmts = 32

staticconstexpr

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

/home/runner/work/amrex/amrex/Src/LinearSolvers/AMReX_AlgMG.H

Public Types

Public Member Functions

Static Public Member Functions

Static Public Attributes

Detailed Description

Member Typedef Documentation

◆ BottomSolver

◆ InterpType

◆ KrylovSolver

◆ Smoother

◆ TS

Constructor & Destructor Documentation

◆ AlgMG() [1/2]

◆ AlgMG() [2/2]

Member Function Documentation

◆ build_interp()

◆ build_interp_mm()

◆ coarse_numbering()

◆ coarsen_level()

◆ create_soc()

◆ define()

◆ getInterp()

◆ getMatrix()

◆ getNumIters()

◆ getNumSetups()

◆ getResidualNorm()

◆ getSetupPeakBytes()

◆ interpolation()

◆ normalize_rows()

◆ numGlobalRows()

◆ numLevels()

◆ pmis()

◆ precond()

◆ project_out()

◆ second_pass_strength()

◆ select_coarse_rows()

◆ setAbsTol()

◆ setAggressiveDirectInterp()

◆ setAggressiveNumLevels()

◆ setBottomMaxIter()

◆ setBottomSmooth()

◆ setBottomSolver()

◆ setBottomTol()

◆ setBottomVerbose()

◆ setChebyshevDegree()

◆ setChebyshevRatio()

◆ setFixedIter()

◆ setInterpType()

◆ setKrylovSolver()

◆ setMaxCoarseSize()

◆ setMaxIter()

◆ setMaxLevels()

◆ setPMaxElmts()

◆ setPostSmooth()

◆ setPreSmooth()

◆ setPrintIndentation()

◆ setRelaxWeight()

◆ setRelTol()

◆ setSingular()

◆ setSmoother()

◆ setStrongThreshold()

◆ setThrowException()

◆ setTruncFactor()

◆ setup()

◆ setVerbose()

◆ solve()

◆ strength_scale()

◆ truncate_interp()

Member Data Documentation

◆ default_p_max_elmts

◆ max_p_elmts