Describe the issue:

I expected np.matvec(A, b, out=b) to behave more or less like np.einsum('...ij,...j', A, b, out=b), but instead b is set to all zeros

It seems to work fine if the output vector is not the same as the input vector.

Reproduce the code example:

>>> import numpy as np
>>> M, N = 100000, 10
... A = np.random.normal(size=(M, N, N))
... b = np.random.normal(size=(M, N))
>>> np.matvec(A, b, out=b)
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]], shape=(100000, 10))

Error message:

Python and NumPy Versions:

2.4.4
3.14.3 (main, Feb 15 2026, 00:49:34) [GCC 15.2.1 20260209]

Runtime Environment:

[{'numpy_version': '2.4.4',
'python': '3.14.3 (main, Feb 15 2026, 00:49:34) [GCC 15.2.1 20260209]',
'uname': uname_result(system='Linux', node='...', release='6.19.10-artix1-1', version='#1 SMP PREEMPT_DYNAMIC Thu, 26 Mar 2026 22:46:44 +0000', machine='x86_64')},
{'simd_extensions': {'baseline': [],
'found': ['X86_V2', 'X86_V3'],
'not_found': ['X86_V4', 'AVX512_ICL', 'AVX512_SPR']}},
{'ignore_floating_point_errors_in_matmul': False},
{'filepath': '/usr/lib/libgomp.so.1.0.0',
'internal_api': 'openmp',
'num_threads': 12,
'prefix': 'libgomp',
'user_api': 'openmp',
'version': None}]

How does this issue affect you or how did you find it:

I needed to apply a linear transformation to many vectors in-place, and it didn't work.