matrix_operations.zig

//! Zigen Example: Solving Real Linear Systems
//! LU, QR, and Cholesky for thermal equilibrium, least-squares, and SPD systems.

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

pub fn main() !void {
    std.debug.print(
        \\
        \\  ╔═══════════════════════════════════════════════╗
        \\  ║    Zigen — Solving Real Linear Systems        ║
        \\  ╚═══════════════════════════════════════════════╝
        \\
        \\
    , .{});

    // ── Problem 1: Thermal Equilibrium (3 nodes) ──────────────────────
    std.debug.print("━━ Problem 1: Thermal Equilibrium (LU) ━━━━━━━━━━━━━━━\n\n", .{});
    std.debug.print("  3-node rod: T_left=100C, T_right=25C, source=50W/m\n\n", .{});

    const K = Zigen.Matrix3f.fromArray(.{
        .{ 2.0, -1.0, 0.0 },
        .{ -1.0, 2.0, -1.0 },
        .{ 0.0, -1.0, 2.0 },
    });
    const f = Zigen.Vector3f.fromArray(.{ 150.0, 50.0, 75.0 });

    const lu = try Zigen.LU(f32, 3).compute(K);
    const T = lu.solve(f);

    std.debug.print("  T1 = {d:.1}C   T2 = {d:.1}C   T3 = {d:.1}C\n", .{
        T.at(0), T.at(1), T.at(2),
    });
    std.debug.print("  det(K) = {d:.1}\n\n", .{lu.determinant()});

    // ── Problem 2: Least-Squares Line Fit ─────────────────────────────
    std.debug.print("━━ Problem 2: Least-Squares Line Fit (QR) ━━━━━━━━━━━━\n\n", .{});
    std.debug.print("  Sensor data: x=[0..4], y=[1.1, 2.9, 5.2, 6.8, 9.1]\n\n", .{});

    const Ao = Zigen.Matrix(f32, 5, 2).fromArray(.{
        .{ 1.0, 0.0 },
        .{ 1.0, 1.0 },
        .{ 1.0, 2.0 },
        .{ 1.0, 3.0 },
        .{ 1.0, 4.0 },
    });
    const yo = Zigen.Vector(f32, 5).fromArray(.{ 1.1, 2.9, 5.2, 6.8, 9.1 });

    const qr = Zigen.QR(f32, 5, 2).compute(Ao);
    const coeff = qr.solve(yo);

    std.debug.print("  Best fit: y = {d:.3} + {d:.3}*x\n\n", .{ coeff.at(0), coeff.at(1) });

    // ── Problem 3: SPD Solve ──────────────────────────────────────────
    std.debug.print("━━ Problem 3: Cholesky Solve (SPD system) ━━━━━━━━━━━━\n\n", .{});

    const S = Zigen.Matrix3f.fromArray(.{
        .{ 6.0, 2.0, 1.0 },
        .{ 2.0, 5.0, 2.0 },
        .{ 1.0, 2.0, 4.0 },
    });
    const b = Zigen.Vector3f.fromArray(.{ 9.0, 9.0, 7.0 });

    const chol = try Zigen.Cholesky(f32, 3).compute(S);
    const x = chol.solve(b);

    std.debug.print("  x = [{d:.4}, {d:.4}, {d:.4}]\n\n", .{ x.at(0), x.at(1), x.at(2) });

    std.debug.print("Done!\n", .{});
}