0
# Relaxation Methods
1

2
Smoothing and relaxation methods for error correction in multigrid cycles. These methods reduce high-frequency error components and are essential for multigrid convergence.
3

4
## Capabilities
5

6
### Point Relaxation Methods
7

8
Classical point-wise relaxation schemes for scalar problems.
9

10
```python { .api }
11
def gauss_seidel(A, x, b, iterations=1, sweep='forward'):
12
    """
13
    Gauss-Seidel relaxation.
14
    
15
    Point-wise relaxation method that updates unknowns sequentially
16
    using most recent values. Excellent smoother for many AMG problems.
17
    
18
    Parameters:
19
    - A: sparse matrix, coefficient matrix
20
    - x: array, solution vector (modified in-place)
21
    - b: array, right-hand side vector
22
    - iterations: int, number of relaxation sweeps (default 1)
23
    - sweep: str, sweep direction:
24
        * 'forward': natural ordering (i = 0, 1, ..., n-1)
25
        * 'backward': reverse ordering (i = n-1, ..., 1, 0)
26
        * 'symmetric': forward then backward sweep
27
    
28
    Returns:
29
    None (x is modified in-place)
30
    
31
    Notes:
32
    - Updates x[i] = (b[i] - sum(A[i,j]*x[j] for j != i)) / A[i,i]
33
    - Requires diagonal dominance for convergence
34
    - More effective than Jacobi for many problems
35
    """
36

37
def jacobi(A, x, b, iterations=1, omega=1.0):
38
    """
39
    Jacobi relaxation.
40
    
41
    Point-wise relaxation using diagonal scaling. Simultaneous
42
    updates make it naturally parallel but may converge slower
43
    than Gauss-Seidel.
44
    
45
    Parameters:
46
    - A: sparse matrix, coefficient matrix
47
    - x: array, solution vector (modified in-place)  
48
    - b: array, right-hand side vector
49
    - iterations: int, number of relaxation sweeps
50
    - omega: float, relaxation parameter (0 < omega <= 1):
51
        * 1.0: standard Jacobi (default)
52
        * 2/3: common damped choice
53
        * optimal: 2/(lambda_max + lambda_min)
54
    
55
    Returns:
56
    None (x is modified in-place)
57
    
58
    Notes:
59
    - Updates x = x + omega * D^{-1} * (b - A*x)
60
    - Parallel-friendly but may need damping
61
    - Good for problems with strong diagonal dominance
62
    """
63
```
64

65
**Usage Examples:**
66

67
```python
68
import pyamg
69
import numpy as np
70

71
# Setup test problem
72
A = pyamg.gallery.poisson((50, 50))
73
b = np.random.rand(A.shape[0])
74
x = np.zeros_like(b)
75

76
# Gauss-Seidel relaxation
77
x_gs = x.copy()
78
pyamg.relaxation.gauss_seidel(A, x_gs, b, iterations=5)
79
residual_gs = np.linalg.norm(b - A*x_gs)
80
print(f"Gauss-Seidel residual: {residual_gs:.2e}")
81

82
# Jacobi relaxation with damping  
83
x_jacobi = x.copy()
84
pyamg.relaxation.jacobi(A, x_jacobi, b, iterations=5, omega=2.0/3.0)
85
residual_jacobi = np.linalg.norm(b - A*x_jacobi)
86
print(f"Jacobi residual: {residual_jacobi:.2e}")
87

88
# Symmetric Gauss-Seidel
89
x_sym = x.copy()
90
pyamg.relaxation.gauss_seidel(A, x_sym, b, iterations=3, sweep='symmetric')
91
```
92

93
### Block Relaxation Methods
94

95
Block-wise relaxation methods for systems of equations and improved robustness.
96

97
```python { .api }
98
def block_gauss_seidel(A, x, b, blocksize=None, iterations=1, sweep='forward'):
99
    """
100
    Block Gauss-Seidel relaxation.
101
    
102
    Relaxation method that solves small blocks of equations exactly
103
    rather than single equations. More robust for systems problems
104
    and strongly coupled equations.
105
    
106
    Parameters:
107
    - A: sparse matrix, coefficient matrix
108
    - x: array, solution vector (modified in-place)
109
    - b: array, right-hand side vector  
110
    - blocksize: int or None, block size:
111
        * None: automatic detection from matrix structure
112
        * int: fixed block size (e.g., 2 for 2D elasticity)
113
    - iterations: int, number of block relaxation sweeps
114
    - sweep: str, sweep direction ('forward', 'backward', 'symmetric')
115
    
116
    Returns:
117
    None (x is modified in-place)
118
    
119
    Notes:
120
    - Solves A_ii * x_i = b_i - sum(A_ij * x_j) for each block i
121
    - More expensive per iteration but more effective
122
    - Essential for systems like elasticity equations
123
    """
124

125
def block_jacobi(A, x, b, blocksize=None, iterations=1, omega=1.0):
126
    """
127
    Block Jacobi relaxation.
128
    
129
    Block version of Jacobi relaxation that solves diagonal blocks
130
    exactly while treating off-diagonal blocks explicitly.
131
    
132
    Parameters:
133
    - A: sparse matrix, coefficient matrix (preferably block structured)
134
    - x: array, solution vector (modified in-place)
135
    - b: array, right-hand side vector
136
    - blocksize: int or None, block size for decomposition
137
    - iterations: int, number of block relaxation sweeps  
138
    - omega: float, block relaxation parameter
139
    
140
    Returns:
141
    None (x is modified in-place)
142
    
143
    Notes:  
144
    - More robust than point Jacobi for systems
145
    - Naturally parallel within each block solve
146
    - Good for problems with physical coupling
147
    """
148
```
149

150
**Usage Examples:**
151

152
```python
153
# Block relaxation for elasticity
154
A_elastic = pyamg.gallery.linear_elasticity((30, 30))
155
b_elastic = np.random.rand(A_elastic.shape[0])
156
x_elastic = np.zeros_like(b_elastic)
157

158
# Block Gauss-Seidel with 2x2 blocks (displacement components)
159
pyamg.relaxation.block_gauss_seidel(A_elastic, x_elastic, b_elastic, 
160
                                   blocksize=2, iterations=3)
161

162
# Block Jacobi
163
x_block_jacobi = np.zeros_like(b_elastic)
164
pyamg.relaxation.block_jacobi(A_elastic, x_block_jacobi, b_elastic,
165
                             blocksize=2, omega=0.8)
166
```
167

168
### Specialized Relaxation Methods
169

170
Advanced relaxation methods for specific problem types and performance optimization.
171

172
```python { .api }
173
def sor_gauss_seidel(A, x, b, omega=1.0, iterations=1, sweep='forward'):
174
    """
175
    SOR (Successive Over-Relaxation) Gauss-Seidel method.
176
    
177
    Accelerated version of Gauss-Seidel with over-relaxation parameter.
178
    Can improve convergence rate when properly tuned.
179
    
180
    Parameters:
181
    - A: sparse matrix, coefficient matrix
182
    - x: array, solution vector (modified in-place)
183
    - b: array, right-hand side vector
184
    - omega: float, over-relaxation parameter:
185
        * 1.0: standard Gauss-Seidel
186
        * 1 < omega < 2: over-relaxation (can accelerate)
187
        * 0 < omega < 1: under-relaxation (more stable)
188
    - iterations: int, number of SOR sweeps
189
    - sweep: str, sweep direction
190
    
191
    Returns:
192
    None (x is modified in-place)
193
    
194
    Notes:
195
    - Optimal omega depends on spectral properties of A
196
    - Can dramatically improve or worsen convergence
197
    - Requires careful parameter tuning
198
    """
199

200
def chebyshev(A, x, b, omega=1.0, degree=1):
201
    """
202
    Chebyshev polynomial relaxation.
203
    
204
    Uses Chebyshev polynomials to accelerate basic relaxation methods.
205
    Provides optimal acceleration when eigenvalue bounds are known.
206
    
207
    Parameters:
208
    - A: sparse matrix, coefficient matrix
209
    - x: array, solution vector (modified in-place)
210
    - b: array, right-hand side vector
211
    - omega: float, relaxation parameter
212
    - degree: int, Chebyshev polynomial degree
213
    
214
    Returns:
215
    None (x is modified in-place)
216
    
217
    Notes:
218
    - Requires eigenvalue estimates for optimal performance
219
    - Higher degree can improve convergence
220
    - More complex than basic methods but potentially more effective
221
    """
222
```
223

224
**Usage Examples:**
225

226
```python
227
# SOR with over-relaxation
228
A = pyamg.gallery.poisson((40, 40))
229
b = np.random.rand(A.shape[0])
230
x_sor = np.zeros_like(b)
231
pyamg.relaxation.sor_gauss_seidel(A, x_sor, b, omega=1.2, iterations=5)
232

233
# Chebyshev acceleration
234
x_cheby = np.zeros_like(b)
235
pyamg.relaxation.chebyshev(A, x_cheby, b, degree=3)
236
```
237

238
### Krylov Relaxation Methods
239

240
Krylov subspace methods used as smoothers within multigrid cycles.
241

242
```python { .api }
243
def krylov_relaxation(A, x, b, method='cg', iterations=1, **kwargs):
244
    """
245
    Krylov method as relaxation smoother.
246
    
247
    Uses truncated Krylov methods as smoothers, providing
248
    more sophisticated error reduction than classical methods.
249
    
250
    Parameters:
251
    - A: sparse matrix, coefficient matrix
252
    - x: array, solution vector (modified in-place)
253
    - b: array, right-hand side vector
254
    - method: str, Krylov method:
255
        * 'cg': Conjugate Gradient (for SPD)
256
        * 'gmres': GMRES (general matrices)
257
        * 'bicgstab': BiCGSTAB (nonsymmetric)
258
    - iterations: int, number of Krylov iterations
259
    - **kwargs: Krylov method parameters
260
    
261
    Returns:
262
    None (x is modified in-place)
263
    
264
    Notes:
265
    - More expensive per iteration than classical methods
266
    - Can be very effective for difficult problems
267
    - Requires careful iteration count selection
268
    """
269
```
270

271
**Usage Examples:**
272

273
```python
274
# CG as smoother for SPD problems
275
A_spd = pyamg.gallery.poisson((30, 30))
276
b_spd = np.random.rand(A_spd.shape[0])
277
x_cg_smooth = np.zeros_like(b_spd)
278
pyamg.relaxation.krylov_relaxation(A_spd, x_cg_smooth, b_spd, 
279
                                  method='cg', iterations=3)
280
```
281

282
## Relaxation Constants and Configurations
283

284
### Default Parameters
285

286
```python { .api }
287
# Default relaxation parameters
288
DEFAULT_SWEEP = 'forward'
289
DEFAULT_NITER = 1
290
SYMMETRIC_RELAXATION = ['jacobi', 'block_jacobi', 'chebyshev']
291
KRYLOV_RELAXATION = ['cg', 'gmres', 'bicgstab']
292
```
293

294
### Relaxation Tuples
295

296
Methods can be specified as tuples with custom parameters:
297

298
```python
299
# Custom relaxation configurations
300
presmoother = ('gauss_seidel', {'sweep': 'forward', 'iterations': 2})
301
postsmoother = ('jacobi', {'omega': 0.67, 'iterations': 1})
302
block_smoother = ('block_gauss_seidel', {'blocksize': 2, 'iterations': 1})
303

304
# Use in solver construction
305
ml = pyamg.ruge_stuben_solver(A, 
306
                             presmoother=presmoother,
307
                             postsmoother=postsmoother)
308
```
309

310
## Selection Guidelines
311

312
### Method Selection by Problem Type
313

314
- **Scalar PDEs**: Gauss-Seidel (forward/backward)
315
- **Systems (elasticity)**: Block Gauss-Seidel with appropriate block size
316
- **SPD problems**: Symmetric Gauss-Seidel or Jacobi
317
- **Indefinite problems**: Jacobi with damping
318
- **Parallel computing**: Jacobi or block Jacobi
319

320
### Performance Considerations
321

322
```python
323
# Comparison of relaxation effectiveness
324
def compare_relaxation_methods(A, b, methods):
325
    x0 = np.zeros_like(b)
326
    results = {}
327
    
328
    for name, (method, params) in methods.items():
329
        x = x0.copy()
330
        method(A, x, b, **params)
331
        residual = np.linalg.norm(b - A*x)
332
        results[name] = residual
333
    
334
    return results
335

336
# Example comparison
337
A = pyamg.gallery.poisson((40, 40))
338
b = np.random.rand(A.shape[0])
339

340
methods = {
341
    'Gauss-Seidel': (pyamg.relaxation.gauss_seidel, {'iterations': 3}),
342
    'Jacobi': (pyamg.relaxation.jacobi, {'iterations': 3, 'omega': 2/3}),
343
    'Symmetric GS': (pyamg.relaxation.gauss_seidel, 
344
                     {'iterations': 2, 'sweep': 'symmetric'})
345
}
346

347
residuals = compare_relaxation_methods(A, b, methods)
348
for method, residual in residuals.items():
349
    print(f"{method}: {residual:.2e}")
350
```
351

352
### Smoothing Strategy
353

354
- **Pre-smoothing**: Forward sweep to reduce high-frequency error
355
- **Post-smoothing**: Backward sweep for complementary error reduction  
356
- **Symmetric**: Forward + backward for symmetric problems
357
- **Multiple iterations**: 2-3 iterations often optimal balance
358

359
### Parameter Tuning
360

361
```python
362
# Conservative (stable) parameters
363
conservative_jacobi = ('jacobi', {'omega': 0.5, 'iterations': 2})
364
conservative_gs = ('gauss_seidel', {'sweep': 'symmetric', 'iterations': 1})
365

366
# Aggressive (faster but less stable) parameters  
367
aggressive_sor = ('sor_gauss_seidel', {'omega': 1.5, 'iterations': 1})
368
aggressive_jacobi = ('jacobi', {'omega': 1.0, 'iterations': 1})
369
```
370

371
### Common Issues
372

373
- **Slow convergence**: Try block methods or increase iterations
374
- **Instability**: Reduce omega parameter or use symmetric sweeps
375
- **Systems problems**: Use block methods with proper block size
376
- **Parallel efficiency**: Prefer Jacobi variants over Gauss-Seidel
377
- **Memory constraints**: Point methods use less memory than block methods