Banded and Sparse Solvers

Demonstrates banded matrix storage utilities and all sparse / iterative linear-system solvers on a 10-node 1D Poisson discretisation.

Functions covered: linearsystems.to_banded(), linearsystems.from_banded(), linearsystems.lslbd(), linearsystems.sparse_solve(), linearsystems.sparse_cg(), linearsystems.sparse_gmres(), linearsystems.lsiccg(), linearsystems.lsisor().

Example Code

"""Banded and Sparse Solvers example.

Demonstrates:
- to_banded / from_banded: convert 10×10 tridiagonal (1D Poisson) to banded
  storage and back, verify round-trip
- lslbd: solve from banded storage
- sparse_solve, sparse_cg, sparse_gmres: sparse matrix solvers
- lsiccg: CG with incomplete-Cholesky preconditioner
- lsisor: SOR iteration (omega=1.2)

Problem: -u'' = f on [0,1] with Dirichlet BCs, n=10 interior points.

Outputs:
- Solution and residuals printed to stdout
- SVG: banded storage visualisation + solution comparison + convergence
  saved to test_output/example_imsl_banded_sparse.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as ssp

from linearsystems import (
    to_banded,
    from_banded,
    lslbd,
    sparse_solve,
    sparse_cg,
    sparse_gmres,
    lsiccg,
    lsisor,
)


def run_demo_imsl_banded_sparse() -> Dict[str, object]:
    """Run banded and sparse solvers demo on 1D Poisson system.

    Builds a 10×10 tridiagonal matrix for the 1D Poisson equation, exercises
    banded storage utilities and all sparse/iterative solvers, and compares
    their solutions.

    Args:
        None

    Returns:
        Dict[str, object]: Keys: ``residuals`` (dict solver→float),
            ``round_trip_error`` (float), ``plot_path`` (str).
    """
    n = 10
    h = 1.0 / (n + 1)
    x_nodes = np.linspace(h, 1.0 - h, n)

    # Exact solution: u(x) = sin(π x),  rhs: f(x) = π² sin(π x)
    u_exact = np.sin(np.pi * x_nodes)
    f = (np.pi ** 2) * np.sin(np.pi * x_nodes)
    rhs = f * h ** 2        # scaled for the (2 -1 -1)/h² system → (2 -1 -1) form
    # Build system A u = rhs where A = tridiag(-1, 2, -1) (h² absorbed into rhs)
    main = 2.0 * np.ones(n)
    off = -1.0 * np.ones(n - 1)
    A_full = np.diag(main) + np.diag(off, 1) + np.diag(off, -1)

    # ------------------------------------------------------------------
    # 1. Banded storage round-trip
    # ------------------------------------------------------------------
    lower_bw = 1
    upper_bw = 1
    ab = to_banded(A_full, lower_bw, upper_bw)
    A_reconstructed = from_banded(ab, lower_bw, upper_bw, n)
    round_trip_error = float(np.linalg.norm(A_full - A_reconstructed))

    print("\nBanded and Sparse Solvers  (1D Poisson, n=10)")
    print("=" * 60)
    print(f"\n1. Banded storage round-trip error: {round_trip_error:.4e}")
    print(f"   ab shape: {ab.shape}  (lower+upper+1=3 rows × n=10 cols)")

    # ------------------------------------------------------------------
    # 2. lslbd — solve from banded storage
    # ------------------------------------------------------------------
    result_bd = lslbd(ab, rhs, lower_bw, upper_bw)
    res_bd = result_bd.residual
    print(f"\n2. lslbd  residual: {res_bd:.4e}")

    # ------------------------------------------------------------------
    # 3. Sparse solvers
    # ------------------------------------------------------------------
    A_sp = ssp.csr_matrix(A_full)

    x_sp_direct = sparse_solve(A_sp, rhs)
    res_sp_direct = float(np.linalg.norm(A_full @ x_sp_direct - rhs))

    result_sp_cg = sparse_cg(A_sp, rhs)
    res_sp_cg = result_sp_cg.residual

    result_sp_gmres = sparse_gmres(A_sp, rhs)
    res_sp_gmres = result_sp_gmres.residual

    print(f"\n3. sparse_solve  residual: {res_sp_direct:.4e}")
    print(f"   sparse_cg    residual: {res_sp_cg:.4e}  "
          f"(iters={result_sp_cg.n_iter})")
    print(f"   sparse_gmres residual: {res_sp_gmres:.4e}  "
          f"(iters={result_sp_gmres.n_iter})")

    # ------------------------------------------------------------------
    # 4. lsiccg — IC-preconditioned CG
    # ------------------------------------------------------------------
    result_iccg = lsiccg(A_full, rhs)
    res_iccg = result_iccg.residual
    print(f"\n4. lsiccg  residual: {res_iccg:.4e}  "
          f"(iters={result_iccg.n_iter})")

    # ------------------------------------------------------------------
    # 5. lsisor — SOR iteration
    # ------------------------------------------------------------------
    result_sor = lsisor(A_full, rhs, omega=1.2, tol=1e-10, maxiter=2000)
    res_sor = result_sor.residual
    print(f"\n5. lsisor  residual: {res_sor:.4e}  "
          f"(iters={result_sor.n_iter}, success={result_sor.success})")

    residuals = {
        "lslbd": res_bd,
        "sparse_solve": res_sp_direct,
        "sparse_cg": res_sp_cg,
        "sparse_gmres": res_sp_gmres,
        "lsiccg": res_iccg,
        "lsisor": res_sor,
    }

    # ------------------------------------------------------------------
    # Plots
    # ------------------------------------------------------------------
    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_banded_sparse.svg"

    fig, axes = plt.subplots(1, 3, figsize=(18, 6))

    # (0) Banded storage visualisation
    ax = axes[0]
    im = ax.imshow(ab, cmap="viridis", aspect="auto")
    ax.set_title("Banded storage ab\n(3 rows × 10 cols)")
    ax.set_xlabel("Column j")
    ax.set_ylabel("Band row")
    ax.set_yticks([0, 1, 2])
    ax.set_yticklabels(["super-diag", "main diag", "sub-diag"])
    fig.colorbar(im, ax=ax, shrink=0.7)

    # (1) Solution comparison
    ax2 = axes[1]
    ax2.plot(x_nodes, u_exact, "k-", linewidth=2.5, label="Exact sin(πx)")
    ax2.plot(x_nodes, result_bd.x, "b--o", ms=5, linewidth=1.5, label="lslbd")
    ax2.plot(x_nodes, x_sp_direct, "r:s", ms=4, linewidth=1.5, label="sparse_solve")
    ax2.plot(x_nodes, result_sp_cg.x, "g-.^", ms=4, linewidth=1.2, label="sparse_cg")
    ax2.set_xlabel("x")
    ax2.set_ylabel("u(x)")
    ax2.set_title("Solution comparison (1D Poisson, n=10)")
    ax2.legend(fontsize=8)
    ax2.grid(True, alpha=0.3)

    # (2) Residuals bar chart
    ax3 = axes[2]
    labels = list(residuals.keys())
    vals = [max(v, 1e-18) for v in residuals.values()]
    colors = ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd", "#8c564b"]
    ax3.bar(range(len(labels)), vals, color=colors, alpha=0.85)
    ax3.set_xticks(range(len(labels)))
    ax3.set_xticklabels(labels, rotation=30, ha="right", fontsize=9)
    ax3.set_yscale("log")
    ax3.set_title("Solver residuals ‖Ax−b‖₂")
    ax3.set_ylabel("Residual (log scale)")
    ax3.grid(True, axis="y", alpha=0.3)

    fig.suptitle("Banded and Sparse Solvers — 1D Poisson (n=10)", fontsize=13)
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved to: {plot_path}")

    return {
        "residuals": residuals,
        "round_trip_error": round_trip_error,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_imsl_banded_sparse()

Plot Output

Left: banded storage visualisation (3 rows × 10 cols). Centre: solution comparison across all solvers vs exact sin(πx). Right: residual norms for all six solver methods.

Console Output

Banded and Sparse Solvers  (1D Poisson, n=10)
============================================================

1. Banded storage round-trip error: 0.0000e+00
   ab shape: (3, 10)  (lower+upper+1=3 rows × n=10 cols)

2. lslbd  residual: 5.4949e-16

3. sparse_solve  residual: 5.1174e-16
   sparse_cg    residual: 1.0122e-15  (iters=1)
   sparse_gmres residual: 1.4513e-15  (iters=1)

4. lsiccg  residual: 5.9525e-16  (iters=1)

5. lsisor  residual: 8.8765e-11  (iters=170, success=True)

Plot saved to: test_output\example_imsl_banded_sparse.svg