Repository: numpy/numpy
Pull request: BUG: Use output when given on numpy.dot C-API branch
Updated the numpy.dot C-API so that it now copies the result to the out parameter and returns it on branch of the function where the correct behaviour was not in place. Added two tests regarding this issue. Closes #21081.

Problem statement:
BUG: scalar.dot(vector) doesn't fill out-array of dtype float16
### Describe the issue:

The output array in the code snippet below should be getting filled with the result of the computation, but it remains unchanged.

### Reproduce the code example:

```python
import numpy as np
x = np.ones((), dtype=np.float16)
y = np.ones((5,), dtype=np.float16)
z = np.zeros((5,), dtype=np.float16)
res = x.dot(y, out=z)
assert np.array_equal(res, y)  # succeeds
assert np.array_equal(z, y)  # fails
```


### Error message:

```shell
Traceback (most recent call last):
  File "/home/mpapadakis/out-half.py", line 7, in <module>
    assert np.array_equal(z, y)  # fails
AssertionError
```


### NumPy/Python version information:

1.22.2 3.10.2 | packaged by conda-forge | (main, Feb  1 2022, 19:29:00) [GCC 9.4.0]

Hints from issue discussion:
It looks like there are other unsupported types `np.int8/16/32/64`.

Tested on numpy 1.22.2, Python 3.10.2 (v3.10.2:a58ebcc701, Jan 13 2022, 14:50:16) [Clang 13.0.0 (clang-1300.0.29.30)].
I'm fairly certain this issue is dtype independent. The relevant function has a dedicated branch for dealing with scalar dot products, but it completely ignores the `out` parameter doing so:

https://github.com/numpy/numpy/blob/a33a10a5f494057b9ddcb9efe650306fccf8eba3/numpy/core/src/multiarray/multiarraymodule.c#L1018-L1025
@BvB93 It makes a lot more sense now, there is a branch that works only for `float32/64`.

https://github.com/numpy/numpy/blob/a33a10a5f494057b9ddcb9efe650306fccf8eba3/numpy/core/src/multiarray/multiarraymodule.c#L1010-L1016

Could this be fixed just by simply copying memory to `out` in case when it is specified?

Commit message: MAINT: fixed test for miscellaneous where non-broadcastable output operand is passed

Patch to review:
diff --git a/numpy/core/src/multiarray/multiarraymodule.c b/numpy/core/src/multiarray/multiarraymodule.c
index adc1558da7db..fc88746e0b71 100644
--- a/numpy/core/src/multiarray/multiarraymodule.c
+++ b/numpy/core/src/multiarray/multiarraymodule.c
@@ -1042,9 +1042,7 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
 #endif
 
     if (PyArray_NDIM(ap1) == 0 || PyArray_NDIM(ap2) == 0) {
-        result = (PyArray_NDIM(ap1) == 0 ? ap1 : ap2);
-        result = (PyArrayObject *)Py_TYPE(result)->tp_as_number->nb_multiply(
-                                        (PyObject *)ap1, (PyObject *)ap2);
+        result = (PyArrayObject *)PyObject_CallFunctionObjArgs(n_ops.multiply, ap1, ap2, out, NULL);
         Py_DECREF(ap1);
         Py_DECREF(ap2);
         return (PyObject *)result;
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index 3f3c1de8ea20..4a064827dd0f 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -6675,6 +6675,22 @@ def test_dot_3args_errors(self):
         r = np.empty((1024, 32), dtype=int)
         assert_raises(ValueError, dot, f, v, r)
 
+    def test_dot_out_result(self):
+        x = np.ones((), dtype=np.float16)
+        y = np.ones((5,), dtype=np.float16)
+        z = np.zeros((5,), dtype=np.float16)
+        res = x.dot(y, out=z)
+        assert np.array_equal(res, y)
+        assert np.array_equal(z, y)
+
+    def test_dot_out_aliasing(self):
+        x = np.ones((), dtype=np.float16)
+        y = np.ones((5,), dtype=np.float16)
+        z = np.zeros((5,), dtype=np.float16)
+        res = x.dot(y, out=z)
+        z[0] = 2
+        assert np.array_equal(res, z)
+
     def test_dot_array_order(self):
         a = np.array([[1, 2], [3, 4]], order='C')
         b = np.array([[1, 2], [3, 4]], order='F')
diff --git a/numpy/typing/tests/data/pass/ndarray_misc.py b/numpy/typing/tests/data/pass/ndarray_misc.py
index 19a1af9e292c..6beacc5d7cc9 100644
--- a/numpy/typing/tests/data/pass/ndarray_misc.py
+++ b/numpy/typing/tests/data/pass/ndarray_misc.py
@@ -150,7 +150,7 @@ class SubClass(np.ndarray): ...
 A.diagonal()
 
 A.dot(1)
-A.dot(1, out=B0)
+A.dot(1, out=B2)
 
 A.nonzero()
 


Code review comments:
diff_hunk: @@ -1042,9 +1042,7 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
 #endif
 
     if (PyArray_NDIM(ap1) == 0 || PyArray_NDIM(ap2) == 0) {
-        result = (PyArray_NDIM(ap1) == 0 ? ap1 : ap2);
-        result = (PyArrayObject *)Py_TYPE(result)->tp_as_number->nb_multiply(
-                                        (PyObject *)ap1, (PyObject *)ap2);
+        result = (PyArrayObject *)PyObject_CallFunctionObjArgs(n_ops.multiply, ap1, ap2, out, NULL);
         Py_DECREF(ap1);
         Py_DECREF(ap2);
         return (PyObject *)result;
path: numpy/core/src/multiarray/multiarraymodule.c
text: @user1:
```suggestion
        return mul_res;
```

diff_hunk: @@ -1042,9 +1042,7 @@ PyArray_MatrixProduct2(PyObject *op1, PyObject *op2, PyArrayObject* out)
 #endif
 
     if (PyArray_NDIM(ap1) == 0 || PyArray_NDIM(ap2) == 0) {
-        result = (PyArray_NDIM(ap1) == 0 ? ap1 : ap2);
-        result = (PyArrayObject *)Py_TYPE(result)->tp_as_number->nb_multiply(
-                                        (PyObject *)ap1, (PyObject *)ap2);
+        result = (PyArrayObject *)PyObject_CallFunctionObjArgs(n_ops.multiply, ap1, ap2, out, NULL);
path: numpy/core/src/multiarray/multiarraymodule.c
text: @user1:
```suggestion
        PyObject *mul_res = PyObject_CallFunctionObjArgs(
                n_ops.multiply, ap1, ap2, out, NULL);
```

@user1:
Just a nit, but I find it awkward that we cast back to array (and it may not even be true in corner cases).  This was already the case before of course.