0
# Utilities and Linear Algebra
1

2
Helper functions for matrix operations, norms, and linear algebra utilities supporting AMG algorithms. These functions provide essential computational building blocks for PyAMG.
3

4
## Capabilities
5

6
### Linear System Utilities
7

8
Functions for creating and manipulating linear systems.
9

10
```python { .api }
11
def make_system(A, b, x0=None, formats=None):
12
    """
13
    Create consistent linear system with proper formats.
14
    
15
    Ensures that matrix and vectors are in compatible formats
16
    and have consistent dimensions for PyAMG solvers.
17
    
18
    Parameters:
19
    - A: array-like or sparse matrix, coefficient matrix
20
    - b: array-like, right-hand side vector
21
    - x0: array-like, initial guess (optional)
22
    - formats: dict, desired matrix formats:
23
        * {'A': 'csr', 'b': 'array'}: format specifications
24
        * None: use default formats (CSR for A, array for vectors)
25
    
26
    Returns:
27
    tuple: (A, b, x0) in consistent, compatible formats
28
        - A: sparse matrix in specified format (default CSR)
29
        - b: numpy array, right-hand side
30
        - x0: numpy array, initial guess (zero if not provided)
31
    
32
    Raises:
33
    ValueError: if dimensions are incompatible
34
    TypeError: if formats cannot be converted
35
    """
36

37
def upcast(*args):
38
    """
39
    Type upcasting for numerical arrays.
40
    
41
    Determines common floating point type for multiple arrays,
42
    ensuring numerical computations use adequate precision.
43
    
44
    Parameters:
45
    - *args: array-like objects to upcast
46
    
47
    Returns:
48
    numpy.dtype: common floating point type (float32, float64, 
49
        complex64, or complex128)
50
    
51
    Notes:
52
    - Promotes to higher precision when mixed types present
53
    - Handles real/complex promotion appropriately
54
    - Used internally by PyAMG for type consistency
55
    """
56
```
57

58
**Usage Examples:**
59

60
```python
61
import pyamg
62
import numpy as np
63
from scipy import sparse
64

65
# Create consistent linear system
66
A_dense = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
67
b_list = [1, 0, 1]
68
A, b, x0 = pyamg.util.make_system(A_dense, b_list)
69
print(f"A type: {type(A)}, b type: {type(b)}")
70

71
# Specify formats
72
A, b, x0 = pyamg.util.make_system(A_dense, b_list, 
73
                                 formats={'A': 'csc', 'b': 'array'})
74

75
# Type upcasting
76
x_float32 = np.array([1.0, 2.0], dtype=np.float32)
77
x_float64 = np.array([3.0, 4.0], dtype=np.float64)
78
common_type = pyamg.util.upcast(x_float32, x_float64)
79
print(f"Common type: {common_type}")
80
```
81

82
### Vector and Matrix Norms
83

84
Comprehensive norm computations for vectors and matrices.
85

86
```python { .api }
87
def norm(x, ord=None):
88
    """
89
    Compute vector or matrix norm.
90
    
91
    Unified interface for computing various norms of vectors
92
    and matrices, with automatic handling of sparse matrices.
93
    
94
    Parameters:
95
    - x: array-like or sparse matrix, input object
96
    - ord: norm type:
97
        * None: 2-norm for vectors, Frobenius norm for matrices
98
        * 'fro': Frobenius norm (matrices)
99
        * 1, 2, np.inf: vector p-norms  
100
        * -1, -2, -np.inf: negative vector p-norms
101
        * 'nuc': nuclear norm (matrices)
102
    
103
    Returns:
104
    float: computed norm value
105
    
106
    Notes:
107
    - Handles both dense and sparse matrices efficiently
108
    - Sparse matrices use optimized algorithms
109
    - Compatible with NumPy norm interface
110
    """
111

112
def infinity_norm(x):
113
    """
114
    Compute infinity norm (maximum absolute value).
115
    
116
    Specialized function for infinity norm computation,
117
    optimized for both dense and sparse arrays.
118
    
119
    Parameters:
120
    - x: array-like or sparse matrix
121
    
122
    Returns:
123
    float: infinity norm ||x||_∞ = max|x_i|
124
    """
125
```
126

127
**Usage Examples:**
128

129
```python
130
# Vector norms
131
x = np.array([3, 4, 0])
132
norm_2 = pyamg.util.norm(x)  # 2-norm
133
norm_1 = pyamg.util.norm(x, ord=1)  # 1-norm
134
norm_inf = pyamg.util.infinity_norm(x)  # infinity norm
135
print(f"2-norm: {norm_2}, 1-norm: {norm_1}, inf-norm: {norm_inf}")
136

137
# Matrix norms
138
A = pyamg.gallery.poisson((20, 20))
139
frobenius_norm = pyamg.util.norm(A, ord='fro')
140
print(f"Frobenius norm: {frobenius_norm}")
141
```
142

143
### Spectral Analysis
144

145
Functions for eigenvalue estimation and spectral properties.
146

147
```python { .api }
148
def approximate_spectral_radius(A, tol=0.1, maxiter=10):
149
    """
150
    Estimate spectral radius of matrix.
151
    
152
    Uses power iteration to estimate the largest eigenvalue
153
    magnitude (spectral radius) without computing all eigenvalues.
154
    
155
    Parameters:
156
    - A: sparse matrix, input matrix
157
    - tol: float, tolerance for power iteration convergence
158
    - maxiter: int, maximum power iterations
159
    
160
    Returns:
161
    float: estimated spectral radius ρ(A) ≈ |λ_max|
162
    
163
    Notes:
164
    - Much faster than full eigenvalue computation
165
    - Accuracy depends on eigenvalue separation
166
    - Used internally for relaxation parameter optimization
167
    """
168

169
def condest(A, t=2):
170
    """
171
    Estimate condition number of matrix.
172
    
173
    Computes 1-norm condition number estimate using
174
    block 1-norm estimation algorithm.
175
    
176
    Parameters:
177
    - A: sparse matrix, input matrix (should be square)
178
    - t: int, number of columns in estimation (higher = more accurate)
179
    
180
    Returns:
181
    float: estimated condition number κ₁(A) = ||A||₁ ||A⁻¹||₁
182
    
183
    Notes:
184
    - Much cheaper than computing exact condition number
185
    - Provides reasonable estimate for most matrices
186
    - Used for solver diagnostics and parameter selection
187
    """
188

189
def cond(A, p=None):
190
    """
191
    Compute exact condition number.
192
    
193
    Computes condition number using matrix factorization
194
    or singular value decomposition.
195
    
196
    Parameters:
197
    - A: array-like, input matrix
198
    - p: norm type (1, 2, 'fro', np.inf)
199
    
200
    Returns:
201
    float: exact condition number
202
    
203
    Notes:
204
    - Expensive for large matrices
205
    - Provides exact result unlike condest
206
    - Use condest for large sparse matrices
207
    """
208
```
209

210
**Usage Examples:**
211

212
```python
213
# Spectral radius estimation
214
A = pyamg.gallery.poisson((50, 50))
215
rho = pyamg.util.approximate_spectral_radius(A)
216
print(f"Spectral radius: {rho:.3f}")
217

218
# Condition number estimation
219
kappa_est = pyamg.util.condest(A)
220
print(f"Estimated condition number: {kappa_est:.2e}")
221

222
# Use for relaxation parameter selection
223
omega_optimal = 2.0 / (1.0 + rho)  # Optimal Jacobi parameter
224
print(f"Optimal Jacobi omega: {omega_optimal:.3f}")
225
```
226

227
### Basic Linear Algebra Operations
228

229
Fundamental linear algebra operations optimized for sparse matrices.
230

231
```python { .api }
232
def axpy(a, x, y):
233
    """
234
    AXPY operation: y = a*x + y.
235
    
236
    Efficient implementation of scaled vector addition,
237
    fundamental building block for iterative methods.
238
    
239
    Parameters:
240
    - a: scalar, scaling factor
241
    - x: array-like, input vector
242
    - y: array-like, accumulation vector (modified in-place)
243
    
244
    Returns:
245
    None (y is modified in-place)
246
    
247
    Notes:
248
    - Optimized for cache efficiency
249
    - Used extensively in Krylov methods
250
    - Handles both real and complex arithmetic
251
    """
252

253
def ishermitian(A, tol=1e-8):
254
    """
255
    Check if matrix is Hermitian (or symmetric for real matrices).
256
    
257
    Tests whether A = A^H within specified tolerance,
258
    important for selecting appropriate algorithms.
259
    
260
    Parameters:
261
    - A: sparse or dense matrix
262
    - tol: float, tolerance for symmetry test
263
    
264
    Returns:
265
    bool: True if matrix is Hermitian/symmetric
266
    
267
    Notes:
268
    - Returns True for symmetric real matrices
269
    - Handles complex matrices (Hermitian test)
270
    - Uses efficient sparse matrix operations
271
    """
272

273
def pinv_array(A, tol=None):
274
    """
275
    Pseudoinverse of dense arrays.
276
    
277
    Computes Moore-Penrose pseudoinverse using SVD,
278
    handling singular and rectangular matrices.
279
    
280
    Parameters:
281
    - A: array-like, input matrix (can be rectangular)
282
    - tol: float, tolerance for rank determination (optional)
283
    
284
    Returns:
285
    array: pseudoinverse A⁺ such that AA⁺A = A
286
    
287
    Notes:
288
    - Handles rank-deficient matrices robustly
289
    - Used for coarse grid solvers in AMG
290
    - More expensive than direct solvers for nonsingular matrices
291
    """
292
```
293

294
**Usage Examples:**
295

296
```python
297
# AXPY operation
298
x = np.array([1.0, 2.0, 3.0])
299
y = np.array([4.0, 5.0, 6.0])
300
pyamg.util.axpy(2.0, x, y)  # y = 2*x + y
301
print(f"Result: {y}")  # [6, 9, 12]
302

303
# Symmetry test
304
A_sym = pyamg.gallery.poisson((30, 30))
305
A_nonsym = A_sym + 0.1 * sparse.triu(A_sym, k=1)
306
print(f"Poisson symmetric: {pyamg.util.ishermitian(A_sym)}")
307
print(f"Modified symmetric: {pyamg.util.ishermitian(A_nonsym)}")
308

309
# Pseudoinverse for singular matrix
310
A_sing = np.array([[1, 1], [1, 1]])  # Rank 1 matrix
311
A_pinv = pyamg.util.pinv_array(A_sing)
312
print(f"Pseudoinverse: \n{A_pinv}")
313
```
314

315
### Matrix Manipulation Utilities
316

317
Functions for matrix format conversion and manipulation.
318

319
```python { .api }
320
def get_blocksize(A):
321
    """
322
    Determine block size of block-structured matrix.
323
    
324
    Analyzes matrix structure to detect consistent block pattern,
325
    important for block-based AMG methods.
326
    
327
    Parameters:
328
    - A: sparse matrix, input matrix to analyze
329
    
330
    Returns:
331
    int: detected block size (1 if no block structure found)
332
    
333
    Notes:
334
    - Heuristic detection based on nonzero pattern
335
    - Used automatically by block smoothers
336
    - Returns 1 for general unstructured matrices
337
    """
338

339
def scale_rows(A, x):
340
    """
341
    Scale matrix rows by vector elements.
342
    
343
    Multiplies each row i of matrix A by scalar x[i],
344
    creating row-scaled matrix.
345
    
346
    Parameters:
347
    - A: sparse matrix, input matrix
348
    - x: array-like, row scaling factors
349
    
350
    Returns:
351
    sparse matrix: row-scaled matrix where row i is multiplied by x[i]
352
    """
353

354
def scale_columns(A, x):
355
    """
356
    Scale matrix columns by vector elements.
357
    
358
    Multiplies each column j of matrix A by scalar x[j],
359
    creating column-scaled matrix.
360
    
361
    Parameters:
362
    - A: sparse matrix, input matrix  
363
    - x: array-like, column scaling factors
364
    
365
    Returns:
366
    sparse matrix: column-scaled matrix where column j is multiplied by x[j]
367
    """
368

369
def symmetric_rescaling(A, x):
370
    """
371
    Apply symmetric scaling to matrix.
372
    
373
    Computes X^{-1} A X where X = diag(x), preserving
374
    eigenvalue structure while improving conditioning.
375
    
376
    Parameters:
377
    - A: sparse matrix, input matrix (preferably symmetric)
378
    - x: array-like, scaling factors
379
    
380
    Returns:
381
    sparse matrix: symmetrically scaled matrix
382
    
383
    Notes:
384
    - Preserves symmetry and definiteness
385
    - Can improve condition number
386
    - Used in preconditioning strategies
387
    """
388
```
389

390
**Usage Examples:**
391

392
```python
393
# Block size detection
394
A_block = pyamg.gallery.linear_elasticity((20, 20))
395
blocksize = pyamg.util.get_blocksize(A_block)
396
print(f"Detected block size: {blocksize}")
397

398
# Matrix scaling
399
A = pyamg.gallery.poisson((10, 10))
400
row_scales = np.random.uniform(0.5, 2.0, A.shape[0])
401
A_row_scaled = pyamg.util.scale_rows(A, row_scales)
402

403
# Symmetric scaling for preconditioning
404
diag_A = np.array(A.diagonal())
405
sqrt_diag = np.sqrt(np.abs(diag_A))
406
A_scaled = pyamg.util.symmetric_rescaling(A, sqrt_diag)
407
```
408

409
### Diagnostic and Analysis Tools
410

411
Functions for solver performance analysis and debugging.
412

413
```python { .api }
414
def profile_solver(ml, b, cycles=10):
415
    """
416
    Profile solver performance characteristics.
417
    
418
    Analyzes solver performance including timing, convergence,
419
    and memory usage for given problem size.
420
    
421
    Parameters:
422
    - ml: MultilevelSolver, AMG solver to profile
423
    - b: array, right-hand side vector for testing
424
    - cycles: int, number of test cycles
425
    
426
    Returns:
427
    dict: performance profile including:
428
        * 'setup_time': solver construction time
429
        * 'solve_time': average solve time per cycle
430
        * 'convergence_factor': average convergence factor
431
        * 'memory_usage': estimated memory usage
432
    """
433

434
def print_table(data, headers=None, title=None):
435
    """
436
    Print formatted table of data.
437
    
438
    Creates nicely formatted ASCII table for displaying
439
    solver statistics, convergence history, etc.
440
    
441
    Parameters:
442
    - data: list of lists, table data (rows × columns)
443
    - headers: list, column headers (optional)
444
    - title: str, table title (optional)
445
    
446
    Returns:
447
    None (prints formatted table)
448
    """
449

450
def asfptype(x):
451
    """
452
    Convert to floating point type.
453
    
454
    Ensures input is appropriate floating point type
455
    for numerical computations, with proper precision.
456
    
457
    Parameters:
458
    - x: array-like, input to convert
459
    
460
    Returns:
461
    array: converted to appropriate floating point type
462
    
463
    Notes:
464
    - Promotes integers to float64
465
    - Preserves existing floating point types
466
    - Handles complex types appropriately
467
    """
468
```
469

470
**Usage Examples:**
471

472
```python
473
# Solver profiling
474
A = pyamg.gallery.poisson((100, 100))
475
ml = pyamg.smoothed_aggregation_solver(A)
476
b = np.random.rand(A.shape[0])
477

478
profile = pyamg.util.profile_solver(ml, b, cycles=5)
479
for key, value in profile.items():
480
    print(f"{key}: {value}")
481

482
# Formatted output table
483
convergence_data = [
484
    [1, 1.0e-1, 0.15],
485
    [2, 1.5e-3, 0.20], 
486
    [3, 2.3e-5, 0.18]
487
]
488
headers = ['Iteration', 'Residual', 'Rate']
489
pyamg.util.print_table(convergence_data, headers=headers, 
490
                       title='Convergence History')
491
```
492

493
## Usage Guidelines
494

495
### Performance Optimization
496

497
- Use `make_system()` to ensure optimal matrix formats
498
- Check matrix properties with `ishermitian()` before solver selection
499
- Profile solvers with `profile_solver()` for performance analysis
500
- Use `condest()` rather than exact condition numbers for large matrices
501

502
### Numerical Stability
503

504
- Apply `symmetric_rescaling()` for poorly conditioned symmetric problems
505
- Use `upcast()` to ensure adequate precision in mixed-type computations
506
- Monitor spectral radius with `approximate_spectral_radius()` for parameter tuning
507

508
### Common Patterns
509

510
```python
511
# Standard solver setup with utilities
512
def setup_solver(A_input, b_input):
513
    # Ensure consistent formats
514
    A, b, x0 = pyamg.util.make_system(A_input, b_input)
515
    
516
    # Check matrix properties
517
    is_symmetric = pyamg.util.ishermitian(A)
518
    condition_est = pyamg.util.condest(A)
519
    
520
    # Select solver based on properties
521
    if is_symmetric and condition_est < 1e12:
522
        ml = pyamg.smoothed_aggregation_solver(A)
523
    else:
524
        ml = pyamg.ruge_stuben_solver(A)
525
    
526
    return ml, b, x0
527

528
# Example usage
529
A = pyamg.gallery.linear_elasticity((30, 30))
530
b = np.random.rand(A.shape[0])
531
ml, b_formatted, x0 = setup_solver(A, b)
532
x = ml.solve(b_formatted, x0=x0)
533
```