dynamic_decompositions.zig

//! Dynamic Decompositions Example
//! Demonstrates the workspace-reuse pattern: allocate once, compute many.

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

pub fn main() !void {
    const allocator = std.heap.page_allocator;

    std.debug.print("=== Dynamic Decompositions: Workspace-Reuse Pattern ===\n\n", .{});

    // ── LU Decomposition ─────────────────────────────────────────────────
    {
        std.debug.print("── LU Decomposition ──\n", .{});
        const LU = Zigen.LUDynamic(f64);

        // Allocate workspace once
        var lu = try LU.init(allocator, 3);
        defer lu.deinit();

        // Matrix 1: [2 1 0; 1 3 1; 0 1 2]
        var A1 = try Zigen.MatrixXd.init(allocator, 3, 3);
        defer A1.deinit();
        const v1 = [_]f64{ 2, 1, 0, 1, 3, 1, 0, 1, 2 };
        for (0..9) |i| A1.data[i] = v1[i];
        lu.computeFrom(A1);
        std.debug.print("  Matrix A1 det = {d:.4}\n", .{lu.determinant()});

        // Matrix 2, reusing the same workspace (zero extra allocations!)
        var A2 = try Zigen.MatrixXd.init(allocator, 3, 3);
        defer A2.deinit();
        const v2 = [_]f64{ 1, 2, 3, 0, 4, 5, 1, 0, 6 };
        for (0..9) |i| A2.data[i] = v2[i];
        lu.computeFrom(A2);
        std.debug.print("  Matrix A2 det = {d:.4}\n\n", .{lu.determinant()});
    }

    // ── Cholesky Decomposition ───────────────────────────────────────────
    {
        std.debug.print("── Cholesky Decomposition ──\n", .{});
        const Chol = Zigen.CholeskyDynamic(f64);

        var chol = try Chol.init(allocator, 3);
        defer chol.deinit();

        // SPD matrix: [4 2 1; 2 5 3; 1 3 6]
        var A = try Zigen.MatrixXd.init(allocator, 3, 3);
        defer A.deinit();
        const v = [_]f64{ 4, 2, 1, 2, 5, 3, 1, 3, 6 };
        for (0..9) |i| A.data[i] = v[i];
        try chol.computeFrom(A);

        std.debug.print("  L diagonal: ", .{});
        for (0..3) |i| std.debug.print("{d:.4} ", .{chol.L.at(i, i)});
        std.debug.print("\n\n", .{});
    }

    // ── QR Decomposition ─────────────────────────────────────────────────
    {
        std.debug.print("── QR Decomposition ──\n", .{});
        const QR = Zigen.QRDynamic(f64);

        var qr = try QR.init(allocator, 3, 3);
        defer qr.deinit();

        var A = try Zigen.MatrixXd.init(allocator, 3, 3);
        defer A.deinit();
        const v = [_]f64{ 1, 2, 3, 4, 5, 6, 7, 8, 0 };
        for (0..9) |i| A.data[i] = v[i];
        qr.computeFrom(A);

        std.debug.print("  R diagonal: ", .{});
        for (0..3) |i| std.debug.print("{d:.4} ", .{qr.R.at(i, i)});
        std.debug.print("\n\n", .{});
    }

    // ── SVD ──────────────────────────────────────────────────────────────
    {
        std.debug.print("── SVD ──\n", .{});
        const SVD = Zigen.SVDDynamic(f64);

        var svd = try SVD.init(allocator, 3, 3);
        defer svd.deinit();

        var A = try Zigen.MatrixXd.init(allocator, 3, 3);
        defer A.deinit();
        const v = [_]f64{ 2, -1, 0, -1, 2, -1, 0, -1, 2 };
        for (0..9) |i| A.data[i] = v[i];
        svd.computeFrom(A);

        std.debug.print("  Singular values: ", .{});
        for (0..3) |i| std.debug.print("{d:.4} ", .{svd.S.data[i]});
        std.debug.print("\n\n", .{});
    }

    // ── Eigenvalue Decomposition ─────────────────────────────────────────
    {
        std.debug.print("── Eigenvalue Decomposition ──\n", .{});
        const Eig = Zigen.SelfAdjointEigenSolverDynamic(f64);

        var eig = try Eig.init(allocator, 3);
        defer eig.deinit();

        var A = try Zigen.MatrixXd.init(allocator, 3, 3);
        defer A.deinit();
        const v = [_]f64{ 4, 2, 1, 2, 5, 3, 1, 3, 6 };
        for (0..9) |i| A.data[i] = v[i];
        eig.computeFrom(A);

        std.debug.print("  Eigenvalues: ", .{});
        for (0..3) |i| std.debug.print("{d:.4} ", .{eig.eigenvalues.data[i]});
        std.debug.print("\n\n", .{});
    }

    std.debug.print("=== All decompositions completed with zero per-call allocations ===\n", .{});
}