Refactored to avoid duplication of code: matvec now uses the matmul_noblas loops internally (as does the matrix-vector implementation in matmul itself), and vecmat creates a loop that uses the loops in vecdot, letting their outer loop deal with the matrix dimension.

@charris - since you were the one originally asking for it, what do you think of this implementation? It is now quite lightweight, since it reuses inner loops.

macos failure unrelated

Moving milestone to 2.1, if we get this in we should change it again, but at least I don't think it needs to be prioritized.

Yes, not critical for 2.0 for sure, but maybe more logical to clear the milestone rather than set it to 2.1 if there still is a chance. This PR has become fairly simple after my refactor.

Out of curiousity, what is the status of this PR?

I think it would still be nice to get this in! I just rebased to make sure things still pass. @charris - as the person who originally suggested to have matvec and vecmat, what do you think?

The circle-ci failure does not seem related:

/home/circleci/repo/doc/source/reference/arrays.promotion.rst:150: ERROR: Unexpected indentation.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:154: WARNING: Block quote ends without a blank line; unexpected unindent.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:189: ERROR: Unexpected indentation.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:191: WARNING: Block quote ends without a blank line; unexpected unindent.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:210: ERROR: Unexpected indentation.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:212: WARNING: Block quote ends without a blank line; unexpected unindent.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:31: ERROR: Too many autonumbered footnote references: only 0 corresponding footnotes available.
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:31: ERROR: Unknown target name: "hist-reasons".
/home/circleci/repo/doc/source/reference/arrays.promotion.rst:68: ERROR: Unknown target name: "nep50".
source/release/notes-towncrier.rst:119: WARNING: Inline literal start-string without end-string.
/home/circleci/repo/doc/source/release/notes-towncrier.rst:119: WARNING: Inline literal start-string without end-string.
source/release/notes-towncrier.rst:119: WARNING: Inline literal start-string without end-string.

EDIT: also seen in #27091 so indeed unrelated.

One thing that occurred to me, it might also make sense to put these into np.linalg. OTOH, since matmul and vecdot are in the main namespace, it is probably clearer there.

One thing that occurred to me, it might also make sense to put these into np.linalg.

I think the main namespace is good enough, for now at least. It might change if the array API adopts them.

@seberg - shall we get this in? I fixed the error I had made in the change-log entry so now also circle-ci passes.

I am happy to get it in, someone should maybe confirm they like it as well. Am travelling right now, so may just punt for another while before doing another pass and merging, sorry.

Did you compare performance of calls to matmul and vecmat/matvec? I would be very surprised if they were different.

There should be no performance difference - it uses the same code. The advantage of matvec and vecmat is that one can pass in stacks of vectors and stacks of matrices instead of being stuck to a single vector.

ping also @charris who originally suggested we have these functions...

We discussed this at a recent triage meeting and would like to put it in. It needs a final review.

@seberg - I hadn't thought about implementation differences for float32 at all...

I'm fairly sure I initially had vecmat quite similar to matvec here, for the blas call just ensuring the conjugate, but with a separate no-blas implementation that ensured conjugation of the vector (since one couldn't use the matmat one). Then, realizing that vecdot had more or less what was needed already, it just seemed simpler to use it, reducing the amount of extra code...

I do think it is better if the two functions behave the same way. My tendency would be to use the better precision sums for both, which would mean repeating the complex vecdot implementation without conjugation for matvec. That said, it is easier to imaging vector multiplication where the precision matters than matrix-vector, since the matrices would generally become quite huge.

I'll look a bit at whether one can do this fairly easily, but let me know if you have suggestions...

@seberg - after spending most of the day on this (grr... still do not quite get the blas dimension stuff - too much switching between fortran and C...), I ended up with a bit of a frankencode: now the blas matrix-vector gemv will be called for float arguments, but not for complex. In principle, it works fine for complex for matvec, but it cannot handle the conjugation of the vector in vecmat; for symmetry, I thought it might be best to just fall through to my default implementation in that case -- which uses the various *_dot, so can still be blas-optimized in some cases.

What do you think?

Sorry, I didn't mean to push you into a rabbit hole :/, just realized that there likely is a difference due to the discussion about sum and float32 and was wondering if we should worry.
I think you are right, gemm could be used to achieve this probably, but gemv doesn't seem to have the transposes we would need.

I suspect your frankencode is fine (and the original would probably also have been fine). @mattip should we just go ahead with this?

Pushed a commit to fix the __module__ thing that failed CI in main. Points at one other thing: I would be totally fine with only adding these to functions to numpy.linalg (no strong opinion, though).

EDIT: Of course feel free to force push my commits away.

Thanks for pushing the fixes - I did squash them together so it is one commit again.

Given that vecdot and matmul are in the main namespace, I think these might as well too...

Yeah, I think we can just put this in for 2.2. @charris wanted to branch this weekend. I pinged the mailing list; but if anyone shouts there (or here), we can always do a backport (it'll be a few weeks before the actual release anyway).

With a fresh look, I now can get gemm to work. Also noticed a slight omission in the blasable tests. So, wait for merging just a little!

OK, now did the gemm for vecmat complex - should be ready!

Thanks Marten.