LSLRG Example — General Linear System Solver

Demonstrates solving a general real linear system Ax=b using linearsystems.lslrg() and comparing with LU factorization.

Example Code

"""IMSL LSLRG example: solve a 5x5 general real linear system.

Solves a 5x5 random linear system (seed=42) using both direct solve (lslrg)
and LU factorization (lufact / lusolve), then compares residuals.

Outputs:
- Table printed to stdout
- SVG bar chart of solution components saved to test_output/example_imsl_lslrg.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 lslrg, lufact, lusolve


def run_demo_imsl_lslrg() -> Dict[str, object]:
    """Run IMSL LSLRG example: solve a 5x5 general linear system.

    Args:
        None

    Returns:
        Dict[str, object]: Result dict with keys ``x`` (np.ndarray),
            ``residual_direct`` (float), ``residual_lu`` (float),
            and ``plot_path`` (str).
    """
    rng = np.random.default_rng(42)
    n = 5
    a = rng.standard_normal((n, n))
    a = a @ a.T + n * np.eye(n)          # make well-conditioned SPD
    x_true = rng.standard_normal(n)
    b = a @ x_true

    # Direct solve
    result = lslrg(a, b)
    x = result.x

    # LU solve
    lu_piv = lufact(a)
    x_lu = lusolve(lu_piv, b)
    res_lu = float(np.linalg.norm(a @ x_lu - b))

    print("\nIMSL LSLRG Example: 5×5 random linear system (seed=42)")
    print("-" * 55)
    print(f"{'Component':<12} {'x_true':>12} {'x_solved':>12} {'error':>12}")
    print("-" * 55)
    for i in range(n):
        print(f"{'x[' + str(i) + ']':<12} {x_true[i]:>12.6f} {x[i]:>12.6f} {abs(x[i]-x_true[i]):>12.2e}")
    print("-" * 55)
    print(f"{'‖Ax-b‖ direct':<30} {result.residual:.4e}")
    print(f"{'‖Ax-b‖ LU':<30} {res_lu:.4e}")
    print("-" * 55)

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

    fig, axes = plt.subplots(1, 2, figsize=(12, 5))

    # Left: bar chart of solution components
    ax = axes[0]
    idx = np.arange(n)
    ax.bar(idx - 0.2, x_true, width=0.35, label="x_true", color="#1f77b4", alpha=0.85)
    ax.bar(idx + 0.2, x, width=0.35, label="x_solved", color="#ff7f0e", alpha=0.85)
    ax.set_xticks(idx)
    ax.set_xticklabels([f"x[{i}]" for i in idx])
    ax.set_title("LSLRG: Solution Components (5×5 system)")
    ax.set_ylabel("Value")
    ax.legend()
    ax.grid(True, axis="y", alpha=0.3)

    # Right: residual comparison
    ax2 = axes[1]
    labels = ["Direct (lslrg)", "LU (lufact/lusolve)"]
    residuals = [result.residual, res_lu]
    colors = ["#2ca02c", "#d62728"]
    ax2.bar(labels, residuals, color=colors, alpha=0.85)
    ax2.set_title("Residual ‖Ax−b‖ comparison")
    ax2.set_ylabel("Residual norm")
    ax2.set_yscale("log")
    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_direct": result.residual, "residual_lu": res_lu, "plot_path": str(plot_path)}


if __name__ == "__main__":
    run_demo_imsl_lslrg()

Plot Output

Solution components x[i] of the 5×5 linear system.

Console Output

IMSL LSLRG Example: 5×5 random linear system (seed=42)
-------------------------------------------------------
Component          x_true     x_solved        error
-------------------------------------------------------
x[0]            -0.352134    -0.352134     1.11e-16
x[1]             0.532309     0.532309     1.11e-16
x[2]             0.365444     0.365444     0.00e+00
x[3]             0.412733     0.412733     1.67e-16
x[4]             0.430821     0.430821     5.55e-17
-------------------------------------------------------
‖Ax-b‖ direct                  1.3688e-15
‖Ax-b‖ LU                      1.3688e-15
-------------------------------------------------------

Plot saved to: test_output\example_imsl_lslrg.svg