LSLPB Example — Banded Positive Definite Solver

Demonstrates solving a banded symmetric positive definite linear system using linearsystems.lslpb().

Example Code

"""IMSL LSLPB example: solve a 6x6 positive definite banded system.

Creates a 6x6 tridiagonal positive definite matrix (2 on diagonal, -1 off),
solves with lslpb, and compares solution to right-hand side.

Outputs:
- Table printed to stdout
- SVG comparison plot saved to test_output/example_imsl_lslpb.svg
"""

from __future__ import annotations

from pathlib import Path
from typing import Dict

import matplotlib.pyplot as plt
import numpy as np

from linearsystems import lslpb


def run_demo_imsl_lslpb() -> Dict[str, object]:
    """Run IMSL LSLPB example: 6x6 banded positive definite system.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``x`` (np.ndarray),
            ``residual`` (float), and ``plot_path`` (str).
    """
    n = 6
    # Tridiagonal: 2 on main, -1 on sub and super diagonals
    a = 2.0 * np.eye(n) - np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)
    b = np.arange(1.0, n + 1)

    result = lslpb(a, b)
    x = result.x

    print("\nIMSL LSLPB Example: 6×6 tridiagonal positive definite system")
    print("-" * 55)
    print(f"Matrix A (tridiagonal, 2 on diagonal, -1 off-diagonal):")
    print(f"  Size: {n}×{n}")
    print(f"  Right-hand side b = {b}")
    print("-" * 55)
    print(f"{'i':<5} {'b[i]':>10} {'x[i]':>12}")
    print("-" * 55)
    for i in range(n):
        print(f"{i:<5} {b[i]:>10.4f} {x[i]:>12.6f}")
    print("-" * 55)
    print(f"{'‖Ax-b‖':<30} {result.residual:.4e}")
    print("-" * 55)

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_lslpb.svg"

    idx = np.arange(n)
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Left: b vs x comparison
    ax = axes[0]
    ax.plot(idx, b, "o-", label="b (RHS)", color="#1f77b4", markersize=8, linewidth=2)
    ax.plot(idx, x, "s--", label="x (solution)", color="#ff7f0e", markersize=8, linewidth=2)
    ax.set_xlabel("Component index")
    ax.set_ylabel("Value")
    ax.set_title("LSLPB: Solution vs Right-Hand Side")
    ax.set_xticks(idx)
    ax.legend()
    ax.grid(True, alpha=0.3)

    # Right: residual per component
    ax2 = axes[1]
    residuals_comp = np.abs(a @ x - b)
    ax2.bar(idx, residuals_comp, color="#2ca02c", alpha=0.85)
    ax2.set_xlabel("Component index")
    ax2.set_ylabel("|Ax - b|[i]")
    ax2.set_title("Component-wise Residual")
    ax2.set_xticks(idx)
    ax2.grid(True, axis="y", alpha=0.3)

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved to: {plot_path}")

    return {"x": x, "residual": result.residual, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_lslpb()

Plot Output

Solution components x[i] of the 6×6 tridiagonal positive definite system.

Console Output

IMSL LSLPB Example: 6×6 tridiagonal positive definite system
-------------------------------------------------------
Matrix A (tridiagonal, 2 on diagonal, -1 off-diagonal):
  Size: 6×6
  Right-hand side b = [1. 2. 3. 4. 5. 6.]
-------------------------------------------------------
i           b[i]         x[i]
-------------------------------------------------------
0         1.0000     8.000000
1         2.0000    15.000000
2         3.0000    20.000000
3         4.0000    22.000000
4         5.0000    20.000000
5         6.0000    13.000000
-------------------------------------------------------
‖Ax-b‖                         5.0243e-15
-------------------------------------------------------

Plot saved to: test_output\example_imsl_lslpb.svg