ENH: Support the dot family for user-dtypes - #31574

np.dot must make out= as C-style array, so I think doing ascontiguousarray might be OK?
Ideally the matmul/BLAS should handle any layout (for perfomance cases we can revert back to let dtype's BLAS handle this)

If we revert then, this will require new release from numpy-quaddtype (PR numpy/numpy-quaddtype#108 is already landed)

By default, any ufunc produces a C-contiguous array, so I don't think you'd need to do anything.

Here, you could just do return np.matmul(a, b, out=out) - by default, output is contiguous.

By default, any ufunc produces a C-contiguous array, so I don't think you'd need to do anything.

Are you sure the detailed handling of the axes is necessary? Remember that with np.matmul, like with all gufuncs, you have the option of explicitly passing in the axes that are contracted; see https://numpy.org/doc/stable/reference/ufuncs.html#optional-keyword-arguments

Unfortunately, for ndim > 2 dot has a funny behavior that the result shape is a.shape + b.shape minus the reduction axes.
I.e. it doesn't actually do a proper broadcasting/stacked operation, it is the same as tensordot(..., axes=[-1, -2]).

However, tensordot actually does this precise dance already. So we could consider being a bit annoying and just rejecting it for this path (we should deprecate it on the other path I guess) with a note of "please use tensordot instead".

There is another good argument for this, the code has:

#if defined(HAVE_CBLAS)
    if (PyArray_NDIM(ap1) <= 2 && PyArray_NDIM(ap2) <= 2 &&

so we don't even use BLAS for this case, which is ridiculous and unless arrays are small users are better of to use the clearer tensordot() anyway.

(If we do that this might be simple enough to just move it to C.)

fallback's _dot_contract is the tensordot-style reshape=>matmul, faster than legacy N-D dot, so rejecting it isn't a perf win, it's purely about not carrying the deprecation-worthy N-D dot semantics. I'm fine either way.

Not sure I am following precisely, but I am leaning to just reject it with a "use tensordot" message. Let's ping @mhvk, though.
(That would mean we should aim for a deprecation as a tendency, I guess.)

Did you check all these ascontiguousarray calls are in fact necessary? I rather doubt it, since matmul produces a contiguous array.

I think you need vecdot here to ensure complex numbers are properly conjugated (do add at test! also for dot, where complex conjugation is not done [incorrectly, but the fallback needs to mimic that mistake...])

You can get the shape check by setting strict=True in assert_array_equal.

Also, I'm actually not 100% sure we should care about contiguity, but if we do, we should test it to ensure we don't regress, by adding something like assert got.strides == ref.strides (same in other tests).

Did you mean to put MIN_VERSION to 1.1 since I think we are at 1.0 (and we need to reject it).
I guess we can make it customizable if needed. The current min-version is vital because even importing is wrong, but after 1.1 is out that isn't the case anymore.

Yeah that also works, just weird habit of having even jumps for sub-versions :)

I think this comment still needs to be addressed...

I guess if the next release is 1.2 then it doesn't matter. Although, I don't really see why numpy-quaddtype should use a different version scheme (distinguising odd from even versions) than NumPy.

This can be a sliding_window_view I think. I think matmul will be OK with core dimensions that have self overlap (if not, we could pass axes to make it a non-core dimension). Although that is of course not BLAS compatible (for our code, I think it might even go back into einsum).

Otherwise, you would bloat memory usage like crazy here. (I suspect we can do this, but I am slightly unsure and i tmay make sense to split the PR into two, since I think dot is a bit clearer.)

I am too concerned about the memory bloat here + sliding_window_view user-dtypes don't support.

Right to split this, I will do another PR for correlate and with the current and we can make further changes there

sliding_window_view user-dtypes don't support.

Ah, that is true and annoying :/. Could possibly not go via __array_interface__ there to make it work, but I am not sure that is better than trying a fix for __array_interface__ (heck even adding a dtype field).

I think this comment still needs to be addressed...

I guess we have a bit of a chicken-egg problem in that this release doesn't exist yet. Is there any chance you can make the test work (with xfails maybe?) against 1.0 so can be pretty confident it works? Right now this whole test is dead code.

Right now v1.0 of numpy-quaddtype has the bug of registering broad promoter/loops and importing it can interfere with other unrelated tests (those issues are fixed upstream and planned to be out in v1.1) so technically we cannot even import it in the current released version.

So this guard is technically saving us from multiple issues (indeed a bad practice) hence I will make sure that next release of quaddtype rightly solves such problems

OK, then what do you think about deleting these dead tests in this PR and then re-adding them after you do a numpy-quaddtype release and bump the pin here?

Yeah good to me.

I dunno, I honestly prefer putting this in with the dead tests (that work and were tested locally), then with a bunch of unit-tests for the _dot_fallback().

(I don't like such unit-tests if it's avoidable. For one they don't age well if we delete _dot_fallback() and then the person who does that work has to now do the double work of also fixing all tests.)

Of course that relies on actually getting that quaddtype release out very soon (and this bug feels like almost enough of a reason to do so?)

should this specify order="C"? Otherwise stuff like this happens:

import numpy as np
from numpy._core.numeric import _dot_fallback

a = np.ones((2, 3, 4, 5)).swapaxes(0, 1)
b = np.ones((5, 6))

expected = np.dot(a, b)
actual = _dot_fallback(a, b)

print(expected.flags.c_contiguous, actual.flags.c_contiguous)
# True False

Matmul and the ufunc machinery have different out compatibility rules. I'm not sure whether this is worth worrying about, but there's also this incosistency:

import numpy as np
from numpy._core.numeric import _dot_fallback

a, b = np.ones((2, 3)), np.ones((3, 4))
outputs = {
    "extra batch": np.empty((5, 2, 4)),
    "Fortran": np.empty((2, 4), order="F"),
    "float32": np.empty((2, 4), dtype=np.float32),
    "strided": np.empty((2, 8))[:, ::2],
}

for label, out in outputs.items():
    for fn in (np.dot, _dot_fallback):
        try:
            result = fn(a, b, out=out)
            print(label, fn.__name__, "accepted", result.shape, result.dtype)
        except ValueError:
            print(label, fn.__name__, "ValueError")

On this PR I get:

extra batch dot ValueError
extra batch _dot_fallback accepted (5, 2, 4) float64
Fortran dot ValueError
Fortran _dot_fallback accepted (2, 4) float64
float32 dot ValueError
float32 _dot_fallback accepted (2, 4) float32
strided dot ValueError
strided _dot_fallback accepted (2, 4) float64

e.g. the "real" dot implementation rejects all these out shapes, dtypes, C/F order, and strides

Fixed, the fallback now passes order="C".
dot c_contiguous True | fallback c_contiguous True

and regarding dot's out= contract, earlier silently depends on the dtype. The fallback now validates out itself (ndim, C-contiguity, dtype) and reuses dot's own error message, so all four of your cases match:

extra batch  dot=ValueError  fallback=ValueError
Fortran      dot=ValueError  fallback=ValueError
float32      dot=ValueError  fallback=ValueError
strided      dot=ValueError  fallback=ValueError

I checked this with other possible permutation/combinations and all matches with the usual dot behaviour

It's not just np.dot though. Lots of other code paths go through MatrixProduct2. Maybe you could say "behaviour of numpy.dot, numpy.inner, numpy.tensordot, and other code paths like numpy.cov that rely on dot products"

Fixed

This still isn't quite right: it should also flag an incorrect output shape with the same number of dimensions and an un-aligned but C-contiguous output.

Rather than duplicating the rules in Python, can you instead extract the output validation from new_array_for_sum into a private C helper that this function can call?

Do you have an opinion on this?

I feel like we either give up and pass out as allocated and checked by dot(). Or we don't worry about any of this and just use the gufunc machinery as is (as it is the "gold standard" and what dot could transition to as well).

Here, we actually check that out.dtype is identical to the result_type, but we don't enforce it in the ufunc when out=None, where it isn't strictly guaranteed. At that point we are almost at casting="no" into the ufunc, but also that I would prefer not to.

(I feel like all of this is churn for fairly weird corner cases -- yes, I don't like the divergence much either, but...)

FWIW, when it comes to the order="C" that seems reasonable, because it is also a reasonable idea if we want to rewire things (unfortunately!).

In general, I don't want us to start to guarantee order, but maybe dot is special enough (to at least not rush changing it).

Do you have an opinion on this?

I feel like we either give up and pass out as allocated and checked by dot(). Or we don't worry about any of this and just use the gufunc machinery as is (as it is the "gold standard" and what dot could transition to as well).

Here, we actually check that out.dtype is identical to the result_type, but we don't enforce it in the ufunc when out=None, where it isn't strictly guaranteed. At that point we are almost at casting="no" into the ufunc, but also that I would prefer not to.

(I feel like all of this is churn for fairly weird corner cases -- yes, I don't like the divergence much either, but...)

I dunno, I honestly prefer putting this in with the dead tests (that work and were tested locally), then with a bunch of unit-tests for the _dot_fallback().

(I don't like such unit-tests if it's avoidable. For one they don't age well if we delete _dot_fallback() and then the person who does that work has to now do the double work of also fixing all tests.)

Of course that relies on actually getting that quaddtype release out very soon (and this bug feels like almost enough of a reason to do so?)

FWIW, when it comes to the order="C" that seems reasonable, because it is also a reasonable idea if we want to rewire things (unfortunately!).

In general, I don't want us to start to guarantee order, but maybe dot is special enough (to at least not rush changing it).