Chebyshev accelerated Jacobi iterative solver for linear systems. More...
#include <chebyshev.hpp>
Public Types | |
| using | OperatorType = OperatorT |
| Operator type to be solved. | |
| using | SolutionVectorType = SrcOf< OperatorType > |
| Solution vector type. | |
| using | RHSVectorType = DstOf< OperatorType > |
| Right-hand side vector type. | |
| using | ScalarType = SolutionVectorType::ScalarType |
| Scalar type for computations. | |
Public Member Functions | |
| Chebyshev (const int order, const SolutionVectorType &inverse_diagonal, const std::vector< SolutionVectorType > &tmps, const int iterations=1, const int max_ev_power_iterations=50) | |
| Construct a Chebyshev accelerated Jacobi solver. | |
| void | solve_impl (OperatorType &A, SolutionVectorType &x, const RHSVectorType &b) |
| Solve the linear system using the Chebyshev iteration. Applies the update rule for the specified number of iterations. | |
| SolutionVectorType & | get_inverse_diagonal () |
| void | refresh_max_eigenvalue_estimate_in_next_solve () |
Chebyshev accelerated Jacobi iterative solver for linear systems.
Computes (given an order \(p \geq 1, \ p \in \mathbb{N}\))
\[ x^{k+1} = x^k + \alpha_k \, D^{-1}(b - A x^k) + \beta_k \, d^k, \quad d^0 = 0 \]
for \(k = 0\) up to \(k+1 = p\), where \(D = \mathrm{diag}(A)\), and the coefficients \( \alpha_k \) and \( \beta_k \) are computed recursively from
\[ \begin{aligned} \rho_0 &= \frac{\delta}{\theta}, \quad \\ \rho_{k+1} &= (2(\theta/\delta) - \rho_k)^{-1}, \quad \\ \alpha_0 &= \theta^{-1}, \\ \alpha_k &= \frac{2 \, \rho_{k+1}}{\delta}, \quad \\ \beta_0 &= 0, \\ \beta_k &= \rho_{k+1} \rho_k, \\ d_0 &= 0, \\ d_k &= x^{k+1} - x^{k} = \alpha_k D^{-1}(b - A x^k) + \beta_k d_{k-1} \end{aligned} \]
with \( \theta = (\lambda_{\max} + \lambda_{\min})/2 \) and \( \delta = (\lambda_{\max} - \lambda_{\min})/2\) , and \( \lambda_{\min}, \lambda_{\max} \) the eigenvalue bounds of \( D^{-1}A \).
The max eigenvalue \(\lambda_{\max}\) is estimated before the first solve (or if requested via refresh_max_eigenvalue_estimate_in_next_solve()) applying a specified number of power iterations. A safety margin is added to hopefully ensure that the estimate is not too far off. The minimum eigenvalue is the obtained by scaling the maximum eigenvalue.
Each iteration executes \(p\) matrix-vector products (plus some vector-vector operations). Setting the order \(p = 1\) recovers the standard weighted Jacobi relaxation. If unsure, try low-order settings with \(p = 2\) or \(p = 3\).
The number of iterations repeats the entire process. For instance, setting the iterations \(\text{it} = 10\) will perform \(\text{it} \times p\) updates.
Satisfies the SolverLike concept (see solver.hpp).
| OperatorT | Operator type (must satisfy OperatorLike). |
| using terra::linalg::solvers::Chebyshev< OperatorT >::OperatorType = OperatorT |
Operator type to be solved.
| using terra::linalg::solvers::Chebyshev< OperatorT >::RHSVectorType = DstOf< OperatorType > |
Right-hand side vector type.
| using terra::linalg::solvers::Chebyshev< OperatorT >::ScalarType = SolutionVectorType::ScalarType |
Scalar type for computations.
| using terra::linalg::solvers::Chebyshev< OperatorT >::SolutionVectorType = SrcOf< OperatorType > |
Solution vector type.
|
inline |
Construct a Chebyshev accelerated Jacobi solver.
| order | Chebyshev order \(p >= 1\). Equivalent to weighted Jacobi for \(p = 1\). |
| inverse_diagonal | Inverse of the diagonal of the operator (D^-1). |
| tmps | Temporary vectors for workspace. We need 2 here. |
| iterations | Number of Chebyshev iterations to perform. |
| max_ev_power_iterations | Number of power iterations executed to estimate the max eigenvalue of ((D^-1)A). An additional safety margin is added internally during the setup of the Chebyshev solver. |
|
inline |
|
inline |
|
inline |
Solve the linear system using the Chebyshev iteration. Applies the update rule for the specified number of iterations.
| A | Operator (matrix). |
| x | Solution vector (output). |
| b | Right-hand side vector (input). |
1.9.8 