0
# PyAMG: Algebraic Multigrid Solvers
1

2
PyAMG is a comprehensive library of Algebraic Multigrid (AMG) solvers with a convenient Python interface. It provides optimal or near-optimal efficiency for solving large sparse linear systems, particularly excelling with problems discretized on unstructured meshes and irregular grids where geometric information is limited or unavailable.
3

4
The library features implementations of Ruge-Stuben (Classical AMG), Smoothed Aggregation (SA), Adaptive Smoothed Aggregation (αSA), and Compatible Relaxation (CR) methods, combining Python accessibility with performance-critical C++ extensions.
5

6
## Package Information
7

8
- **Package Name**: pyamg
9
- **Language**: Python (with C++ extensions)
10
- **Installation**: `pip install pyamg`
11
- **Documentation**: http://pyamg.readthedocs.io/en/latest/
12

13
## Core Imports
14

15
```python
16
import pyamg
17
```
18

19
Common usage imports:
20

21
```python
22
import pyamg
23
import numpy as np
24
import scipy.sparse as sp
25
```
26

27
## Basic Usage
28

29
```python
30
import pyamg
31
import numpy as np
32

33
# Create a 2D Poisson problem
34
A = pyamg.gallery.poisson((100, 100), format='csr')
35

36
# Construct multigrid hierarchy using Classical AMG
37
ml = pyamg.ruge_stuben_solver(A)
38

39
# Create random right-hand side vector
40
b = np.random.rand(A.shape[0])
41

42
# Solve the linear system Ax = b
43
x = ml.solve(b, tol=1e-10)
44

45
# Check residual
46
residual = np.linalg.norm(b - A*x)
47
print(f"Residual: {residual}")
48

49
# Alternative: use high-level solve interface  
50
x = pyamg.solve(A, b)
51
```
52

53
## Architecture
54

55
PyAMG is built around the **MultilevelSolver** class, which manages a hierarchy of progressively coarser grids and operators. The architecture enables efficient solution of large linear systems through:
56

57
- **Multilevel Structure**: Automatic generation of coarse grid hierarchies
58
- **Flexible Solvers**: Classical AMG, Smoothed Aggregation, and hybrid methods
59
- **Performance Core**: C++ implementations for computationally intensive operations
60
- **Modular Design**: Separate modules for different AMG components (aggregation, relaxation, strength measures, etc.)
61

62
Key components:
63
- **MultilevelSolver**: Main solver class managing the hierarchy
64
- **Strength of Connection**: Methods to determine grid connectivity
65
- **Coarsening**: Algorithms to select coarse grid points
66
- **Interpolation**: Operators to transfer between grid levels
67
- **Relaxation**: Smoothing methods for error correction
68

69
## Capabilities
70

71
### High-Level Solving Interface
72

73
Simple interface for solving linear systems without detailed AMG knowledge. Automatically selects appropriate solver methods and parameters based on problem characteristics.
74

75
```python { .api }
76
def solve(A, b, x0=None, tol=1e-5, maxiter=400, return_solver=False,
77
          existing_solver=None, verb=True, residuals=None):
78
    """
79
    Solve the linear system Ax = b using AMG.
80
    
81
    Parameters:
82
    - A: sparse matrix, coefficient matrix
83
    - b: array, right-hand side vector
84
    - x0: array, initial guess
85
    - tol: float, convergence tolerance (default 1e-5)
86
    - maxiter: int, maximum iterations (default 400)
87
    - return_solver: bool, return solver if True
88
    - existing_solver: MultilevelSolver, reuse existing solver
89
    - verb: bool, verbose output if True
90
    - residuals: list, store residual norms
91
    
92
    Returns:
93
    array or tuple: solution vector x, optionally with solver
94
    """
95

96
def solver(A, config):
97
    """
98
    Create AMG solver for matrix A using configuration.
99
    
100
    Parameters:
101
    - A: sparse matrix, coefficient matrix
102
    - config: dict, solver configuration dictionary
103
    
104
    Returns:
105
    MultilevelSolver: configured AMG solver
106
    """
107

108
def solver_configuration(A, B=None, verb=True):
109
    """
110
    Generate solver configuration dictionary.
111
    
112
    Parameters:
113
    - A: sparse matrix, coefficient matrix  
114
    - B: array, near null-space modes (default None)
115
    - verb: bool, verbose output if True
116
    
117
    Returns:
118
    dict: solver configuration parameters
119
    """
120
```
121

122
[High-Level Interface](./high-level-interface.md)
123

124
### AMG Solver Constructors
125

126
Factory functions for creating specific types of AMG solvers with detailed control over algorithm parameters and behavior.
127

128
```python { .api }
129
def ruge_stuben_solver(A, strength=('classical', {'theta': 0.25}),
130
                      CF=('RS', {'second_pass': False}),
131
                      interpolation='classical',
132
                      presmoother=('gauss_seidel', {'sweep': 'symmetric'}),
133
                      postsmoother=('gauss_seidel', {'sweep': 'symmetric'}),
134
                      max_levels=30, max_coarse=10, keep=False, **kwargs):
135
    """
136
    Create Classical (Ruge-Stuben) AMG solver.
137
    
138
    Parameters:
139
    - A: sparse matrix, coefficient matrix
140
    - strength: str or tuple, strength of connection measure
141
    - CF: str or tuple, coarse/fine splitting method  
142
    - interpolation: str, interpolation method
143
    - presmoother: str or tuple, pre-smoothing method
144
    - postsmoother: str or tuple, post-smoothing method
145
    - max_levels: int, maximum number of levels
146
    - max_coarse: int, maximum coarse grid size
147
    - keep: bool, keep extra operators for debugging
148
    
149
    Returns:
150
    MultilevelSolver: Classical AMG solver
151
    """
152

153
def smoothed_aggregation_solver(A, B=None, BH=None, symmetry='hermitian',
154
                               strength='symmetric', aggregate='standard',
155
                               smooth=('jacobi', {'omega': 4.0/3.0}),
156
                               presmoother=('block_gauss_seidel', {'sweep': 'symmetric'}),
157
                               postsmoother=('block_gauss_seidel', {'sweep': 'symmetric'}),
158
                               improve_candidates=(('block_gauss_seidel', 
159
                                                   {'sweep': 'symmetric', 'iterations': 4}), None),
160
                               max_levels=10, max_coarse=10, diagonal_dominance=False,
161
                               keep=False, **kwargs):
162
    """
163
    Create Smoothed Aggregation AMG solver.
164
    
165
    Parameters:
166
    - A: sparse matrix, coefficient matrix
167
    - B: array, near null-space modes
168
    - BH: array, hermitian transpose of B
169
    - symmetry: str, matrix symmetry assumption
170
    - strength: str, strength of connection measure
171
    - aggregate: str, aggregation method
172
    - smooth: tuple, prolongation smoothing method
173
    - presmoother: tuple, pre-smoothing method
174
    - postsmoother: tuple, post-smoothing method
175
    - improve_candidates: tuple, candidate improvement parameters
176
    - max_levels: int, maximum number of levels
177
    - max_coarse: int, maximum coarse grid size
178
    - diagonal_dominance: bool, diagonal dominance flag
179
    - keep: bool, keep extra operators for debugging
180
    
181
    Returns:
182
    MultilevelSolver: Smoothed Aggregation AMG solver
183
    """
184

185
def air_solver(A, strength='classical', presmoother='l1_gauss_seidel',
186
               postsmoother='l1_gauss_seidel', max_levels=10, max_coarse=500,
187
               coarse_solver='pinv2', **kwargs):
188
    """
189
    Create Approximate Ideal Restriction AMG solver.
190
    
191
    Parameters:
192
    - A: sparse matrix, coefficient matrix
193
    - strength: str, strength of connection measure
194
    - presmoother: str, pre-smoothing method
195
    - postsmoother: str, post-smoothing method
196
    - max_levels: int, maximum number of levels
197
    - max_coarse: int, maximum coarse grid size  
198
    - coarse_solver: str, coarse grid solver method
199
    
200
    Returns:
201
    MultilevelSolver: AIR AMG solver
202
    """
203

204
def adaptive_sa_solver(A, initial_candidates=None, symmetry='hermitian',
205
                       num_candidates=1, candidate_iters=5, epsilon=0.1,
206
                       max_levels=10, max_coarse=10, **kwargs):
207
    """
208
    Create Adaptive Smoothed Aggregation AMG solver.
209
    
210
    Parameters:
211
    - A: sparse matrix, coefficient matrix
212
    - initial_candidates: array, initial candidate basis
213
    - symmetry: str, matrix symmetry assumption
214
    - num_candidates: int, number of candidates to generate
215
    - candidate_iters: int, smoothing iterations per candidate
216
    - epsilon: float, target convergence factor
217
    - max_levels: int, maximum number of levels
218
    - max_coarse: int, maximum coarse grid size
219
    
220
    Returns:
221
    MultilevelSolver: Adaptive SA AMG solver
222
    """
223
```
224

225
[AMG Solver Constructors](./solver-constructors.md)
226

227
### MultilevelSolver Class
228

229
The core solver class that manages AMG hierarchies and provides solve methods with extensive control over solution process.
230

231
```python { .api }
232
class MultilevelSolver:
233
    """
234
    Multilevel AMG solver managing hierarchy of operators and grids.
235
    """
236
    
237
    def solve(self, b, x0=None, tol=1e-5, maxiter=None, 
238
              cycle='V', accel=None, **kwargs):
239
        """
240
        Solve Ax = b using multigrid cycles.
241
        
242
        Parameters:
243
        - b: array, right-hand side vector
244
        - x0: array, initial guess
245
        - tol: float, convergence tolerance
246
        - maxiter: int, maximum iterations
247
        - cycle: str, multigrid cycle type ('V', 'W', 'F')
248
        - accel: str, acceleration method
249
        
250
        Returns:
251
        array: solution vector
252
        """
253
    
254
    def aspreconditioner(self, cycle='V'):
255
        """
256
        Return linear operator for use as preconditioner.
257
        
258
        Parameters:
259
        - cycle: str, multigrid cycle type
260
        
261
        Returns:
262
        LinearOperator: preconditioner operator
263
        """
264
```
265

266
[MultilevelSolver](./multilevel-solver.md)
267

268
### Test Problem Gallery
269

270
Comprehensive collection of test matrices and problems for AMG development and evaluation, including PDE discretizations and finite element problems.
271

272
```python { .api }
273
def poisson(grid, format='csr', dtype=float):
274
    """
275
    Generate discrete Poisson operator.
276
    
277
    Parameters:
278
    - grid: tuple, grid dimensions (nx,) or (nx, ny) or (nx, ny, nz)
279
    - format: str, sparse matrix format
280
    - dtype: data type
281
    
282
    Returns:
283
    sparse matrix: Poisson operator
284
    """
285

286
def linear_elasticity(grid, format='csr', dtype=float):
287
    """
288
    Generate linear elasticity operator.
289
    
290
    Parameters:
291
    - grid: tuple, grid dimensions
292
    - format: str, sparse matrix format
293
    - dtype: data type
294
    
295
    Returns:
296
    sparse matrix: elasticity operator
297
    """
298

299
def demo(**kwargs):
300
    """
301
    Run PyAMG demonstration showing solver capabilities.
302
    """
303
```
304

305
[Test Problem Gallery](./gallery.md)
306

307
### Krylov Iterative Methods
308

309
Complete suite of Krylov subspace methods for iterative solution of linear systems, optimized for use with AMG preconditioning.
310

311
```python { .api }
312
def gmres(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
313
    """
314
    Generalized Minimal Residual method.
315
    
316
    Parameters:
317
    - A: sparse matrix or LinearOperator
318
    - b: array, right-hand side
319
    - x0: array, initial guess
320
    - tol: float, convergence tolerance
321
    - maxiter: int, maximum iterations
322
    - M: preconditioner
323
    
324
    Returns:
325
    tuple: (solution, info)
326
    """
327

328
def cg(A, b, x0=None, tol=1e-5, maxiter=None, M=None, **kwargs):
329
    """
330
    Conjugate Gradient method.
331
    
332
    Parameters:
333
    - A: sparse matrix or LinearOperator
334
    - b: array, right-hand side  
335
    - x0: array, initial guess
336
    - tol: float, convergence tolerance
337
    - maxiter: int, maximum iterations
338
    - M: preconditioner
339
    
340
    Returns:
341
    tuple: (solution, info)
342
    """
343
```
344

345
[Krylov Methods](./krylov-methods.md)
346

347
### Strength of Connection
348

349
Methods for determining the strength of connections between unknowns, fundamental to AMG coarsening strategies.
350

351
```python { .api }
352
def classical_strength_of_connection(A, theta=0.25):
353
    """
354
    Classical strength of connection measure.
355
    
356
    Parameters:
357
    - A: sparse matrix, coefficient matrix
358
    - theta: float, strength threshold
359
    
360
    Returns:
361
    sparse matrix: strength matrix
362
    """
363

364
def symmetric_strength_of_connection(A, theta=0.0):
365
    """
366
    Symmetric strength of connection measure.
367
    
368
    Parameters:
369
    - A: sparse matrix, coefficient matrix
370
    - theta: float, strength threshold
371
    
372
    Returns:
373
    sparse matrix: strength matrix
374
    """
375

376
def energy_based_strength_of_connection(A, theta=0.0, **kwargs):
377
    """
378
    Energy-based strength of connection measure.
379
    
380
    Parameters:
381
    - A: sparse matrix, coefficient matrix
382
    - theta: float, strength threshold
383
    
384
    Returns:
385
    sparse matrix: strength matrix
386
    """
387
```
388

389
[Strength of Connection](./strength-of-connection.md)
390

391
### Aggregation Methods
392

393
Algorithms for grouping fine grid points into coarse grid aggregates for Smoothed Aggregation AMG.
394

395
```python { .api }
396
def standard_aggregation(C, **kwargs):
397
    """
398
    Standard aggregation algorithm.
399
    
400
    Parameters:
401
    - C: sparse matrix, strength of connection matrix
402
    
403
    Returns:
404
    array: aggregation array mapping fine to coarse
405
    """
406

407
def naive_aggregation(C, **kwargs):
408
    """
409
    Naive aggregation algorithm.
410
    
411
    Parameters:
412
    - C: sparse matrix, strength of connection matrix
413
    
414
    Returns:
415
    array: aggregation array
416
    """
417

418
def lloyd_aggregation(C, ratio=0.03, distance='unit', **kwargs):
419
    """
420
    Lloyd clustering-based aggregation.
421
    
422
    Parameters:
423
    - C: sparse matrix, strength of connection matrix
424
    - ratio: float, desired coarsening ratio
425
    - distance: str, distance measure
426
    
427
    Returns:
428
    array: aggregation array
429
    """
430
```
431

432
[Aggregation Methods](./aggregation-methods.md)
433

434
### Relaxation Methods
435

436
Smoothing and relaxation methods for error correction in multigrid cycles.
437

438
```python { .api }
439
def gauss_seidel(A, x, b, iterations=1, sweep='forward'):
440
    """
441
    Gauss-Seidel relaxation.
442
    
443
    Parameters:
444
    - A: sparse matrix, coefficient matrix
445
    - x: array, solution vector (modified in-place)
446
    - b: array, right-hand side
447
    - iterations: int, number of relaxation sweeps
448
    - sweep: str, sweep direction ('forward', 'backward', 'symmetric')
449
    """
450

451
def jacobi(A, x, b, iterations=1, omega=1.0):
452
    """
453
    Jacobi relaxation.
454
    
455
    Parameters:
456
    - A: sparse matrix, coefficient matrix
457
    - x: array, solution vector (modified in-place)
458
    - b: array, right-hand side
459
    - iterations: int, number of relaxation sweeps
460
    - omega: float, relaxation parameter
461
    """
462
```
463

464
[Relaxation Methods](./relaxation-methods.md)
465

466
### Utilities and Linear Algebra
467

468
Helper functions for matrix operations, norms, and linear algebra utilities supporting AMG algorithms.
469

470
```python { .api }
471
def norm(x, ord=None):
472
    """
473
    Compute vector or matrix norm.
474
    
475
    Parameters:
476
    - x: array-like, input vector or matrix
477
    - ord: norm type
478
    
479
    Returns:
480
    float: computed norm
481
    """
482

483
def approximate_spectral_radius(A, tol=0.1, maxiter=10):
484
    """
485
    Estimate spectral radius of matrix.
486
    
487
    Parameters:
488
    - A: sparse matrix
489
    - tol: float, tolerance for estimation
490
    - maxiter: int, maximum iterations
491
    
492
    Returns:
493
    float: estimated spectral radius
494
    """
495

496
def make_system(A, b, x0=None, formats=None):
497
    """
498
    Create consistent linear system with proper formats.
499
    
500
    Parameters:
501
    - A: matrix, coefficient matrix
502
    - b: array, right-hand side
503
    - x0: array, initial guess
504
    - formats: dict, desired matrix formats
505
    
506
    Returns:
507
    tuple: (A, b, x0) in consistent formats
508
    """
509
```
510

511
[Utilities](./utilities.md)
512

513
## Types
514

515
```python { .api }
516
# Type aliases for common PyAMG objects
517
MultilevelSolver = pyamg.multilevel.MultilevelSolver
518
```
519

520
## Error Handling
521

522
PyAMG raises standard Python exceptions:
523

524
- **ValueError**: Invalid parameters or incompatible matrix/vector dimensions
525
- **TypeError**: Incorrect input types
526
- **RuntimeError**: Convergence failures or algorithmic issues
527
- **MemoryError**: Insufficient memory for large problems
528

529
Common error patterns:
530
- Matrix format compatibility issues (use `.tocsr()` for CSR format)
531
- Singular or near-singular matrices in coarse grid solvers
532
- Memory limitations with very large problems or dense coarse operators