Repository: numpy/numpy
Pull request: BUG: Fix float16 einsum fastpaths using wrong tempvar
Fix for gh-20305 @Qiyu8 maybe you know where to best add some tests?  I was hoping to add some parametrization and be done with it, but didn't really find the right test to parametrize for float16 that would fail without the fixes.

EDIT: In case anyone ever stumbles on this.  The previous tests probably did run all the code paths, the problem is the missing fractional part.  Using the wrong tempvar here, meant the accumulation was done in `(u)int16`.  Without a fractional part, that ends up with the correct result though.

Problem statement:
BUG: einsum all zero for float16
### Describe the issue:

Starting with numpy-1.21.0 (at least some) einsum operations for arrays with dtype=np.float16 yield all zero results.

### Reproduce the code example:

```python
import numpy as np

a = (np.linspace(0, 301, 300).reshape(10, 10, 3) / 301).astype(np.float16)
mat = np.array([[0.39, 0.54, 0.0089], [0.07, 0.96, 0.001], [0.02, 0.13, 0.94]], dtype=np.float16)
assert np.all((a @ mat.T) != 0)
assert np.all(np.einsum("ij, ...j", mat.astype(np.float32), a.astype(np.float32)) != 0)
assert np.all(np.einsum("ij, ...j", mat, a) != 0)
```


### Error message:

_No response_

### NumPy/Python version information:

1.21.2 3.8.12 (default, Oct 12 2021, 13:49:34) 
[GCC 7.5.0]

Hints from issue discussion:
@Qiyu8 could  you have a look at this?  This seems possibly related to the einsum SIMD refactor. 
Here is a more minimal reproducer (executed on 1.21.2):

```python
>>> import numpy as np
>>> a = np.array([1.764 , 0.4001], dtype=np.float16)
>>> np.einsum("i,i->", a, a)
1.0
>>> (a ** 2).sum()
3.271
```

the same code works fine on numpy 1.20.3.

I ran the above code on macos/arm64 with numpy from conda-forge.
Turns out the `temptype` usage is spotty in the new code, which would affect only float16.   Have to look into the tests it feels like they should have noticed at least some of this...

Commit message: BUG: Fix float16 einsum fastpaths using wrong tempvar

Patch to review:
diff --git a/numpy/core/src/multiarray/einsum_sumprod.c.src b/numpy/core/src/multiarray/einsum_sumprod.c.src
index 29ceabd71d58..3114a58960ef 100644
--- a/numpy/core/src/multiarray/einsum_sumprod.c.src
+++ b/numpy/core/src/multiarray/einsum_sumprod.c.src
@@ -337,13 +337,13 @@ static NPY_GCC_OPT_3 void
         /**begin repeat2
          * #i = 0, 1, 2, 3#
          */
-        const @type@ b@i@ = @from@(data[@i@]);
-        const @type@ c@i@ = @from@(data_out[@i@]);
+        const @temptype@ b@i@ = @from@(data[@i@]);
+        const @temptype@ c@i@ = @from@(data_out[@i@]);
         /**end repeat2**/
         /**begin repeat2
          * #i = 0, 1, 2, 3#
          */
-        const @type@ abc@i@ = scalar * b@i@ + c@i@;
+        const @temptype@ abc@i@ = scalar * b@i@ + c@i@;
         /**end repeat2**/
         /**begin repeat2
          * #i = 0, 1, 2, 3#
@@ -353,8 +353,8 @@ static NPY_GCC_OPT_3 void
     }
 #endif // !NPY_DISABLE_OPTIMIZATION
     for (; count > 0; --count, ++data, ++data_out) {
-        const @type@ b = @from@(*data);
-        const @type@ c = @from@(*data_out);
+        const @temptype@ b = @from@(*data);
+        const @temptype@ c = @from@(*data_out);
         *data_out = @to@(scalar * b + c);
     }
 #endif // NPYV check for @type@
@@ -417,14 +417,14 @@ static void
         /**begin repeat2
          * #i = 0, 1, 2, 3#
          */
-        const @type@ a@i@ = @from@(data0[@i@]);
-        const @type@ b@i@ = @from@(data1[@i@]);
-        const @type@ c@i@ = @from@(data_out[@i@]);
+        const @temptype@ a@i@ = @from@(data0[@i@]);
+        const @temptype@ b@i@ = @from@(data1[@i@]);
+        const @temptype@ c@i@ = @from@(data_out[@i@]);
         /**end repeat2**/
         /**begin repeat2
          * #i = 0, 1, 2, 3#
          */
-        const @type@ abc@i@ = a@i@ * b@i@ + c@i@;
+        const @temptype@ abc@i@ = a@i@ * b@i@ + c@i@;
         /**end repeat2**/
         /**begin repeat2
          * #i = 0, 1, 2, 3#
@@ -434,9 +434,9 @@ static void
     }
 #endif // !NPY_DISABLE_OPTIMIZATION
     for (; count > 0; --count, ++data0, ++data1, ++data_out) {
-        const @type@ a = @from@(*data0);
-        const @type@ b = @from@(*data1);
-        const @type@ c = @from@(*data_out);
+        const @temptype@ a = @from@(*data0);
+        const @temptype@ b = @from@(*data1);
+        const @temptype@ c = @from@(*data_out);
         *data_out = @to@(a * b + c);
     }
 #endif // NPYV check for @type@
@@ -521,14 +521,14 @@ static NPY_GCC_OPT_3 void
         /**begin repeat2
          * #i = 0, 1, 2, 3#
          */
-        const @type@ ab@i@ = @from@(data0[@i@]) * @from@(data1[@i@]);
+        const @temptype@ ab@i@ = @from@(data0[@i@]) * @from@(data1[@i@]);
         /**end repeat2**/
         accum += ab0 + ab1 + ab2 + ab3;
     }
 #endif // !NPY_DISABLE_OPTIMIZATION
     for (; count > 0; --count, ++data0, ++data1) {
-        const @type@ a = @from@(*data0);
-        const @type@ b = @from@(*data1);
+        const @temptype@ a = @from@(*data0);
+        const @temptype@ b = @from@(*data1);
         accum += a * b;
     }
 #endif // NPYV check for @type@


Code review comments:
diff_hunk: @@ -744,6 +746,53 @@ def test_einsum_all_contig_non_contig_output(self):
         np.einsum('ij,jk->ik', x, x, out=out)
         assert_array_equal(out.base, correct_base)
 
+    @pytest.mark.parametrize("dtype",
+             np.typecodes["AllFloat"] + np.typecodes["AllInteger"])
+    def test_different_paths(self, dtype):
+        # Test originally added to cover broken float16 path: gh-20305
+        # Likely most are covered elsewhere, at least partially.
+        # Note that it is plausible that this test needs to be relaxed for
+        # floating point round-off errors.
+        dtype = np.dtype(dtype)
+        # Simple test, designed to excersize most specialized code paths,
+        # note the +0.1 for floats.
+        arr = (np.arange(7) + 0.1).astype(dtype)
+        scalar = np.array(2, dtype=dtype)
+
+        # contig -> scalar:
+        res = np.einsum('i->', arr)
+        assert_allclose(res, dtype.type(21.7))
+        # contig, contig -> contig:
+        res = np.einsum('i,i->i', arr, arr)
+        assert_array_equal(res, arr * arr)
+        # noncontig, noncontig -> contig:
+        res = np.einsum('i,i->i', arr.repeat(2)[::2], arr.repeat(2)[::2])
+        assert_array_equal(res, arr * arr)
+        # contig + contig -> scalar
+        assert np.einsum('i,i->', arr, arr) == (arr * arr).sum()
path: numpy/core/tests/test_einsum.py
text: @user1:
This is failing with `AssertionError: assert 95.270004 == 95.27`

@user1:
Also with `AssertionError: assert 95.27 == 95.26999999999998`''

The problems are on arm64 and s390x.

@author:
arg, yeah it is a bit annoying: To trigger the issue I have to use non-integers but these rounding issues are tricky with pairwise summations.

But: I just realized using `+0.5` should work, should give exact results but trigger the issue well enough.