Fix utils.matmul

[r-pascua]

It turns out that there's a silent failure mode in the jitted version of utils.matmul. Apparently, numba is completely fine with running something that looks like this:

@numba.njit
def test(left, right):
    out = np.zeros_like(right)
    for i in range(out.shape[0]):
        for j in range(out.shape[1]):
            out[i,j] = left[i,j] * right[i,j]
    return out

# Try taking the Hadamard product of arrays with different shapes
left = np.random.normal(size=(5,3))
right = np.random.normal(size=(5,5))
out = test(left, right)

This code block surprisingly doesn't error. At a glance, the output looks somewhat normal, but a closer look will reveal that the output array has garbage values like 2.17e-320. This is a relatively common error (in my experience) when wrapping C code with Python and providing the C code with arrays that have the wrong data type, and unfortunately it can be pretty tricky to catch if you're not paying close attention to the code outputs.

Anyway, currently the utils.matmul function incorrectly assumes that providing arrays with bad shapes will error, so a sanity check should be performed prior to running the jitted matrix multiplication routines to catch this kind of edge case.