High-Level Solving Interface

def solve(A, b, x0=None, tol=1e-5, maxiter=400, return_solver=False,
          existing_solver=None, verb=True, residuals=None):
    """
    Solve the linear system Ax = b using AMG.
    
    Automatically selects optimal solver method and parameters for the given
    matrix. Uses robust smoothed aggregation settings by default.
    
    Parameters:
    - A: sparse matrix, coefficient matrix (scipy.sparse format)
    - b: array-like, right-hand side vector
    - x0: array-like, initial guess (default: random vector)
    - tol: float, convergence tolerance (default 1e-5)
    - maxiter: int, maximum iterations (default 400)
    - return_solver: bool, if True return solver along with solution
    - existing_solver: MultilevelSolver, reuse existing solver to save setup cost
    - verb: bool, print verbose output during solve
    - residuals: list, stores residual norms if provided
    
    Returns:
    array or tuple: solution vector x, or (x, ml) if return_solver=True
    
    Raises:
    ValueError: if A and b have incompatible dimensions
    RuntimeError: if solver fails to converge
    """

import pyamg
import numpy as np
from scipy.sparse import csr_matrix

# Simple usage with automatic method selection
A = pyamg.gallery.poisson((100, 100))
b = np.random.rand(A.shape[0])
x = pyamg.solve(A, b)

# Control convergence and verbosity
x = pyamg.solve(A, b, tol=1e-8, maxiter=100, verb=False)

# Return solver for reuse
x, ml = pyamg.solve(A, b, return_solver=True)

# Reuse existing solver (efficient for multiple solves)
b2 = np.random.rand(A.shape[0])
x2 = pyamg.solve(A, b2, existing_solver=ml)

# Monitor convergence with residuals
residuals = []
x = pyamg.solve(A, b, residuals=residuals)
print(f"Converged in {len(residuals)} iterations")

# Provide initial guess
x0 = np.ones(A.shape[0])
x = pyamg.solve(A, b, x0=x0)

def solver(A, config):
    """
    Create AMG solver for matrix A using provided configuration.
    
    Creates smoothed aggregation solver using the provided configuration
    dictionary. Configuration should be generated using solver_configuration().
    
    Parameters:
    - A: sparse matrix, coefficient matrix
    - config: dict, solver configuration parameters
    
    Returns:
    MultilevelSolver: configured smoothed aggregation solver
    
    Raises:
    TypeError: if solver creation fails
    ValueError: if matrix is not square or has invalid properties
    """

# Create solver with automatic configuration
A = pyamg.gallery.poisson((50, 50))
config = pyamg.solver_configuration(A)
ml = pyamg.solver(A, config)

# Solve multiple systems with same solver
b1 = np.random.rand(A.shape[0])
b2 = np.random.rand(A.shape[0])
x1 = ml.solve(b1)
x2 = ml.solve(b2)

# Generate configuration with verbose output disabled
config = pyamg.solver_configuration(A, verb=False)
ml = pyamg.solver(A, config)

# Provide near null-space for elasticity problems
A = pyamg.gallery.linear_elasticity((10, 10))
B = np.ones((A.shape[0], 1))  # Constant vector
config = pyamg.solver_configuration(A, B=B)
ml = pyamg.solver(A, config)

def solver_configuration(A, B=None, verb=True):
    """
    Generate solver configuration dictionary based on matrix analysis.
    
    Analyzes matrix properties and returns dictionary of recommended
    smoothed aggregation solver parameters.
    
    Parameters:
    - A: sparse matrix, coefficient matrix to analyze
    - B: array, near null-space modes (default None for constant vector)
    - verb: bool, verbose output during analysis (default True)
    
    Returns:
    dict: configuration dictionary with keys:
        * 'symmetry': detected matrix symmetry
        * 'strength': strength of connection parameters
        * 'aggregate': aggregation parameters
        * 'smooth': prolongation smoothing parameters
        * 'presmoother': pre-smoothing configuration
        * 'postsmoother': post-smoothing configuration
        * 'max_levels': maximum number of levels
        * 'max_coarse': maximum coarse grid size
        * 'coarse_solver': coarse grid solver method
        * 'B': near null-space modes
        * 'BH': hermitian transpose of B
        * 'keep': flag to keep operators
    
    Raises:
    ValueError: if matrix analysis fails
    """

# Get automatic configuration
A = pyamg.gallery.linear_elasticity((20, 20))
config = pyamg.solver_configuration(A)
print(config)

# Silent analysis
config = pyamg.solver_configuration(A, verb=False)

# Create solver from configuration
ml = pyamg.solver(A, config)

# Provide near null-space modes
B = np.ones((A.shape[0], 3))  # Multiple modes for elasticity
config = pyamg.solver_configuration(A, B=B, verb=False)
ml = pyamg.solver(A, config)

def make_csr(A):
    """
    Convert matrix to CSR format required by PyAMG solvers.
    
    Parameters:
    - A: array-like or sparse matrix, input matrix
    
    Returns:
    scipy.sparse.csr_matrix: matrix in CSR format
    
    Raises:
    ValueError: if conversion fails or matrix is not 2D
    """

import numpy as np
from scipy.sparse import coo_matrix

# Convert dense matrix
A_dense = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
A_csr = pyamg.make_csr(A_dense)

# Convert other sparse formats
A_coo = coo_matrix(A_dense)
A_csr = pyamg.make_csr(A_coo)

# Use in solver
ml = pyamg.solver(A_csr)