Broken matrix-matrix multiplication with single-column matrices

[nvb]

Hi maintainers, thanks for the beautiful library!

I've noticed a bug with matrix-matrix multiplication that materializes only with single-column matrices. Explanation and reproduction below.

Description:

• Single-column matrices (shape [N, 1]) produce all-zero results when used with matrix-matrix multiplication (dot product) in ndarray 0.16.1 with BLAS enabled
• The bug affects matrices with 16+ dimensions but not smaller ones (3D works fine)
• Both column-major and row-major layouts are affected
• Matrix-vector multiplication (SGEMV) works correctly, but matrix-matrix multiplication (SGEMM) with single-column matrices fails

Reproduction:

#[test]
fn test_ndarray_sgemm_bug_simple() {
    // Minimal ndarray bug reproduction: single-column matrices return all zeros
    use ndarray::prelude::*;

    let n = 64; // Bug manifests at 16+ dimensions
    let matrix = Array2::<f32>::eye(n); // Identity matrix
    let vector: Vec<f32> = (1..=n).map(|i| i as f32).collect();

    // This works (SGEMV):
    let correct = matrix.dot(&ArrayView1::from(&vector));

    // This fails (SGEMM with single-column matrix):
    let single_col = ArrayView2::from_shape((n, 1).f(), &vector).unwrap();
    let buggy = matrix.dot(&single_col);

    println!("SGEMV result: {:?}", &correct.to_vec()[..3]);
    println!("SGEMM result: {:?}", &buggy.column(0).to_vec()[..3]);

    // Should be [1.0, 2.0, 3.0, ...] but SGEMM gives [0.0, 0.0, 0.0, ...]
    let all_zeros = buggy.column(0).iter().all(|&x| x == 0.0);
    if all_zeros {
        println!("🐛 BUG: SGEMM with single-column matrix returns all zeros!");
    }
}

[nvb]

I'm only seeing this on my M4 Mac and not on a Linux machine, so it may actually be a problem with my OpenBLAS installation.

[nvb]

This works (with Accelerate BLAS)

RUSTFLAGS="-C link-arg=-framework -C link-arg=Accelerate" cargo test --features blas

but this does not (with OpenBlas on Mac)

RUSTFLAGS="-L/opt/homebrew/opt/openblas/lib -lopenblas" cargo test --features blas

[nvb]

Interestingly, test_mat_mul (on master) also fails with OpenBlas 0.3.30.

[nvb]

Here's a reproduction in C, courtesy of my daily quotas of Cursor Agent and Claude Code:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// Include BLAS headers - OpenBLAS or Accelerate
#ifdef __APPLE__
    #include <Accelerate/Accelerate.h>
#else
    #include <cblas.h>
#endif

int main() {
    // Test case from ndarray test_mat_mul
    // A is 45x33, B is 33x33 identity, C should be 45x33 (same as A)
    
    const int m = 45, n = 33, k = 33;
    const float alpha = 1.0f, beta = 0.0f;
    
    // Allocate matrices
    float *A = (float*)malloc(m * k * sizeof(float));
    float *B = (float*)malloc(k * n * sizeof(float));  
    float *C = (float*)malloc(m * n * sizeof(float));
    
    // Initialize A with linspace(0, 1484, 45*33) like ndarray test
    for (int i = 0; i < m * k; i++) {
        A[i] = (float)i;
    }
    
    // Initialize B as identity matrix 
    memset(B, 0, k * n * sizeof(float));
    for (int i = 0; i < k; i++) {
        B[i * n + i] = 1.0f;  // B[i][i] = 1
    }
    
    // Initialize C to zero
    memset(C, 0, m * n * sizeof(float));
    
    printf("Before SGEMM:\n");
    printf("A[0..2,0..2]:\n");
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            printf("%.2f ", A[i * k + j]);
        }
        printf("\n");
    }
    
    printf("B[0..2,0..2]:\n");
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            printf("%.2f ", B[i * n + j]);
        }
        printf("\n");
    }
    
    // Call OpenBLAS SGEMM
    // C = alpha * A * B + beta * C
    // Row-major layout: CblasRowMajor, no transpose
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                m, n, k,
                alpha,
                A, k,    // lda = k (leading dimension of A)
                B, n,    // ldb = n (leading dimension of B)
                beta,
                C, n);   // ldc = n (leading dimension of C)
    
    printf("\nAfter SGEMM:\n");
    printf("C[0..2,0..2] (should equal A[0..2,0..2]):\n");
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < 3; j++) {
            printf("%.2f ", C[i * n + j]);
        }
        printf("\n");
    }
    
    // Check if result is correct (C should equal A since A * I = A)
    int correct = 1;
    for (int i = 0; i < m * n; i++) {
        if (fabs(C[i] - A[i]) > 1e-5) {
            correct = 0;
            break;
        }
    }
    
    printf("\nResult: %s\n", correct ? "CORRECT" : "INCORRECT");
    
    // Print some diagnostics
    printf("\nDiagnostics:\n");
    printf("All C values zero? %s\n", C[0] == 0.0f && C[1] == 0.0f ? "YES" : "NO");
    
    free(A);
    free(B);
    free(C);
    
    return correct ? 0 : 1;
}

OpenBlas (0.3.30)

~/rust/ndarray % gcc -I/opt/homebrew/opt/openblas/include -L/opt/homebrew/opt/openblas/lib -lopenblas debug_openblas.c -o debug_openblas -Wno-deprecated
./debug_openblas
Before SGEMM:
A[0..2,0..2]:
0.00 1.00 2.00
33.00 34.00 35.00
66.00 67.00 68.00
B[0..2,0..2]:
1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00

After SGEMM:
C[0..2,0..2] (should equal A[0..2,0..2]):
0.00 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.00

Result: INCORRECT

Diagnostics:
All C values zero? YES

Accelerate (15.5)

~/rust/ndarray % gcc -framework Accelerate debug_openblas.c -o debug_accelerate -Wno-deprecated
./debug_accelerate
Before SGEMM:
A[0..2,0..2]:
0.00 1.00 2.00
33.00 34.00 35.00
66.00 67.00 68.00
B[0..2,0..2]:
1.00 0.00 0.00
0.00 1.00 0.00
0.00 0.00 1.00

After SGEMM:
C[0..2,0..2] (should equal A[0..2,0..2]):
0.00 1.00 2.00
33.00 34.00 35.00
66.00 67.00 68.00

Result: CORRECT

Diagnostics:
All C values zero? NO

[nvb]

I opened OpenMathLib/OpenBLAS#5429 in the OpenBLAS repo to track this. I'll keep this issue open (if you don't mind) in case it turns out that we're not using OpenBLAS as expected. Or it case it turns out that someone here has any insight into what might be going wrong.

[nvb]

Closing, as this does seem to be a bug in OpenBLAS on M4 macs, likely due to their new SME-based fast path for small-matrix SGEMM that was added in 0.3.30.