a port of https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/spectralnorm-gcc-5.html to zig

spectralnorm.zig

// $ zig build-exe --release-fast spectralnorm.zig && time ./spectralnorm 5500


const std = @import("std");
const warn = std.debug.warn;

const V2f64 = @Vector(2, f64);

inline fn A(_i: usize, _j: usize) f64 {
    var i: f64 = @intToFloat(f64, _i);
    var j: f64 = @intToFloat(f64, _j);
    return ((i + j) * (i + j + 1) / 2 + i + 1);
}

fn dot(v: []f64, u: []f64, n: usize) f64 {
    var i: usize = 0;
    var sum: f64 = 0;
    while (i < n) : (i += 1)
        sum += v[i] * u[i];
    return sum;
}

fn mult_Av(v: []f64, out: []f64, n: usize) void {
    var i: usize = 0;

    // #pragma omp parallel for
    while (i < n) : (i += 1) {
        var sum: V2f64 = .{ 0, 0 };

        var j: usize = 0;
        while (j < n) : (j += 2) {
            const b: V2f64 = .{ v[j], v[j + 1] };
            const a: V2f64 = .{ A(i, j), A(i, j + 1) };
            sum += b / a;
        }
        out[i] = sum[0] + sum[1];
    }
}

fn mult_Atv(v: []f64, out: []f64, n: usize) void {
    var i: usize = 0;

    // #pragma omp parallel for

    while (i < n) : (i += 1) {
        var sum: V2f64 = .{ 0, 0 };

        var j: usize = 0;
        while (j < n) : (j += 2) {
            const b: V2f64 = .{ v[j], v[j + 1] };
            const a: V2f64 = .{ A(j, i), A(j + 1, i) };
            sum += b / a;
        }
        out[i] = sum[0] + sum[1];
    }
}

fn mult_AtAv(v: []f64, out: []f64, tmp: []f64, n: usize) void {
    mult_Av(v, tmp, n);
    mult_Atv(tmp, out, n);
}

pub fn main() anyerror!void {
    const args = try std.process.argsAlloc(std.heap.page_allocator);
    defer std.process.argsFree(std.heap.page_allocator, args);
    var n = std.fmt.parseInt(usize, args[1], 10) catch 2000;
    // warn("{}\n", .{n});
    if (n & 1 == 1) n += 1;
    const a = std.heap.page_allocator;
    var u = try a.alloc(f64, n);
    var v = try a.alloc(f64, n);
    var tmp = try a.alloc(f64, n);
    for (u) |_, i| u[i] = 1;
    var i: u8 = 0;
    while (i < 10) : (i += 1) {
        mult_AtAv(u, v, tmp, n);
        mult_AtAv(v, u, tmp, n);
    }
    warn("{d:.9}\n", .{std.math.sqrt(dot(u, v, n) / dot(v, v, n))});
}