Implement `gemv` numba dispatch

[jessegrabowski]

Description

When working on #1416 I found that numba has terrible performance on matrix-vector multiplication:

@ricardoV94 thinks this is because numba probably only uses gemm for everything, and never uses gemv. Since we have a GEMV Op already, it will be very easy to follow the pattern I used in #1416 to make a dispatch for GEMV.

We might also have to add the GEMV rewrite to the numba mode as well -- I know they're disabled for jax, for example. Haven't checked for numba, but something to be aware of.

[ricardoV94]

Yes all BLAS rewrites are disabled on NUMBA, those are the ones that add things like GEMM, GEMV, GER

[ricardoV94]

Here is a MRE:

import numpy as np
import pytensor
import pytensor.tensor as pt

A = pt.matrix("A")
x = pt.vector("x")
out = A @ x

c_fn = pytensor.function([A, x], out, trust_input=True)
c_fn.dprint()
# CGemv{inplace} [id A] 2
#  ├─ AllocEmpty{dtype='float64'} [id B] 1
#  │  └─ Shape_i{0} [id C] 0
#  │     └─ A [id D]
#  ├─ 1.0 [id E]
#  ├─ A [id D]
#  ├─ x [id F]
#  └─ 0.0 [id G]

numba_fn = pytensor.function([A, x], out, mode="NUMBA", trust_input=True)
numba_fn.dprint()
# Squeeze{axis=1} [id A] 2
#  └─ dot [id B] 1
#     ├─ A [id C]
#     └─ ExpandDims{axis=1} [id D] 0
#        └─ x [id E]

rng = np.random.default_rng(1)
A_test = rng.normal(size=(1024, 512))
x_test = rng.normal(size=(512,))
np.testing.assert_allclose(c_fn(A_test, x_test), numba_fn(A_test, x_test))
%timeit c_fn(A_test, x_test) # 338 μs ± 8.23 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
%timeit numba_fn(A_test, x_test)  # 6.6 ms ± 978 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)