Repository: numpy/numpy
Pull request: BUG: fix matmul with transposed out arg
Fixes #29164, the fix was suggested by @seberg

From the last comment in the issue:

> Ideally, make sure we have tests for all other cases as well.

I think that input is well tested via `numpy/_core/tests/test_multiarray.py::TestMatmul`'s `test_dot_equivalent`, and there are some tests in `test_out_contiguous` but that the specific case we hit here is missing. I added it: the test crashes before the fix, and passes after.

Problem statement:
BUG: Incorrect result for `numpy.matmul` when `order="F"` is passed
### Describe the issue:

`numpy.matmul` returns incorrect result when `order="F"` with numpy-2.3.0. With numpy-2.2.6, it was working fine.

### Reproduce the code example:

```python
import numpy
print("numpy version: ", numpy.__version__)
a = numpy.arange(24, dtype="f8").reshape((4,2,3))
b = numpy.arange(60, dtype="f8").reshape((4,3,5))
res = numpy.matmul(a, b, order="F")
print(res[0, ...])
# numpy version:  2.3.0
#  ** On entry to DGEMM  parameter number 13 had an illegal value
#  ** On entry to DGEMM  parameter number 13 had an illegal value
#  ** On entry to DGEMM  parameter number 13 had an illegal value
#  ** On entry to DGEMM  parameter number 13 had an illegal value
# [[5.43472924e-310 5.43508168e-310 6.85110253e-310 6.85110253e-310
#   6.85110253e-310]
#  [0.00000000e+000 5.43508169e-310 5.43508169e-310 6.85110251e-310
#   6.85110253e-310]]
```

### Error message:

```shell

```

### Python and NumPy Versions:

2.3.0
3.12.11 | packaged by conda-forge | (main, Jun  4 2025, 14:45:31) [GCC 13.3.0]

### Runtime Environment:

_No response_

### Context for the issue:

_No response_

Hints from issue discussion:
What platform are you on?
> What platform are you on?

Running on Ubuntu 24.04.
And a different but similar in concept failing example

```python
import numpy
print("numpy version: ", numpy.__version__)
A = numpy.arange(100, dtype="f8").reshape((10, 10))
a = A[::-2, ::-2]
OUT = numpy.empty((10, 10), dtype="f8")
out = OUT[::-2, ::-2]
numpy.matmul(a, a, out=out)
print(out[0, ...])
# numpy version:  2.3.0
#  ** On entry to DGEMM  parameter number 13 had an illegal value
# [0. 0. 0. 0. 0.]
```

that again was working fine with numpy-2.2
```python
import numpy
print("numpy version: ", numpy.__version__)
A = numpy.arange(100, dtype="f8").reshape((10, 10))
a = A[::-2, ::-2]
OUT = numpy.empty((10, 10), dtype="f8")
out = OUT[::-2, ::-2]
numpy.matmul(a, a, out=out)
print(out[0, ...])
# numpy version:  2.2.6
# [28425. 27475. 26525. 25575. 24625.]
```
> > What platform are you on?
> 
> Running on Ubuntu 24.04.

What kind of machine?
We did not change OpenBLAS for 2.3.0, both 2.3.0 and 2.2.6 use `OpenBLAS 0.3.29`. I can reproduce the failure of the first script on Python3.11 using NumPy from PyPI on my standard AMD x86_64 Ubuntu 24.04 machine. 

The other reproducers here are missing a definition for `shape`. 

Did we change something in `matmul` or `dot`? The only thing I see in the release notes is #28442.
> The other reproducers here are missing a definition for `shape`.

Fixed.

Running this on a local build of v2.3.0 without OpenBLAS does not fail.
Building with OpenBLAS like this fails
```
pip install "scipy-openblas64==0.3.29.0"
PKG_CONFIG_PATH=$PWD/.openblas spin build --clean --with-scipy-openblas=64
spin python /tmp/issue29164.py
```

but differently. I wonder ~if~ how I catch the assert when the wheel build does not:
```
python3: ../numpy/_core/src/umath/matmul.c.src:176: 
   DOUBLE_matmul_matrixmatrix: Assertion `is_blasable2d(os_m, os_p, m, p, sizeof(npy_double))' failed.
```

In any case, that file was changed as part of #23752, which seems to have touched the non-contiguous matmul, making it quite suspect. @xor2k thoughts?
> but differently. I wonder if how I catch the assert when the wheel build does not:

I think the `spin` build defaults are just different and keep asserts included.  Anyway, seems useful here, since that is surely the bug (very surprising that we don't have tests that caught this :().
Indeed, bisecting around f10429192b points to that being the commit that starts failing.
The code ends up using BLAS in DOUBLE_matmul here https://github.com/numpy/numpy/blob/30a4883d32a720a4af8a273d34673e671819e7b1/numpy/_core/src/umath/matmul.c.src#L608

when `o_f_blasable` is `0`. Previously, the condition `o_f_blasable==0` would have led to using a  `DOUBLE_matmul_inner_noblas`. The logic is pretty tangled, no wonder we got it wrong.

Commit message: BUG: fix matmul with transposed out arg

Patch to review:
diff --git a/numpy/_core/src/umath/matmul.c.src b/numpy/_core/src/umath/matmul.c.src
index d9be7b1d6826..02c4fde56bf2 100644
--- a/numpy/_core/src/umath/matmul.c.src
+++ b/numpy/_core/src/umath/matmul.c.src
@@ -596,7 +596,7 @@ NPY_NO_EXPORT void
              * Use transpose equivalence:
              * matmul(a, b, o) == matmul(b.T, a.T, o.T)
              */
-            if (o_f_blasable) {
+            if (o_transpose) {
                 @TYPE@_matmul_matrixmatrix(
                     ip2_, is2_p_, is2_n_,
                     ip1_, is1_n_, is1_m_,
diff --git a/numpy/_core/tests/test_multiarray.py b/numpy/_core/tests/test_multiarray.py
index 7603449ba28e..34740963f6d8 100644
--- a/numpy/_core/tests/test_multiarray.py
+++ b/numpy/_core/tests/test_multiarray.py
@@ -7272,6 +7272,10 @@ def test_out_contiguous(self):
         assert_array_equal(c, tgt_mv)
         c = self.matmul(v, a.T, out=out[:, 0, 0])
         assert_array_equal(c, tgt_mv)
+        # issue 29164
+        out_f = np.zeros((10, 4), dtype=float)
+        c = self.matmul(a, b, out=out_f[::-2, ::-2])
+        assert_array_equal(c, tgt)
 
         # test out contiguous in only last dim
         out = np.ones((10, 2), dtype=float)


Code review comments:
diff_hunk: @@ -7272,6 +7272,10 @@ def test_out_contiguous(self):
         assert_array_equal(c, tgt_mv)
         c = self.matmul(v, a.T, out=out[:, 0, 0])
         assert_array_equal(c, tgt_mv)
+        # issue 29164
+        out_f = np.zeros((10, 4), dtype=float)
+        c = self.matmul(a, b, out=out_f[::-2, ::-2])
path: numpy/_core/tests/test_multiarray.py
text: @user1:
```suggestion
        out_f = np.zeros((10, 4), dtype=float, order="F")
        c = self.matmul(a, b, out=out_f[::2, ::2])
```
Or did it really fail before?

@user1:
Ohhhhh, the test uses `::-2` so does fail, but that is actually another bug in the code (but not a critical one, just a performance one).
Strides should be compared by their absolute value to decide whether to transpose.

@user1:
@user2 do you want to check on the stride logic here?  Seems a bit awkward to iterate this in fortran order due to negative strides, when the memory order is more like C otherwise.

@user2:
Stride comparisons happen only in lines 557 to 559 of `matmul.c.src`. There all left and right sides should be `abs`, let me check...

@user2:
Looks like I can't make comments on lines to be changed.

@user2:
So this is how it needs to look like
```python
            npy_bool i1_transpose = labs(is1_m) < labs(is1_n),
                     i2_transpose = labs(is2_n) < labs(is2_p),
                     o_transpose = labs(os_m) < labs(os_p);
```
I'll just make a dedicated merge request for that.

@user2:
I mean should I or would it be more convenient to have it in this pull request, too?

Pardon the "merge request" above, I'm using gitlab at work 😅 

@user2:
I have an idea where to move this test to make it more thorough: at the end of function `test_dot_equivalent` in line 7320 (w.r.t. the unmodified file):
```python
        # matrix matrix, issue 29164
        if [len(args[0].shape), len(args[1].shape)] == [2, 2]:
            out_f = np.zeros((r2.shape[0] * 2, r2.shape[1] * 2), order='F')
            r4 = np.matmul(*args, out=out_f[::2, ::2])
            assert_equal(r2, r4)
```

How about that? Then we would have all the combinations of c/f/non-contiguous inputs and outputs if I'm not mistaken and nothing is missing in the `pytest.mark.parametrize` args. Plus, the current suggestion above is basically a case of `test_dot_equivalent`, since it compares an output of `dot` with an output of `matmul`. Default `dtype` is `float64`, so it will trigger BLAS, too.

@author:
OK.