iterative_solvers.zig

//! Zigen Example: Iterative Solvers — CG vs BiCGSTAB
//! Comparing convergence on a sparse SPD system (1D Laplacian).

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 — Iterative Solvers (CG vs BiCGSTAB)   ║
        \\  ╚═══════════════════════════════════════════════╝
        \\
        \\
    , .{});

    const n = 50;
    std.debug.print("━━ System: {d}x{d} 1D Laplacian (SPD) ━━━━━━━━━━━━━━━━\n\n", .{ n, n });

    var b_data: [n]f64 = undefined;
    for (0..n) |i| b_data[i] = 1.0;

    // ── CG ────────────────────────────────────────────────────────────
    std.debug.print("━━ Conjugate Gradient ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var cg = Zigen.ConjugateGradient(f64).init(200, 1e-10);
    var dummy: u8 = 0;
    const cg_result = try cg.solve(allocator, n, &tridiagMatvec, @ptrCast(&dummy), &b_data);
    defer allocator.free(cg_result.x);

    std.debug.print("  Iterations: {d}\n", .{cg_result.iterations});
    std.debug.print("  Residual: {d:.10}\n", .{cg_result.residual});
    std.debug.print("  x[0..2] = [{d:.6}, {d:.6}, {d:.6}]\n\n", .{
        cg_result.x[0], cg_result.x[1], cg_result.x[2],
    });

    // ── BiCGSTAB ──────────────────────────────────────────────────────
    std.debug.print("━━ BiCGSTAB ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var bicg = Zigen.BiCGSTAB(f64).init(200, 1e-10);
    const bicg_result = try bicg.solve(allocator, n, &tridiagMatvec, @ptrCast(&dummy), &b_data);
    defer allocator.free(bicg_result.x);

    std.debug.print("  Iterations: {d}\n", .{bicg_result.iterations});
    std.debug.print("  Residual: {d:.10}\n", .{bicg_result.residual});
    std.debug.print("  x[0..2] = [{d:.6}, {d:.6}, {d:.6}]\n\n", .{
        bicg_result.x[0], bicg_result.x[1], bicg_result.x[2],
    });

    // ── Compare ───────────────────────────────────────────────────────
    std.debug.print("━━ Comparison ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\n\n", .{});

    var max_diff: f64 = 0;
    for (0..n) |i| {
        const d = @abs(cg_result.x[i] - bicg_result.x[i]);
        if (d > max_diff) max_diff = d;
    }
    std.debug.print("  Max |CG - BiCGSTAB| = {d:.10}\n", .{max_diff});
    std.debug.print("  CG iters:      {d}\n", .{cg_result.iterations});
    std.debug.print("  BiCGSTAB iters: {d}\n\n", .{bicg_result.iterations});

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

fn tridiagMatvec(x: []const f64, y: []f64, _: *const anyopaque) void {
    const nn = x.len;
    for (0..nn) |i| {
        var val: f64 = 2.0 * x[i];
        if (i > 0) val -= x[i - 1];
        if (i < nn - 1) val -= x[i + 1];
        y[i] = val;
    }
}