const std = @import("std");
const expectEqual = std.testing.expectEqual;
const Matrix = @import("Matrix.zig").Matrix;

test "Matrix LU decomposition" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();

    // Test 3x3 matrix
    var mat: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    mat.at(0, 0).* = 2.0;
    mat.at(0, 1).* = 1.0;
    mat.at(0, 2).* = 1.0;
    mat.at(1, 0).* = 4.0;
    mat.at(1, 1).* = 3.0;
    mat.at(1, 2).* = 3.0;
    mat.at(2, 0).* = 8.0;
    mat.at(2, 1).* = 7.0;
    mat.at(2, 2).* = 9.0;

    // Compute LU decomposition
    var lu_result = try mat.lu();
    defer lu_result.deinit();

    // Verify dimensions
    try expectEqual(@as(usize, 3), lu_result.l.rows);
    try expectEqual(@as(usize, 3), lu_result.l.cols);
    try expectEqual(@as(usize, 3), lu_result.u.rows);
    try expectEqual(@as(usize, 3), lu_result.u.cols);

    // Verify L is lower triangular with 1s on diagonal
    try expectEqual(@as(f64, 1.0), lu_result.l.at(0, 0).*);
    try expectEqual(@as(f64, 1.0), lu_result.l.at(1, 1).*);
    try expectEqual(@as(f64, 1.0), lu_result.l.at(2, 2).*);
    try expectEqual(@as(f64, 0.0), lu_result.l.at(0, 1).*);
    try expectEqual(@as(f64, 0.0), lu_result.l.at(0, 2).*);
    try expectEqual(@as(f64, 0.0), lu_result.l.at(1, 2).*);

    // Verify U is upper triangular
    try expectEqual(@as(f64, 0.0), lu_result.u.at(1, 0).*);
    try expectEqual(@as(f64, 0.0), lu_result.u.at(2, 0).*);
    try expectEqual(@as(f64, 0.0), lu_result.u.at(2, 1).*);

    // Reconstruct PA = LU
    var pa: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    defer pa.deinit();

    // Apply permutation: PA[i,j] = A[p[i],j]
    for (0..3) |i| {
        for (0..3) |j| {
            pa.at(i, j).* = mat.at(lu_result.p.indices[i], j).*;
        }
    }

    // Compute L * U
    var lu_product: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    defer lu_product.deinit();
    @memset(lu_product.items, 0);

    for (0..3) |i| {
        for (0..3) |j| {
            for (0..3) |k| {
                lu_product.at(i, j).* += lu_result.l.at(i, k).* * lu_result.u.at(k, j).*;
            }
        }
    }

    // Verify PA = LU (within numerical tolerance)
    const eps = 1e-10;
    for (0..3) |i| {
        for (0..3) |j| {
            const diff = @abs(pa.at(i, j).* - lu_product.at(i, j).*);
            try std.testing.expect(diff < eps);
        }
    }

    // Verify permutationMatrix() returns correct P such that PA = LU
    var p_mat = try lu_result.permutationMatrix();
    defer p_mat.deinit();

    // Compute P * A
    var pa_from_mat = try p_mat.dot(mat);
    defer pa_from_mat.deinit();

    // Re-verify that rows of PA match rows of L*U (or permuted rows of A)
    for (0..3) |i| {
        for (0..3) |j| {
            try std.testing.expectApproxEqAbs(pa.at(i, j).*, pa_from_mat.at(i, j).*, eps);
        }
    }
}

test "Matrix QR decomposition simple" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();

    // Simple test matrix where column 2 has largest norm
    var mat: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    mat.at(0, 0).* = 1.0;
    mat.at(0, 1).* = 0.0;
    mat.at(0, 2).* = 3.0;
    mat.at(1, 0).* = 0.0;
    mat.at(1, 1).* = 2.0;
    mat.at(1, 2).* = 0.0;
    mat.at(2, 0).* = 0.0;
    mat.at(2, 1).* = 0.0;
    mat.at(2, 2).* = 4.0;

    var qr_result = try mat.qr();
    defer qr_result.deinit();

    // Remove debug print

    // The largest column (column 2) should be first
    try expectEqual(@as(usize, 2), qr_result.perm.indices[0]);
}

test "Matrix QR decomposition" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();

    // Test 3x3 matrix
    var mat: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    mat.at(0, 0).* = 12.0;
    mat.at(0, 1).* = -51.0;
    mat.at(0, 2).* = 4.0;
    mat.at(1, 0).* = 6.0;
    mat.at(1, 1).* = 167.0;
    mat.at(1, 2).* = -68.0;
    mat.at(2, 0).* = -4.0;
    mat.at(2, 1).* = 24.0;
    mat.at(2, 2).* = -41.0;

    // Compute QR decomposition
    var qr_result = try mat.qr();
    defer qr_result.deinit();

    // Verify dimensions
    try expectEqual(@as(usize, 3), qr_result.q.rows);
    try expectEqual(@as(usize, 3), qr_result.q.cols);
    try expectEqual(@as(usize, 3), qr_result.r.rows);
    try expectEqual(@as(usize, 3), qr_result.r.cols);

    // Verify R is upper triangular
    try expectEqual(@as(f64, 0.0), qr_result.r.at(1, 0).*);
    try expectEqual(@as(f64, 0.0), qr_result.r.at(2, 0).*);
    try expectEqual(@as(f64, 0.0), qr_result.r.at(2, 1).*);

    // Verify Q is orthogonal: Q^T * Q should be identity
    var qtq: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    defer qtq.deinit();
    @memset(qtq.items, 0);

    for (0..3) |i| {
        for (0..3) |j| {
            for (0..3) |k| {
                qtq.at(i, j).* += qr_result.q.at(k, i).* * qr_result.q.at(k, j).*;
            }
        }
    }

    // Check that Q^T * Q is approximately identity
    const eps = 1e-10;
    for (0..3) |i| {
        for (0..3) |j| {
            const expected: f64 = if (i == j) 1.0 else 0.0;
            const diff = @abs(qtq.at(i, j).* - expected);
            try std.testing.expect(diff < eps);
        }
    }

    // Verify A * P = Q * R (with column pivoting)
    var qr_product: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    defer qr_product.deinit();
    @memset(qr_product.items, 0);

    for (0..3) |i| {
        for (0..3) |j| {
            for (0..3) |k| {
                qr_product.at(i, j).* += qr_result.q.at(i, k).* * qr_result.r.at(k, j).*;
            }
        }
    }

    // Apply permutation: A * P should equal Q * R
    var ap: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    defer ap.deinit();

    // Apply permutation: A * P = Q * R
    // perm[j] tells us which original column is now at position j
    // So we directly copy column perm[j] of A to position j of AP
    for (0..3) |i| {
        for (0..3) |j| {
            ap.at(i, j).* = mat.at(i, qr_result.perm.indices[j]).*;
        }
    }

    // Verify A*P = Q*R by checking Frobenius norm of difference
    var total_error: f64 = 0;
    for (0..3) |i| {
        for (0..3) |j| {
            const diff = ap.at(i, j).* - qr_product.at(i, j).*;
            total_error += diff * diff;
        }
    }
    const frobenius_error = @sqrt(total_error);

    // Also compute Frobenius norm of A for relative error
    var a_norm: f64 = 0;
    for (0..3) |i| {
        for (0..3) |j| {
            const val = mat.at(i, j).*;
            a_norm += val * val;
        }
    }
    a_norm = @sqrt(a_norm);

    const relative_error = frobenius_error / a_norm;
    if (relative_error >= 1e-8) {
        std.debug.print("\nQR test failing with relative error: {}\n", .{relative_error});
        std.debug.print("Permutation: {} {} {}\n", .{ qr_result.perm.indices[0], qr_result.perm.indices[1], qr_result.perm.indices[2] });

        // Let's check if it's an issue with our test by computing Q*R directly
        std.debug.print("\nDirect check - Q*R:\n", .{});
        for (0..3) |i| {
            for (0..3) |j| {
                std.debug.print("{d:8.2} ", .{qr_product.at(i, j).*});
            }
            std.debug.print("\n", .{});
        }

        std.debug.print("\nA permuted:\n", .{});
        for (0..3) |i| {
            for (0..3) |j| {
                std.debug.print("{d:8.2} ", .{ap.at(i, j).*});
            }
            std.debug.print("\n", .{});
        }
    }
    try std.testing.expect(relative_error < 1e-8);

    // Verify rank is computed correctly (should be 3 for this full-rank matrix)
    try expectEqual(@as(usize, 3), qr_result.rank);

    // Verify permutationMatrix() returns correct P such that AP = QR
    var p_mat = try qr_result.permutationMatrix();
    defer p_mat.deinit();

    // Compute A * P
    var ap_from_mat = try mat.dot(p_mat);
    defer ap_from_mat.deinit();

    // Re-verify that columns of AP match permuted columns of A
    for (0..3) |i| {
        for (0..3) |j| {
            try std.testing.expectApproxEqAbs(ap.at(i, j).*, ap_from_mat.at(i, j).*, 1e-10);
        }
    }

    // Verify columns are ordered by decreasing diagonal values in R (skip for now due to permutation complexity)
    // try std.testing.expect(@abs(qr_result.r.at(0, 0).*) >= @abs(qr_result.r.at(1, 1).*));
    // try std.testing.expect(@abs(qr_result.r.at(1, 1).*) >= @abs(qr_result.r.at(2, 2).*));

    // Test rectangular matrix (4x3) with linearly independent columns
    var rect_mat: Matrix(f64) = try .init(arena.allocator(), 4, 3);
    rect_mat.at(0, 0).* = 1.0;
    rect_mat.at(0, 1).* = 0.0;
    rect_mat.at(0, 2).* = 0.0;
    rect_mat.at(1, 0).* = 1.0;
    rect_mat.at(1, 1).* = 1.0;
    rect_mat.at(1, 2).* = 0.0;
    rect_mat.at(2, 0).* = 1.0;
    rect_mat.at(2, 1).* = 1.0;
    rect_mat.at(2, 2).* = 1.0;
    rect_mat.at(3, 0).* = 1.0;
    rect_mat.at(3, 1).* = 1.0;
    rect_mat.at(3, 2).* = 2.0;

    var rect_qr = try rect_mat.qr();
    defer rect_qr.deinit();

    // Verify dimensions for rectangular matrix
    try expectEqual(@as(usize, 4), rect_qr.q.rows);
    try expectEqual(@as(usize, 3), rect_qr.q.cols);
    try expectEqual(@as(usize, 3), rect_qr.r.rows);
    try expectEqual(@as(usize, 3), rect_qr.r.cols);

    // Verify A * P = Q * R for rectangular matrix
    var rect_product: Matrix(f64) = try .init(arena.allocator(), 4, 3);
    defer rect_product.deinit();
    @memset(rect_product.items, 0);

    for (0..4) |i| {
        for (0..3) |j| {
            for (0..3) |k| {
                rect_product.at(i, j).* += rect_qr.q.at(i, k).* * rect_qr.r.at(k, j).*;
            }
        }
    }

    // Apply permutation to columns of rectangular matrix
    var rect_ap: Matrix(f64) = try .init(arena.allocator(), 4, 3);
    defer rect_ap.deinit();

    for (0..4) |i| {
        for (0..3) |j| {
            rect_ap.at(i, j).* = rect_mat.at(i, rect_qr.perm.indices[j]).*;
        }
    }

    // Verify using relative Frobenius norm
    var rect_error: f64 = 0;
    var rect_norm: f64 = 0;
    for (0..4) |i| {
        for (0..3) |j| {
            const diff = rect_ap.at(i, j).* - rect_product.at(i, j).*;
            rect_error += diff * diff;
            const val = rect_mat.at(i, j).*;
            rect_norm += val * val;
        }
    }
    const rect_relative_error = @sqrt(rect_error) / @sqrt(rect_norm);
    try std.testing.expect(rect_relative_error < 1e-10);

    // Verify rank is 3 for this full-rank rectangular matrix
    try expectEqual(@as(usize, 3), rect_qr.rank);
}

test "Matrix QR decomposition with rank-deficient matrix" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();

    // Create a rank-deficient 4x3 matrix (rank 2)
    // Third column is exactly the sum of first two columns
    var mat: Matrix(f64) = try .init(arena.allocator(), 4, 3);
    mat.at(0, 0).* = 1.0;
    mat.at(0, 1).* = 2.0;
    mat.at(0, 2).* = 3.0; // 1 + 2
    mat.at(1, 0).* = 2.0;
    mat.at(1, 1).* = 3.0;
    mat.at(1, 2).* = 5.0; // 2 + 3
    mat.at(2, 0).* = 3.0;
    mat.at(2, 1).* = 4.0;
    mat.at(2, 2).* = 7.0; // 3 + 4
    mat.at(3, 0).* = 4.0;
    mat.at(3, 1).* = 5.0;
    mat.at(3, 2).* = 9.0; // 4 + 5

    var qr_result = try mat.qr();
    defer qr_result.deinit();

    // Verify rank is 2
    try expectEqual(@as(usize, 2), qr_result.rank);

    // Verify that R has a zero diagonal element at position (2,2)
    const eps = 1e-10;
    try std.testing.expect(@abs(qr_result.r.at(2, 2).*) < eps);

    // Verify A * P = Q * R still holds
    var qr_product: Matrix(f64) = try .init(arena.allocator(), 4, 3);
    defer qr_product.deinit();
    @memset(qr_product.items, 0);

    for (0..4) |i| {
        for (0..3) |j| {
            for (0..3) |k| {
                qr_product.at(i, j).* += qr_result.q.at(i, k).* * qr_result.r.at(k, j).*;
            }
        }
    }

    // Apply permutation
    var ap: Matrix(f64) = try .init(arena.allocator(), 4, 3);
    defer ap.deinit();

    for (0..4) |i| {
        for (0..3) |j| {
            ap.at(i, j).* = mat.at(i, qr_result.perm.indices[j]).*;
        }
    }

    // Verify A*P = Q*R using relative Frobenius norm
    var deficient_error: f64 = 0;
    var deficient_norm: f64 = 0;
    for (0..4) |i| {
        for (0..3) |j| {
            const diff = ap.at(i, j).* - qr_product.at(i, j).*;
            deficient_error += diff * diff;
            const val = mat.at(i, j).*;
            deficient_norm += val * val;
        }
    }
    const deficient_relative_error = @sqrt(deficient_error) / @sqrt(deficient_norm);
    if (deficient_relative_error >= 1e-10) {
        std.debug.print("Deficient test: relative_error = {}, error = {}, norm = {}\n", .{ deficient_relative_error, @sqrt(deficient_error), @sqrt(deficient_norm) });
    }
    try std.testing.expect(deficient_relative_error < 1e-10);

    // Test with zero matrix (rank 0)
    const zero_mat: Matrix(f64) = try .initAll(arena.allocator(), 3, 3, 0);

    var zero_qr = try zero_mat.qr();
    defer zero_qr.deinit();

    // Verify rank is 0
    try expectEqual(@as(usize, 0), zero_qr.rank);
}

test "Matrix rank computation" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();

    // Test 1: Full rank matrix
    var full_rank: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    full_rank.at(0, 0).* = 1.0;
    full_rank.at(0, 1).* = 2.0;
    full_rank.at(0, 2).* = 3.0;
    full_rank.at(1, 0).* = 4.0;
    full_rank.at(1, 1).* = 5.0;
    full_rank.at(1, 2).* = 6.0;
    full_rank.at(2, 0).* = 7.0;
    full_rank.at(2, 1).* = 8.0;
    full_rank.at(2, 2).* = 10.0; // Made it 10 instead of 9 to ensure full rank

    try expectEqual(@as(usize, 3), try full_rank.rank());

    // Test 2: Rank deficient matrix (rank 1)
    // All columns are multiples of the first column
    var rank_1: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    rank_1.at(0, 0).* = 1.0;
    rank_1.at(0, 1).* = 2.0; // 2 * col0
    rank_1.at(0, 2).* = 3.0; // 3 * col0
    rank_1.at(1, 0).* = 2.0;
    rank_1.at(1, 1).* = 4.0; // 2 * col0
    rank_1.at(1, 2).* = 6.0; // 3 * col0
    rank_1.at(2, 0).* = 3.0;
    rank_1.at(2, 1).* = 6.0; // 2 * col0
    rank_1.at(2, 2).* = 9.0; // 3 * col0

    try expectEqual(@as(usize, 1), try rank_1.rank());

    // Test 2b: Rank 2 matrix
    var rank_2: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    rank_2.at(0, 0).* = 1.0;
    rank_2.at(0, 1).* = 0.0;
    rank_2.at(0, 2).* = 1.0; // col2 = col0
    rank_2.at(1, 0).* = 0.0;
    rank_2.at(1, 1).* = 1.0;
    rank_2.at(1, 2).* = 0.0; // col2 = col0
    rank_2.at(2, 0).* = 0.0;
    rank_2.at(2, 1).* = 0.0;
    rank_2.at(2, 2).* = 0.0; // col2 = col0

    try expectEqual(@as(usize, 2), try rank_2.rank());

    // Test 3: Zero matrix (rank 0)
    const zero_mat: Matrix(f64) = try .initAll(arena.allocator(), 4, 3, 0);
    try expectEqual(@as(usize, 0), try zero_mat.rank());

    // Test 4: Rectangular matrix with rank deficiency
    var rect_mat: Matrix(f64) = try .init(arena.allocator(), 5, 3);
    // Make columns 0 and 1 independent, column 2 = column0 + column1
    rect_mat.at(0, 0).* = 1.0;
    rect_mat.at(0, 1).* = 0.0;
    rect_mat.at(0, 2).* = 1.0; // col0 + col1
    rect_mat.at(1, 0).* = 0.0;
    rect_mat.at(1, 1).* = 1.0;
    rect_mat.at(1, 2).* = 1.0; // col0 + col1
    rect_mat.at(2, 0).* = 1.0;
    rect_mat.at(2, 1).* = 1.0;
    rect_mat.at(2, 2).* = 2.0; // col0 + col1
    rect_mat.at(3, 0).* = 0.0;
    rect_mat.at(3, 1).* = 2.0;
    rect_mat.at(3, 2).* = 2.0; // col0 + col1
    rect_mat.at(4, 0).* = 2.0;
    rect_mat.at(4, 1).* = 1.0;
    rect_mat.at(4, 2).* = 3.0; // col0 + col1

    try expectEqual(@as(usize, 2), try rect_mat.rank());

    // Test 5: Single element matrix
    const single: Matrix(f64) = try .initAll(arena.allocator(), 1, 1, 5.0);
    try expectEqual(@as(usize, 1), try single.rank());

    // Test 6: Column vector
    var col_vec: Matrix(f64) = try .init(arena.allocator(), 5, 1);
    for (0..5) |i| {
        col_vec.at(i, 0).* = @floatFromInt(i + 1);
    }
    try expectEqual(@as(usize, 1), try col_vec.rank());

    // Test 7: Row vector
    var row_vec: Matrix(f64) = try .init(arena.allocator(), 1, 5);
    for (0..5) |i| {
        row_vec.at(0, i).* = @floatFromInt(i + 1);
    }
    try expectEqual(@as(usize, 1), try row_vec.rank());
}

test "Matrix Cholesky decomposition" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();
    const MatrixError = @import("Matrix.zig").MatrixError;

    // Test 3x3 symmetric positive-definite matrix
    var mat: Matrix(f64) = try .init(arena.allocator(), 3, 3);
    // A = L * L^T
    // L = [[2, 0, 0], [1, 2, 0], [1, 1, 2]]
    // A = [[4, 2, 2], [2, 5, 3], [2, 3, 6]]
    mat.at(0, 0).* = 4.0;
    mat.at(0, 1).* = 2.0;
    mat.at(0, 2).* = 2.0;
    mat.at(1, 0).* = 2.0;
    mat.at(1, 1).* = 5.0;
    mat.at(1, 2).* = 3.0;
    mat.at(2, 0).* = 2.0;
    mat.at(2, 1).* = 3.0;
    mat.at(2, 2).* = 6.0;

    var chol = try mat.chol();
    defer chol.deinit();

    // Check L
    const eps = 1e-10;
    const expected_l_data = [_]f64{
        2.0, 0.0, 0.0,
        1.0, 2.0, 0.0,
        1.0, 1.0, 2.0,
    };
    var expected_l: Matrix(f64) = try .fromSlice(arena.allocator(), 3, 3, &expected_l_data);
    defer expected_l.deinit();

    for (0..3) |i| {
        for (0..3) |j| {
            try std.testing.expectApproxEqAbs(expected_l.at(i, j).*, chol.at(i, j).*, eps);
        }
    }

    // Verify L * L^T = A
    var lt = try chol.transpose();
    defer lt.deinit();
    var recon = try chol.dot(lt);
    defer recon.deinit();

    for (0..3) |i| {
        for (0..3) |j| {
            try std.testing.expectApproxEqAbs(mat.at(i, j).*, recon.at(i, j).*, eps);
        }
    }

    // Test non-positive definite matrix
    var non_spd: Matrix(f64) = try .init(arena.allocator(), 2, 2);
    non_spd.at(0, 0).* = 1.0;
    non_spd.at(0, 1).* = 2.0;
    non_spd.at(1, 0).* = 2.0;
    non_spd.at(1, 1).* = 1.0; // Det = 1 - 4 = -3, not positive definite

    try std.testing.expectError(MatrixError.NotPositiveDefinite, non_spd.chol());
}

test "Matrix solve single right-hand side" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();
    const alloc = arena.allocator();
    const eps = 1e-12;

    // Solution is x = [1, 2, 3].
    const a: Matrix(f64) = try .fromSlice(alloc, 3, 3, &[_]f64{
        2, 1, 1,
        4, 3, 3,
        8, 7, 9,
    });
    const b: Matrix(f64) = try .fromSlice(alloc, 3, 1, &[_]f64{ 7, 19, 49 });

    const x = try a.solve(b);
    try std.testing.expectApproxEqAbs(@as(f64, 1), x.at(0, 0).*, eps);
    try std.testing.expectApproxEqAbs(@as(f64, 2), x.at(1, 0).*, eps);
    try std.testing.expectApproxEqAbs(@as(f64, 3), x.at(2, 0).*, eps);

    const recon = try a.dot(x);
    for (0..3) |i| try std.testing.expectApproxEqAbs(b.at(i, 0).*, recon.at(i, 0).*, eps);
}

test "Matrix solve multiple right-hand sides via reused factorization" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();
    const alloc = arena.allocator();
    const eps = 1e-12;

    const a: Matrix(f64) = try .fromSlice(alloc, 2, 2, &[_]f64{
        3, 2,
        1, 4,
    });
    // RHS = I, so the solution is A^-1 and A*x == I.
    const b: Matrix(f64) = try .fromSlice(alloc, 2, 2, &[_]f64{
        1, 0,
        0, 1,
    });

    var lu_result = try a.lu();
    defer lu_result.deinit();
    const x = try lu_result.solve(b);

    const recon = try a.dot(x);
    for (0..2) |i| {
        for (0..2) |j| {
            const expected: f64 = if (i == j) 1 else 0;
            try std.testing.expectApproxEqAbs(expected, recon.at(i, j).*, eps);
        }
    }
}

test "Matrix solve reports singular systems" {
    var arena: std.heap.ArenaAllocator = .init(std.testing.allocator);
    defer arena.deinit();
    const alloc = arena.allocator();
    const MatrixError = @import("Matrix.zig").MatrixError;

    // Second row is twice the first: singular.
    const a: Matrix(f64) = try .fromSlice(alloc, 2, 2, &[_]f64{
        1, 2,
        2, 4,
    });
    const b: Matrix(f64) = try .fromSlice(alloc, 2, 1, &[_]f64{ 1, 2 });
    try std.testing.expectError(MatrixError.Singular, a.solve(b));

    const b_bad: Matrix(f64) = try .fromSlice(alloc, 3, 1, &[_]f64{ 1, 2, 3 });
    try std.testing.expectError(MatrixError.DimensionMismatch, a.solve(b_bad));
}
