zigen/examples/linear_algebra.zig at main · coderonion/zigen · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
//! Zigen Example: Decompositions Showcase
//! LU, QR, SVD, Cholesky, and Eigenvalues with factorization verification.

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

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

    const A = Zigen.Matrix3f.fromArray(.{
        .{ 4.0, 2.0, 1.0 },
        .{ 2.0, 5.0, 3.0 },
        .{ 1.0, 3.0, 6.0 },
    });
    std.debug.print("Matrix A:\n", .{});
    printMat3(A);

    // ── 1. LU ─────────────────────────────────────────────────────────
    std.debug.print("━━ 1. LU: A = L*U ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const lu = try Zigen.LU(f32, 3).compute(A);
    std.debug.print("L:\n", .{});
    printMat3(lu.L);
    std.debug.print("U:\n", .{});
    printMat3(lu.U);

    const LU_prod = lu.L.mul(3, lu.U);
    std.debug.print("  ||L*U - A|| = {d:.6}\n", .{LU_prod.sub(A).norm()});
    std.debug.print("  det(A) = {d:.2}\n\n", .{lu.determinant()});

    // ── 2. QR ─────────────────────────────────────────────────────────
    std.debug.print("━━ 2. QR: A = Q*R ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const qr = Zigen.QR(f32, 3, 3).compute(A);
    std.debug.print("Q:\n", .{});
    printMat3(qr.Q);
    std.debug.print("R:\n", .{});
    printMat3(qr.R);

    const QR_prod = qr.Q.mul(3, qr.R);
    std.debug.print("  ||Q*R - A|| = {d:.6}\n", .{QR_prod.sub(A).norm()});
    const QtQ = qr.Q.transpose().mul(3, qr.Q);
    std.debug.print("  ||Q^T*Q - I|| = {d:.6}\n\n", .{QtQ.sub(Zigen.Matrix3f.identity()).norm()});

    // ── 3. SVD ────────────────────────────────────────────────────────
    std.debug.print("━━ 3. SVD: A = U*S*V^T ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const svd = Zigen.SVD(f32, 3, 3).compute(A);
    std.debug.print("  Singular values: {d:.4}  {d:.4}  {d:.4}\n", .{
        svd.S.at(0), svd.S.at(1), svd.S.at(2),
    });
    const cond = svd.S.at(0) / svd.S.at(2);
    std.debug.print("  Condition number: {d:.2}\n\n", .{cond});

    // ── 4. Cholesky ───────────────────────────────────────────────────
    std.debug.print("━━ 4. Cholesky: A = L*L^T ━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const chol = try Zigen.Cholesky(f32, 3).compute(A);
    std.debug.print("L:\n", .{});
    printMat3(chol.L);

    const LLt = chol.L.mul(3, chol.L.transpose());
    std.debug.print("  ||L*L^T - A|| = {d:.6}\n\n", .{LLt.sub(A).norm()});

    // ── 5. Eigenvalues ────────────────────────────────────────────────
    std.debug.print("━━ 5. Eigenvalues (Symmetric A) ━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    const eig = Zigen.SelfAdjointEigenSolver(f32, 3).compute(A);
    const ev = eig.eigenvalues();
    std.debug.print("  Eigenvalues: {d:.4}  {d:.4}  {d:.4}\n", .{
        ev.at(0), ev.at(1), ev.at(2),
    });

    const eig_sum = ev.at(0) + ev.at(1) + ev.at(2);
    std.debug.print("  Sum = {d:.4} (trace = {d:.4})\n", .{ eig_sum, A.trace() });

    const eig_prod = ev.at(0) * ev.at(1) * ev.at(2);
    std.debug.print("  Product = {d:.4} (det = {d:.4})\n\n", .{ eig_prod, lu.determinant() });

    std.debug.print("All decompositions verified!\n", .{});
}

fn printMat3(m: Zigen.Matrix3f) void {
    var i: usize = 0;
    while (i < 3) : (i += 1) {
        std.debug.print("  | {d:>9.4} {d:>9.4} {d:>9.4} |\n", .{ m.at(i, 0), m.at(i, 1), m.at(i, 2) });
    }
    std.debug.print("\n", .{});
}