MultilevelSolver Class

class MultilevelSolver:
    """
    Multilevel AMG solver managing hierarchy of operators and grids.
    
    Manages the complete AMG hierarchy including:
    - Level operators (matrices)  
    - Prolongation and restriction operators
    - Pre and post smoothers
    - Coarse grid solver
    
    Provides methods for solving linear systems using multigrid cycles
    and accessing hierarchy information.
    
    Attributes:
    - levels: list of level dictionaries containing operators
    - A: list of matrices at each level
    - P: list of prolongation operators  
    - R: list of restriction operators
    - presmoother: list of pre-smoothing operators
    - postsmoother: list of post-smoothing operators
    """
    
    def solve(self, b, x0=None, tol=1e-5, maxiter=None, 
              cycle='V', accel=None, callback=None, residuals=None, **kwargs):
        """
        Solve Ax = b using multigrid cycles.
        
        Main solution method supporting various multigrid cycles,
        convergence criteria, and Krylov acceleration.
        
        Parameters:
        - b: array-like, right-hand side vector
        - x0: array-like, initial guess (default: zero vector)
        - tol: float, convergence tolerance (default 1e-5)
        - maxiter: int, maximum iterations (default: automatic)
        - cycle: str, multigrid cycle type:
            * 'V': V-cycle (default, most efficient)
            * 'W': W-cycle (more robust, higher cost)  
            * 'F': F-cycle (flexible cycle)
        - accel: str, Krylov acceleration method:
            * None: no acceleration (default)
            * 'cg': conjugate gradient acceleration
            * 'gmres': GMRES acceleration
            * 'bicgstab': BiCGSTAB acceleration
            * 'fgmres': flexible GMRES
        - callback: callable, function called after each iteration
        - residuals: list, stores residual norms if provided
        - return_residuals: bool, return residual history
        
        Returns:
        array: solution vector x such that ||b - Ax|| < tol * ||b||
        
        Raises:
        ValueError: if b has wrong dimensions or parameters are invalid
        RuntimeError: if solver fails to converge within maxiter
        """
    
    def aspreconditioner(self, cycle='V'):
        """
        Return linear operator for use as preconditioner.
        
        Creates scipy LinearOperator that applies one multigrid
        cycle, suitable for use with Krylov methods.
        
        Parameters:
        - cycle: str, multigrid cycle type ('V', 'W', 'F')
        
        Returns:
        scipy.sparse.linalg.LinearOperator: preconditioner operator
            with matvec method applying multigrid cycle
        
        Raises:
        ValueError: if cycle type is invalid
        """
    
    def __repr__(self):
        """
        String representation showing hierarchy information.
        
        Returns:
        str: formatted hierarchy summary including:
            - Number of levels
            - Grid complexity (total unknowns / fine grid unknowns)  
            - Operator complexity (total nonzeros / fine grid nonzeros)
            - Coarse solver information
            - Per-level statistics (unknowns, nonzeros, percentages)
        """
    
    def __len__(self):
        """
        Number of levels in the hierarchy.
        
        Returns:
        int: number of levels (fine grid is level 0)
        """

import pyamg
import numpy as np
from scipy.sparse.linalg import gmres

# Create solver and examine hierarchy
A = pyamg.gallery.poisson((100, 100))
ml = pyamg.ruge_stuben_solver(A)
print(ml)  # Display hierarchy information
print(f"Number of levels: {len(ml)}")

# Basic solve with V-cycle
b = np.random.rand(A.shape[0])
x = ml.solve(b)

# Solve with W-cycle and custom tolerance
x = ml.solve(b, cycle='W', tol=1e-8)

# Solve with initial guess
x0 = np.ones(A.shape[0])
x = ml.solve(b, x0=x0, maxiter=50)

# Use as preconditioner with Krylov method
M = ml.aspreconditioner()
x, info = gmres(A, b, M=M, tol=1e-8)

# Monitor convergence
residuals = []
x = ml.solve(b, residuals=residuals)
print(f"Convergence in {len(residuals)} iterations")

# Custom callback function
def callback(x):
    r = np.linalg.norm(b - A*x)
    print(f"Residual: {r:.2e}")

x = ml.solve(b, callback=callback)

class MultilevelSolver:
    
    @property 
    def A(self):
        """
        List of matrices at each level.
        
        Returns:
        list: sparse matrices [A0, A1, ..., A_{L-1}]
            where A0 is the fine grid matrix
        """
    
    @property
    def P(self):
        """  
        List of prolongation operators.
        
        Returns:
        list: prolongation matrices [P0, P1, ..., P_{L-2}]
            where Pi maps from level i+1 to level i
        """
    
    @property
    def R(self):
        """
        List of restriction operators.
        
        Returns: 
        list: restriction matrices [R0, R1, ..., R_{L-2}]
            where Ri maps from level i to level i+1
        """
    
    @property
    def levels(self):
        """
        Complete level information.
        
        Returns:
        list: level dictionaries containing:
            - 'A': operator matrix
            - 'P': prolongation operator (if not coarsest)
            - 'R': restriction operator (if not coarsest) 
            - 'presmoother': pre-smoothing operator
            - 'postsmoother': post-smoothing operator
            - Additional level-specific data
        """
        
    def grid_complexity(self):
        """
        Compute grid complexity of hierarchy.
        
        Grid complexity measures memory overhead as ratio of
        total unknowns across all levels to fine grid unknowns.
        
        Returns:
        float: grid complexity = sum(Ni) / N0
            where Ni is number of unknowns at level i
        """
    
    def operator_complexity(self):
        """
        Compute operator complexity of hierarchy.
        
        Operator complexity measures storage overhead as ratio of
        total nonzeros across all levels to fine grid nonzeros.
        
        Returns:
        float: operator complexity = sum(nnz(Ai)) / nnz(A0)
        """

# Examine hierarchy structure
A = pyamg.gallery.linear_elasticity((50, 50))
ml = pyamg.smoothed_aggregation_solver(A)

print(f"Grid complexity: {ml.grid_complexity():.2f}")
print(f"Operator complexity: {ml.operator_complexity():.2f}")

# Access level operators
print(f"Fine grid size: {ml.A[0].shape}")
print(f"Coarse grid size: {ml.A[-1].shape}")

# Examine prolongation operators
for i, P in enumerate(ml.P):
    print(f"P{i} maps {P.shape[1]} -> {P.shape[0]} points")

# Access complete level information
for i, level in enumerate(ml.levels):
    A_level = level['A']
    print(f"Level {i}: {A_level.shape[0]} unknowns, "
          f"{A_level.nnz} nonzeros")

class MultilevelSolver:
    
    def cycle(self, level, x, b, cycle='V'):
        """
        Apply one multigrid cycle starting at given level.
        
        Low-level method for applying multigrid cycles with
        explicit level control. Used internally by solve().
        
        Parameters:
        - level: int, starting level (0 = fine grid)
        - x: array, solution vector (modified in-place)
        - b: array, right-hand side vector  
        - cycle: str, cycle type ('V', 'W', 'F')
        
        Returns:
        None (x is modified in-place)
        """
        
    def _solve_level0(self, x, b, **kwargs):
        """
        Apply smoothing at finest level.
        
        Internal method applying pre/post smoothing at level 0.
        
        Parameters:
        - x: array, solution vector (modified in-place)
        - b: array, right-hand side
        - **kwargs: smoothing parameters
        """

def multilevel_solver(A, **kwargs):
    """
    Alias for MultilevelSolver constructor.
    
    Provided for backward compatibility and alternative naming.
    
    Parameters:
    - A: sparse matrix, coefficient matrix
    - **kwargs: solver construction parameters
    
    Returns:
    MultilevelSolver: AMG solver instance
    """