sparse_systems.zig

//! Zigen Example: Sparse Systems — 1D Heat Equation
//! Solves steady-state heat equation -d²T/dx² = source(x)
//! with Dirichlet BCs using SparseLU.

const std = @import("std");
const Zigen = @import("zigen");

pub fn main() !void {
    var gpa = std.heap.GeneralPurposeAllocator(.{}){};
    defer _ = gpa.deinit();
    const allocator = gpa.allocator();

    std.debug.print(
        \\
        \\  ╔═══════════════════════════════════════════════╗
        \\  ║    Zigen — Sparse: 1D Heat Equation           ║
        \\  ╚═══════════════════════════════════════════════╝
        \\
        \\
    , .{});

    const n = 10;
    const T_left: f64 = 100.0;
    const T_right: f64 = 20.0;
    const source: f64 = 50.0;
    const dx: f64 = 1.0 / @as(f64, @floatFromInt(n + 1));

    std.debug.print("  Grid: {d} interior nodes, dx = {d:.4}\n", .{ n, dx });
    std.debug.print("  Boundary: T_left={d:.0}°C, T_right={d:.0}°C\n\n", .{ T_left, T_right });

    // ── Build tridiagonal sparse Laplacian ────────────────────────────
    std.debug.print("━━ Building Sparse Matrix ━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var K = Zigen.SparseMatrix(f64).init(allocator, n, n);
    defer K.deinit();

    for (0..n) |i| {
        try K.insert(i, i, 2.0 / (dx * dx));
        if (i > 0) try K.insert(i, i - 1, -1.0 / (dx * dx));
        if (i < n - 1) try K.insert(i, i + 1, -1.0 / (dx * dx));
    }
    try K.finalize();

    std.debug.print("  Size: {d}x{d}, nnz={d}\n\n", .{ n, n, K.nonZeros() });

    // ── Build RHS ─────────────────────────────────────────────────────
    var b = try Zigen.VectorXd.init(allocator, n);
    defer b.deinit();
    for (0..n) |i| b.data[i] = source;
    b.data[0] += T_left / (dx * dx);
    b.data[n - 1] += T_right / (dx * dx);

    // ── Solve ─────────────────────────────────────────────────────────
    std.debug.print("━━ Solving with SparseLU ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var solver = try Zigen.SparseLU(f64).compute(allocator, K);
    defer solver.deinit();

    var T = try solver.solve(b);
    defer T.deinit();

    std.debug.print("  Position     Temperature\n", .{});
    std.debug.print("  x = 0.0000   T = {d:>7.2}°C  (boundary)\n", .{T_left});
    for (0..n) |i| {
        const x = @as(f64, @floatFromInt(i + 1)) * dx;
        std.debug.print("  x = {d:.4}   T = {d:>7.2}°C\n", .{ x, T.data[i] });
    }
    std.debug.print("  x = 1.0000   T = {d:>7.2}°C  (boundary)\n\n", .{T_right});

    // ── MatrixMarket round-trip ───────────────────────────────────────
    std.debug.print("━━ MatrixMarket Round-Trip ━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const mtx_data = try Zigen.saveMatrixMarket(f64, allocator, K);
    defer allocator.free(mtx_data);

    var K2 = try Zigen.loadMatrixMarket(f64, allocator, mtx_data);
    defer K2.deinit();

    std.debug.print("  Saved {d} bytes → loaded {d}x{d} nnz={d}\n\n", .{
        mtx_data.len, K2.rows, K2.cols, K2.nonZeros(),
    });

    std.debug.print("✅ Example complete!\n", .{});
}