Repository navigation
[BUG] matmul produces incorrect results for F-contiguous / non-row-major inputs #89
Description
np.matmul on QuadPrecision arrays returns incorrect results whenever either input is F-contiguous (or otherwise not row-major-contiguous in its last two axes).
Reproducer
click to expand
import numpy as np
from numpy_quaddtype import QuadPrecision, QuadPrecDType
dt = QuadPrecDType(backend='sleef')
def Q(v): return QuadPrecision(str(float(v)), backend='sleef')
A = np.asfortranarray(np.array([[Q(1), Q(2), Q(3)],
[Q(4), Q(5), Q(6)]], dtype=dt))
B = np.asfortranarray(np.array([[Q(1), Q(2)],
[Q(3), Q(4)],
[Q(5), Q(6)]], dtype=dt))
# A is F-contig: strides=(16, 32); B is F-contig: strides=(16, 48)
print("quad F@F:", [[float(v) for v in r] for r in (A @ B)])
# -> [[23.0, 27.0], [35.0, 32.0]] <- WRONG
# Same call against native dtypes works fine — proving the bug is on our side:
A32 = np.asfortranarray(np.array([[1., 2., 3.], [4., 5., 6.]], dtype=np.float32))
B32 = np.asfortranarray(np.array([[1., 2.], [3., 4.], [5., 6.]], dtype=np.float32))
print("float32 F@F:", (A32 @ B32).tolist())
# -> [[22.0, 28.0], [49.0, 64.0]] <- correct
# Workaround: force C-contig before matmul
print("quad (C@C):",
[[float(v) for v in r] for r in (np.ascontiguousarray(A) @ np.ascontiguousarray(B))])
# -> [[22.0, 28.0], [49.0, 64.0]] <- correctRoot cause
quad_matmul_strided_loop_aligned in src/csrc/umath/matmul.cpp unconditionally:
size_t lda = A_row_stride / sizeof(Sleef_quad);
// ...
qblas_gemm('R', 'N', 'N', m, p, n, &alpha, A_ptr, lda, B_ptr, ldb,
&beta, C_ptr, ldc_numpy);This silently assumes:
- the input is row-major contiguous in its core dimensions and the right layout argument to QBLAS is always
'R'.
Possible fixes
- Inspect each operand's last-two-axis strides: dispatch directly to
qblas_gemm('R', ...)or('C', ...)when both operands share a contiguous layout, otherwise memcpy the offending operand(s) into a row-major temp buffer (same trick the unaligned path already uses) and callqblas_gemm('R', ...). - Performance (needs QBLAS-side change too). For mixed-layout inputs (R @ C and C @ R) the copy above can be skipped by passing
transa/transb='T'instead, which would requireqblas_gemmto gain transpose support. See [FEAT] Support transpose flagstransa, transb != 'N'in GEMM SwayamInSync/QBLAS#13
Related tests
On the #88 two tests already exercise this bug and are marked xfail:
tests/test_dot.py::TestMatmulBatched::test_batched_fortran_ordertests/test_dot.py::TestMatmulBatched::test_batched_transposed_view
Activity
SwayamInSync commentedon May 12, 2026
If not urgent, then I will be happy to add the slow copy path for correctness, and test the transpose kernel perf as supported added in QBLAS