Utilities and Linear Algebra

def make_system(A, b, x0=None, formats=None):
    """
    Create consistent linear system with proper formats.
    
    Ensures that matrix and vectors are in compatible formats
    and have consistent dimensions for PyAMG solvers.
    
    Parameters:
    - A: array-like or sparse matrix, coefficient matrix
    - b: array-like, right-hand side vector
    - x0: array-like, initial guess (optional)
    - formats: dict, desired matrix formats:
        * {'A': 'csr', 'b': 'array'}: format specifications
        * None: use default formats (CSR for A, array for vectors)
    
    Returns:
    tuple: (A, b, x0) in consistent, compatible formats
        - A: sparse matrix in specified format (default CSR)
        - b: numpy array, right-hand side
        - x0: numpy array, initial guess (zero if not provided)
    
    Raises:
    ValueError: if dimensions are incompatible
    TypeError: if formats cannot be converted
    """

def upcast(*args):
    """
    Type upcasting for numerical arrays.
    
    Determines common floating point type for multiple arrays,
    ensuring numerical computations use adequate precision.
    
    Parameters:
    - *args: array-like objects to upcast
    
    Returns:
    numpy.dtype: common floating point type (float32, float64, 
        complex64, or complex128)
    
    Notes:
    - Promotes to higher precision when mixed types present
    - Handles real/complex promotion appropriately
    - Used internally by PyAMG for type consistency
    """

import pyamg
import numpy as np
from scipy import sparse

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

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

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

def norm(x, ord=None):
    """
    Compute vector or matrix norm.
    
    Unified interface for computing various norms of vectors
    and matrices, with automatic handling of sparse matrices.
    
    Parameters:
    - x: array-like or sparse matrix, input object
    - ord: norm type:
        * None: 2-norm for vectors, Frobenius norm for matrices
        * 'fro': Frobenius norm (matrices)
        * 1, 2, np.inf: vector p-norms  
        * -1, -2, -np.inf: negative vector p-norms
        * 'nuc': nuclear norm (matrices)
    
    Returns:
    float: computed norm value
    
    Notes:
    - Handles both dense and sparse matrices efficiently
    - Sparse matrices use optimized algorithms
    - Compatible with NumPy norm interface
    """

def infinity_norm(x):
    """
    Compute infinity norm (maximum absolute value).
    
    Specialized function for infinity norm computation,
    optimized for both dense and sparse arrays.
    
    Parameters:
    - x: array-like or sparse matrix
    
    Returns:
    float: infinity norm ||x||_∞ = max|x_i|
    """

# Vector norms
x = np.array([3, 4, 0])
norm_2 = pyamg.util.norm(x)  # 2-norm
norm_1 = pyamg.util.norm(x, ord=1)  # 1-norm
norm_inf = pyamg.util.infinity_norm(x)  # infinity norm
print(f"2-norm: {norm_2}, 1-norm: {norm_1}, inf-norm: {norm_inf}")

# Matrix norms
A = pyamg.gallery.poisson((20, 20))
frobenius_norm = pyamg.util.norm(A, ord='fro')
print(f"Frobenius norm: {frobenius_norm}")

def approximate_spectral_radius(A, tol=0.1, maxiter=10):
    """
    Estimate spectral radius of matrix.
    
    Uses power iteration to estimate the largest eigenvalue
    magnitude (spectral radius) without computing all eigenvalues.
    
    Parameters:
    - A: sparse matrix, input matrix
    - tol: float, tolerance for power iteration convergence
    - maxiter: int, maximum power iterations
    
    Returns:
    float: estimated spectral radius ρ(A) ≈ |λ_max|
    
    Notes:
    - Much faster than full eigenvalue computation
    - Accuracy depends on eigenvalue separation
    - Used internally for relaxation parameter optimization
    """

def condest(A, t=2):
    """
    Estimate condition number of matrix.
    
    Computes 1-norm condition number estimate using
    block 1-norm estimation algorithm.
    
    Parameters:
    - A: sparse matrix, input matrix (should be square)
    - t: int, number of columns in estimation (higher = more accurate)
    
    Returns:
    float: estimated condition number κ₁(A) = ||A||₁ ||A⁻¹||₁
    
    Notes:
    - Much cheaper than computing exact condition number
    - Provides reasonable estimate for most matrices
    - Used for solver diagnostics and parameter selection
    """

def cond(A, p=None):
    """
    Compute exact condition number.
    
    Computes condition number using matrix factorization
    or singular value decomposition.
    
    Parameters:
    - A: array-like, input matrix
    - p: norm type (1, 2, 'fro', np.inf)
    
    Returns:
    float: exact condition number
    
    Notes:
    - Expensive for large matrices
    - Provides exact result unlike condest
    - Use condest for large sparse matrices
    """

# Spectral radius estimation
A = pyamg.gallery.poisson((50, 50))
rho = pyamg.util.approximate_spectral_radius(A)
print(f"Spectral radius: {rho:.3f}")

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

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

def axpy(a, x, y):
    """
    AXPY operation: y = a*x + y.
    
    Efficient implementation of scaled vector addition,
    fundamental building block for iterative methods.
    
    Parameters:
    - a: scalar, scaling factor
    - x: array-like, input vector
    - y: array-like, accumulation vector (modified in-place)
    
    Returns:
    None (y is modified in-place)
    
    Notes:
    - Optimized for cache efficiency
    - Used extensively in Krylov methods
    - Handles both real and complex arithmetic
    """

def ishermitian(A, tol=1e-8):
    """
    Check if matrix is Hermitian (or symmetric for real matrices).
    
    Tests whether A = A^H within specified tolerance,
    important for selecting appropriate algorithms.
    
    Parameters:
    - A: sparse or dense matrix
    - tol: float, tolerance for symmetry test
    
    Returns:
    bool: True if matrix is Hermitian/symmetric
    
    Notes:
    - Returns True for symmetric real matrices
    - Handles complex matrices (Hermitian test)
    - Uses efficient sparse matrix operations
    """

def pinv_array(A, tol=None):
    """
    Pseudoinverse of dense arrays.
    
    Computes Moore-Penrose pseudoinverse using SVD,
    handling singular and rectangular matrices.
    
    Parameters:
    - A: array-like, input matrix (can be rectangular)
    - tol: float, tolerance for rank determination (optional)
    
    Returns:
    array: pseudoinverse A⁺ such that AA⁺A = A
    
    Notes:
    - Handles rank-deficient matrices robustly
    - Used for coarse grid solvers in AMG
    - More expensive than direct solvers for nonsingular matrices
    """

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

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

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

def get_blocksize(A):
    """
    Determine block size of block-structured matrix.
    
    Analyzes matrix structure to detect consistent block pattern,
    important for block-based AMG methods.
    
    Parameters:
    - A: sparse matrix, input matrix to analyze
    
    Returns:
    int: detected block size (1 if no block structure found)
    
    Notes:
    - Heuristic detection based on nonzero pattern
    - Used automatically by block smoothers
    - Returns 1 for general unstructured matrices
    """

def scale_rows(A, x):
    """
    Scale matrix rows by vector elements.
    
    Multiplies each row i of matrix A by scalar x[i],
    creating row-scaled matrix.
    
    Parameters:
    - A: sparse matrix, input matrix
    - x: array-like, row scaling factors
    
    Returns:
    sparse matrix: row-scaled matrix where row i is multiplied by x[i]
    """

def scale_columns(A, x):
    """
    Scale matrix columns by vector elements.
    
    Multiplies each column j of matrix A by scalar x[j],
    creating column-scaled matrix.
    
    Parameters:
    - A: sparse matrix, input matrix  
    - x: array-like, column scaling factors
    
    Returns:
    sparse matrix: column-scaled matrix where column j is multiplied by x[j]
    """

def symmetric_rescaling(A, x):
    """
    Apply symmetric scaling to matrix.
    
    Computes X^{-1} A X where X = diag(x), preserving
    eigenvalue structure while improving conditioning.
    
    Parameters:
    - A: sparse matrix, input matrix (preferably symmetric)
    - x: array-like, scaling factors
    
    Returns:
    sparse matrix: symmetrically scaled matrix
    
    Notes:
    - Preserves symmetry and definiteness
    - Can improve condition number
    - Used in preconditioning strategies
    """

# Block size detection
A_block = pyamg.gallery.linear_elasticity((20, 20))
blocksize = pyamg.util.get_blocksize(A_block)
print(f"Detected block size: {blocksize}")

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

# Symmetric scaling for preconditioning
diag_A = np.array(A.diagonal())
sqrt_diag = np.sqrt(np.abs(diag_A))
A_scaled = pyamg.util.symmetric_rescaling(A, sqrt_diag)

def profile_solver(ml, b, cycles=10):
    """
    Profile solver performance characteristics.
    
    Analyzes solver performance including timing, convergence,
    and memory usage for given problem size.
    
    Parameters:
    - ml: MultilevelSolver, AMG solver to profile
    - b: array, right-hand side vector for testing
    - cycles: int, number of test cycles
    
    Returns:
    dict: performance profile including:
        * 'setup_time': solver construction time
        * 'solve_time': average solve time per cycle
        * 'convergence_factor': average convergence factor
        * 'memory_usage': estimated memory usage
    """

def print_table(data, headers=None, title=None):
    """
    Print formatted table of data.
    
    Creates nicely formatted ASCII table for displaying
    solver statistics, convergence history, etc.
    
    Parameters:
    - data: list of lists, table data (rows × columns)
    - headers: list, column headers (optional)
    - title: str, table title (optional)
    
    Returns:
    None (prints formatted table)
    """

def asfptype(x):
    """
    Convert to floating point type.
    
    Ensures input is appropriate floating point type
    for numerical computations, with proper precision.
    
    Parameters:
    - x: array-like, input to convert
    
    Returns:
    array: converted to appropriate floating point type
    
    Notes:
    - Promotes integers to float64
    - Preserves existing floating point types
    - Handles complex types appropriately
    """

# Solver profiling
A = pyamg.gallery.poisson((100, 100))
ml = pyamg.smoothed_aggregation_solver(A)
b = np.random.rand(A.shape[0])

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

# Formatted output table
convergence_data = [
    [1, 1.0e-1, 0.15],
    [2, 1.5e-3, 0.20], 
    [3, 2.3e-5, 0.18]
]
headers = ['Iteration', 'Residual', 'Rate']
pyamg.util.print_table(convergence_data, headers=headers, 
                       title='Convergence History')

# Standard solver setup with utilities
def setup_solver(A_input, b_input):
    # Ensure consistent formats
    A, b, x0 = pyamg.util.make_system(A_input, b_input)
    
    # Check matrix properties
    is_symmetric = pyamg.util.ishermitian(A)
    condition_est = pyamg.util.condest(A)
    
    # Select solver based on properties
    if is_symmetric and condition_est < 1e12:
        ml = pyamg.smoothed_aggregation_solver(A)
    else:
        ml = pyamg.ruge_stuben_solver(A)
    
    return ml, b, x0

# Example usage
A = pyamg.gallery.linear_elasticity((30, 30))
b = np.random.rand(A.shape[0])
ml, b_formatted, x0 = setup_solver(A, b)
x = ml.solve(b_formatted, x0=x0)