# Solvers¶

## Code Structure¶

template<class OperatorType, class VectorType, typename ValueType>
class Solver : public rocalution::RocalutionObj

Base class for all solvers and preconditioners.

Most of the solvers can be performed on linear operators LocalMatrix, LocalStencil and GlobalMatrix - i.e. the solvers can be performed locally (on a shared memory system) or in a distributed manner (on a cluster) via MPI. The only exception is the AMG (Algebraic Multigrid) solver which has two versions (one for LocalMatrix and one for GlobalMatrix class). The only pure local solvers (which do not support global/MPI operations) are the mixed-precision defect-correction solver and all direct solvers.

All solvers need three template parameters - Operators, Vectors and Scalar type.

The Solver class is purely virtual and provides an interface for

• SetOperator() to set the operator $$A$$, i.e. the user can pass the matrix here.

• Build() to build the solver (including preconditioners, sub-solvers, etc.). The user need to specify the operator first before calling Build().

• Solve() to solve the system $$Ax = b$$. The user need to pass a right-hand-side $$b$$ and a vector $$x$$, where the solution will be obtained.

• Print() to show solver information.

• ReBuildNumeric() to only re-build the solver numerically (if possible).

• MoveToHost() and MoveToAccelerator() to offload the solver (including preconditioners and sub-solvers) to the host/accelerator.

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

It provides an interface for

void rocalution::Solver::SetOperator(const OperatorType &op)

Set the Operator of the solver.

virtual void rocalution::Solver::Build(void)

Build the solver (data allocation, structure and numerical computation)

virtual void rocalution::Solver::Clear(void)

Clear (free all local data) the solver.

virtual void rocalution::Solver::Solve(const VectorType &rhs, VectorType *x) = 0

Solve Operator x = rhs.

virtual void rocalution::Solver::Print(void) const = 0

virtual void rocalution::Solver::ReBuildNumeric(void)

Rebuild the solver only with numerical computation (no allocation or data structure computation)

virtual void rocalution::Solver::MoveToHost(void)

Move all data (i.e. move the solver) to the host.

virtual void rocalution::Solver::MoveToAccelerator(void)

Move all data (i.e. move the solver) to the accelerator.

## Iterative Linear Solvers¶

template<class OperatorType, class VectorType, typename ValueType>
class IterativeLinearSolver : public rocalution::Solver<OperatorType, VectorType, ValueType>

Base class for all linear iterative solvers.

The iterative solvers are controlled by an iteration control object, which monitors the convergence properties of the solver, i.e. maximum number of iteration, relative tolerance, absolute tolerance and divergence tolerance. The iteration control can also record the residual history and store it in an ASCII file.

All iterative solvers are controlled based on

• Absolute stopping criteria, when $$|r_{k}|_{L_{p}} < \epsilon_{abs}$$

• Relative stopping criteria, when $$|r_{k}|_{L_{p}} / |r_{1}|_{L_{p}} \leq \epsilon_{rel}$$

• Divergence stopping criteria, when $$|r_{k}|_{L_{p}} / |r_{1}|_{L_{p}} \geq \epsilon_{div}$$

• Maximum number of iteration $$N$$, when $$k = N$$

where $$k$$ is the current iteration, $$r_{k}$$ the residual for the current iteration $$k$$ (i.e. $$r_{k} = b - Ax_{k}$$) and $$r_{1}$$ the starting residual (i.e. $$r_{1} = b - Ax_{init}$$). In addition, the minimum number of iterations $$M$$ can be specified. In this case, the solver will not stop to iterate, before $$k \geq M$$.

The $$L_{p}$$ norm is used for the computation, where $$p$$ could be 1, 2 and $$\infty$$. The norm computation can be set with SetResidualNorm() with 1 for $$L_{1}$$, 2 for $$L_{2}$$ and 3 for $$L_{\infty}$$. For the computation with $$L_{\infty}$$, the index of the maximum value can be obtained with GetAmaxResidualIndex(). If this function is called and $$L_{\infty}$$ was not selected, this function will return -1.

The reached criteria can be obtained with GetSolverStatus(), returning

• 0, if no criteria has been reached yet

• 1, if absolute tolerance has been reached

• 2, if relative tolerance has been reached

• 3, if divergence tolerance has been reached

• 4, if maximum number of iteration has been reached

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

It provides an interface for

void rocalution::IterativeLinearSolver::Init(double abs_tol, double rel_tol, double div_tol, int max_iter)

Initialize the solver with absolute/relative/divergence tolerance and maximum number of iterations.

void rocalution::IterativeLinearSolver::Init(double abs_tol, double rel_tol, double div_tol, int min_iter, int max_iter)

Initialize the solver with absolute/relative/divergence tolerance and minimum/maximum number of iterations.

void rocalution::IterativeLinearSolver::InitMinIter(int min_iter)

Set the minimum number of iterations.

void rocalution::IterativeLinearSolver::InitMaxIter(int max_iter)

Set the maximum number of iterations.

void rocalution::IterativeLinearSolver::InitTol(double abs, double rel, double div)

Set the absolute/relative/divergence tolerance.

void rocalution::IterativeLinearSolver::RecordResidualHistory(void)

Record the residual history.

void rocalution::IterativeLinearSolver::RecordHistory(const std::string &filename) const

Write the history to file.

virtual void rocalution::IterativeLinearSolver::Verbose(int verb = 1)

Set the solver verbosity output.

virtual void rocalution::IterativeLinearSolver::SetPreconditioner(Solver<OperatorType, VectorType, ValueType> &precond)

Set a preconditioner of the linear solver.

void rocalution::IterativeLinearSolver::SetResidualNorm(int resnorm)

Set the residual norm to $$L_1$$, $$L_2$$ or $$L_\infty$$ norm.

• resnorm = 1 -> $$L_1$$ norm

• resnorm = 2 -> $$L_2$$ norm

• resnorm = 3 -> $$L_\infty$$ norm

virtual int rocalution::IterativeLinearSolver::GetAmaxResidualIndex(void)

Return absolute maximum index of residual vector when using $$L_\infty$$ norm.

virtual int rocalution::IterativeLinearSolver::GetSolverStatus(void)

Return the current status.

## Building and Solving Phase¶

Each iterative solver consists of a building step and a solving step. During the building step all necessary auxiliary data is allocated and the preconditioner is constructed. After that, the user can call the solving procedure, the solving step can be called several times.

When the initial matrix associated with the solver is on the accelerator, the solver will try to build everything on the accelerator. However, some preconditioners and solvers (such as FSAI and AMG) need to be constructed on the host before they can be transferred to the accelerator. If the initial matrix is on the host and we want to run the solver on the accelerator then we need to move the solver to the accelerator as well as the matrix, the right-hand-side and the solution vector.

Note

If you have a preconditioner associate with the solver, it will be moved automatically to the accelerator when you move the solver.

// CG solver
CG<LocalMatrix<ValueType>, LocalVector<ValueType>, ValueType> ls;
// Multi-Colored ILU preconditioner
MultiColoredILU<LocalMatrix<ValueType>, LocalVector<ValueType>, ValueType> p;

// Move matrix and vectors to the accelerator
mat.MoveToAccelerator();
rhs.MoveToAccelerator();
x.MoveToAccelerator();

// Set mat to be the operator
ls.SetOperator(mat);
// Set p as the preconditioner of ls
ls.SetPreconditioner(p);

// Build the solver and preconditioner on the accelerator
ls.Build();

// Compute the solution on the accelerator
ls.Solve(rhs, &x);

// CG solver
CG<LocalMatrix<ValueType>, LocalVector<ValueType>, ValueType> ls;
// Multi-Colored ILU preconditioner
MultiColoredILU<LocalMatrix<ValueType>, LocalVector<ValueType>, ValueType> p;

// Set mat to be the operator
ls.SetOperator(mat);
// Set p as the preconditioner of ls
ls.SetPreconditioner(p);

// Build the solver and preconditioner on the host
ls.Build();

// Move matrix and vectors to the accelerator
mat.MoveToAccelerator();
rhs.MoveToAccelerator();
x.MoveToAccelerator();

// Move linear solver to the accelerator
ls.MoveToAccelerator();

// Compute the solution on the accelerator
ls.Solve(rhs, &x);


## Clear Function and Destructor¶

The rocalution::Solver::Clear() function clears all the data which is in the solver, including the associated preconditioner. Thus, the solver is not anymore associated with this preconditioner.

Note

The preconditioner is not deleted (via destructor), only a rocalution::Preconditioner::Clear() is called.

Note

When the destructor of the solver class is called, it automatically calls the Clear() function. Be careful, when declaring your solver and preconditioner in different places - we highly recommend to manually call the Clear() function of the solver and not to rely on the destructor of the solver.

## Numerical Update¶

Some preconditioners require two phases in the their construction: an algebraic (e.g. compute a pattern or structure) and a numerical (compute the actual values) phase. In cases, where the structure of the input matrix is a constant (e.g. Newton-like methods) it is not necessary to fully re-construct the preconditioner. In this case, the user can apply a numerical update to the current preconditioner and pass the new operator with rocalution::Solver::ReBuildNumeric(). If the preconditioner/solver does not support the numerical update, then a full rocalution::Solver::Clear() and rocalution::Solver::Build() will be performed.

## Fixed-Point Iteration¶

template<class OperatorType, class VectorType, typename ValueType>
class FixedPoint : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Fixed-Point Iteration Scheme.

The Fixed-Point iteration scheme is based on additive splitting of the matrix $$A = M + N$$. The scheme reads

$x_{k+1} = M^{-1} (b - N x_{k}).$
It can also be reformulated as a weighted defect correction scheme
$x_{k+1} = x_{k} - \omega M^{-1} (Ax_{k} - b).$
The inversion of $$M$$ can be performed by preconditioners (Jacobi, Gauss-Seidel, ILU, etc.) or by any type of solvers.

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

void rocalution::FixedPoint::SetRelaxation(ValueType omega)

Set relaxation parameter $$\omega$$.

## Krylov Subspace Solvers¶

### CG¶

template<class OperatorType, class VectorType, typename ValueType>
class CG : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

The Conjugate Gradient method is the best known iterative method for solving sparse symmetric positive definite (SPD) linear systems $$Ax=b$$. It is based on orthogonal projection onto the Krylov subspace $$\mathcal{K}_{m}(r_{0}, A)$$, where $$r_{0}$$ is the initial residual. The method can be preconditioned, where the approximation should also be SPD. SAAD

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### CR¶

template<class OperatorType, class VectorType, typename ValueType>
class CR : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Conjugate Residual Method.

The Conjugate Residual method is an iterative method for solving sparse symmetric semi-positive definite linear systems $$Ax=b$$. It is a Krylov subspace method and differs from the much more popular Conjugate Gradient method that the system matrix is not required to be positive definite. The method can be preconditioned where the approximation should also be SPD or semi-positive definite. SAAD

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### GMRES¶

template<class OperatorType, class VectorType, typename ValueType>
class GMRES : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Generalized Minimum Residual Method.

The Generalized Minimum Residual method (GMRES) is a projection method for solving sparse (non) symmetric linear systems $$Ax=b$$, based on restarting technique. The solution is approximated in a Krylov subspace $$\mathcal{K}=\mathcal{K}_{m}$$ and $$\mathcal{L}=A\mathcal{K}_{m}$$ with minimal residual, where $$\mathcal{K}_{m}$$ is the $$m$$-th Krylov subspace with $$v_{1} = r_{0}/||r_{0}||_{2}$$. SAAD

The Krylov subspace basis size can be set using SetBasisSize(). The default size is 30.

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::GMRES::SetBasisSize(int size_basis)

Set the size of the Krylov subspace basis.

### FGMRES¶

template<class OperatorType, class VectorType, typename ValueType>
class FGMRES : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Flexible Generalized Minimum Residual Method.

The Flexible Generalized Minimum Residual method (FGMRES) is a projection method for solving sparse (non) symmetric linear systems $$Ax=b$$. It is similar to the GMRES method with the only difference, the FGMRES is based on a window shifting of the Krylov subspace and thus allows the preconditioner $$M^{-1}$$ to be not a constant operator. This can be especially helpful if the operation $$M^{-1}x$$ is the result of another iterative process and not a constant operator. SAAD

The Krylov subspace basis size can be set using SetBasisSize(). The default size is 30.

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::FGMRES::SetBasisSize(int size_basis)

Set the size of the Krylov subspace basis.

### BiCGStab¶

template<class OperatorType, class VectorType, typename ValueType>
class BiCGStab : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

The Bi-Conjugate Gradient Stabilized method is a variation of CGS and solves sparse (non) symmetric linear systems $$Ax=b$$. SAAD

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### IDR¶

template<class OperatorType, class VectorType, typename ValueType>
class IDR : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Induced Dimension Reduction Method.

The Induced Dimension Reduction method is a Krylov subspace method for solving sparse (non) symmetric linear systems $$Ax=b$$. IDR(s) generates residuals in a sequence of nested subspaces. IDR1 IDR2

The dimension of the shadow space can be set by SetShadowSpace(). The default size of the shadow space is 4.

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

Set the size of the Shadow Space.

### FCG¶

template<class OperatorType, class VectorType, typename ValueType>
class FCG : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

The Flexible Conjugate Gradient method is an iterative method for solving sparse symmetric positive definite linear systems $$Ax=b$$. It is similar to the Conjugate Gradient method with the only difference, that it allows the preconditioner $$M^{-1}$$ to be not a constant operator. This can be especially helpful if the operation $$M^{-1}x$$ is the result of another iterative process and not a constant operator. fcg

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### QMRCGStab¶

template<class OperatorType, class VectorType, typename ValueType>
class QMRCGStab : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Quasi-Minimal Residual Conjugate Gradient Stabilized Method.

The Quasi-Minimal Residual Conjugate Gradient Stabilized method is a variant of the Krylov subspace BiCGStab method for solving sparse (non) symmetric linear systems $$Ax=b$$. qmrcgstab

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### BiCGStab(l)¶

template<class OperatorType, class VectorType, typename ValueType>
class BiCGStabl : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

The Bi-Conjugate Gradient Stabilized (l) method is a generalization of BiCGStab for solving sparse (non) symmetric linear systems $$Ax=b$$. It minimizes residuals over $$l$$-dimensional Krylov subspaces. The degree $$l$$ can be set with SetOrder(). bicgstabl

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::BiCGStabl::SetOrder(int l)

Set the order.

## Chebyshev Iteration Scheme¶

template<class OperatorType, class VectorType, typename ValueType>
class Chebyshev : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Chebyshev Iteration Scheme.

The Chebyshev Iteration scheme (also known as acceleration scheme) is similar to the CG method but requires minimum and maximum eigenvalues of the operator. templates

tparam OperatorType

- can be LocalMatrix, GlobalMatrix or LocalStencil

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

## Mixed-Precision Defect Correction Scheme¶

template<class OperatorTypeH, class VectorTypeH, typename ValueTypeH, class OperatorTypeL, class VectorTypeL, typename ValueTypeL>
class MixedPrecisionDC : public rocalution::IterativeLinearSolver<OperatorTypeH, VectorTypeH, ValueTypeH>

Mixed-Precision Defect Correction Scheme.

The Mixed-Precision solver is based on a defect-correction scheme. The current implementation of the library is using host based correction in double precision and accelerator computation in single precision. The solver is implemeting the scheme

$x_{k+1} = x_{k} + A^{-1} r_{k},$
where the computation of the residual $$r_{k} = b - Ax_{k}$$ and the update $$x_{k+1} = x_{k} + d_{k}$$ are performed on the host in double precision. The computation of the residual system $$Ad_{k} = r_{k}$$ is performed on the accelerator in single precision. In addition to the setup functions of the iterative solver, the user need to specify the inner ( $$Ad_{k} = r_{k}$$) solver.

tparam OperatorTypeH

- can be LocalMatrix

tparam VectorTypeH

- can be LocalVector

tparam ValueTypeH

- can be double

tparam OperatorTypeL

- can be LocalMatrix

tparam VectorTypeL

- can be LocalVector

tparam ValueTypeL

- can be float

## MultiGrid Solvers¶

The library provides algebraic multigrid as well as a skeleton for geometric multigrid methods. The BaseMultigrid class itself is not constructing the data for the method. It contains the solution procedure for V, W and K-cycles. The AMG has two different versions for Local (non-MPI) and for Global (MPI) type of computations.

template<class OperatorType, class VectorType, typename ValueType>
class BaseMultiGrid : public rocalution::IterativeLinearSolver<OperatorType, VectorType, ValueType>

Base class for all multigrid solvers Trottenberg2003.

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### Geometric MultiGrid¶

template<class OperatorType, class VectorType, typename ValueType>
class MultiGrid : public rocalution::BaseMultiGrid<OperatorType, VectorType, ValueType>

MultiGrid Method.

The MultiGrid method can be used with external data, such as externally computed restriction, prolongation and operator hierarchy. The user need to pass all this information for each level and for its construction. This includes smoothing step, prolongation/restriction, grid traversing and coarse grid solver. This data need to be passed to the solver. Trottenberg2003

• Restriction and prolongation operations can be performed in two ways, based on Restriction() and Prolongation() of the LocalVector class, or by matrix-vector multiplication. This is configured by a set function.

• Smoothers can be of any iterative linear solver. Valid options are Jacobi, Gauss-Seidel, ILU, etc. using a FixedPoint iteration scheme with pre-defined number of iterations. The smoothers could also be a solver such as CG, BiCGStab, etc.

• Coarse grid solver could be of any iterative linear solver type. The class also provides mechanisms to specify, where the coarse grid solver has to be performed, on the host or on the accelerator. The coarse grid solver can be preconditioned.

• Grid scaling based on a $$L_2$$ norm ratio.

• Operator matrices need to be passed on each grid level.

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

### Algebraic MultiGrid¶

template<class OperatorType, class VectorType, typename ValueType>
class BaseAMG : public rocalution::BaseMultiGrid<OperatorType, VectorType, ValueType>

Base class for all algebraic multigrid solvers.

The Algebraic MultiGrid solver is based on the BaseMultiGrid class. The coarsening is obtained by different aggregation techniques. The smoothers can be constructed inside or outside of the class.

All parameters in the Algebraic MultiGrid class can be set externally, including smoothers and coarse grid solver.

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

virtual void rocalution::BaseAMG::BuildHierarchy(void)

Create AMG hierarchy.

virtual void rocalution::BaseAMG::BuildSmoothers(void)

Create AMG smoothers.

void rocalution::BaseAMG::SetCoarsestLevel(int coarse_size)

Set coarsest level for hierarchy creation.

void rocalution::BaseAMG::SetManualSmoothers(bool sm_manual)

Set flag to pass smoothers manually for each level.

void rocalution::BaseAMG::SetManualSolver(bool s_manual)

Set flag to pass coarse grid solver manually.

void rocalution::BaseAMG::SetDefaultSmootherFormat(unsigned int op_format)

Set the smoother operator format.

void rocalution::BaseAMG::SetOperatorFormat(unsigned int op_format)

Set the operator format.

int rocalution::BaseAMG::GetNumLevels(void)

Returns the number of levels in hierarchy.

## Unsmoothed Aggregation AMG¶

template<class OperatorType, class VectorType, typename ValueType>
class UAAMG : public rocalution::BaseAMG<OperatorType, VectorType, ValueType>

Unsmoothed Aggregation Algebraic MultiGrid Method.

The Unsmoothed Aggregation Algebraic MultiGrid method is based on unsmoothed aggregation based interpolation scheme. stueben

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

void rocalution::UAAMG::SetCouplingStrength(ValueType eps)

Set coupling strength.

void rocalution::UAAMG::SetOverInterp(ValueType overInterp)

Set over-interpolation parameter for aggregation.

## Smoothed Aggregation AMG¶

template<class OperatorType, class VectorType, typename ValueType>
class SAAMG : public rocalution::BaseAMG<OperatorType, VectorType, ValueType>

Smoothed Aggregation Algebraic MultiGrid Method.

The Smoothed Aggregation Algebraic MultiGrid method is based on smoothed aggregation based interpolation scheme. vanek

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

void rocalution::SAAMG::SetCouplingStrength(ValueType eps)

Set coupling strength.

void rocalution::SAAMG::SetInterpRelax(ValueType relax)

Set the relaxation parameter.

## Ruge-Stueben AMG¶

template<class OperatorType, class VectorType, typename ValueType>
class RugeStuebenAMG : public rocalution::BaseAMG<OperatorType, VectorType, ValueType>

Ruge-Stueben Algebraic MultiGrid Method.

The Ruge-Stueben Algebraic MultiGrid method is based on the classic Ruge-Stueben coarsening with direct interpolation. The solver provides high-efficiency in terms of complexity of the solver (i.e. number of iterations). However, most of the time it has a higher building step and requires higher memory usage. stueben

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

void rocalution::RugeStuebenAMG::SetCouplingStrength(ValueType eps)

Set coupling strength.

## Pairwise AMG¶

template<class OperatorType, class VectorType, typename ValueType>
class PairwiseAMG : public rocalution::BaseAMG<OperatorType, VectorType, ValueType>

Pairwise Aggregation Algebraic MultiGrid Method.

The Pairwise Aggregation Algebraic MultiGrid method is based on a pairwise aggregation matching scheme. It delivers very efficient building phase which is suitable for Poisson-like equation. Most of the time it requires K-cycle for the solving phase to provide low number of iterations. This version has multi-node support. pairwiseamg

tparam OperatorType

- can be LocalMatrix or GlobalMatrix

tparam VectorType

- can be LocalVector or GlobalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

void rocalution::PairwiseAMG::SetBeta(ValueType beta)

Set beta for pairwise aggregation.

void rocalution::PairwiseAMG::SetOrdering(unsigned int ordering)

Set re-ordering for aggregation.

void rocalution::PairwiseAMG::SetCoarseningFactor(double factor)

Set target coarsening factor.

## Direct Linear Solvers¶

template<class OperatorType, class VectorType, typename ValueType>
class DirectLinearSolver : public rocalution::Solver<OperatorType, VectorType, ValueType>

Base class for all direct linear solvers.

The library provides three direct methods - LU, QR and Inversion (based on QR decomposition). The user can pass a sparse matrix, internally it will be converted to dense and then the selected method will be applied. These methods are not very optimal and due to the fact that the matrix is converted to a dense format, these methods should be used only for very small matrices.

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

template<class OperatorType, class VectorType, typename ValueType>
class LU : public rocalution::DirectLinearSolver<OperatorType, VectorType, ValueType>

LU Decomposition.

Lower-Upper Decomposition factors a given square matrix into lower and upper triangular matrix, such that $$A = LU$$.

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

template<class OperatorType, class VectorType, typename ValueType>
class QR : public rocalution::DirectLinearSolver<OperatorType, VectorType, ValueType>

QR Decomposition.

The QR Decomposition decomposes a given matrix into $$A = QR$$, such that $$Q$$ is an orthogonal matrix and $$R$$ an upper triangular matrix.

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

template<class OperatorType, class VectorType, typename ValueType>
class Inversion : public rocalution::DirectLinearSolver<OperatorType, VectorType, ValueType>

Matrix Inversion.

Full matrix inversion based on QR decomposition.

tparam OperatorType

- can be LocalMatrix

tparam VectorType

- can be LocalVector

tparam ValueType

- can be float, double, std::complex<float> or std::complex<double>

Note

These methods can only be used with local-type problems.