where can I find the numpy.matmul() source code?

I do not obtain the same results when I use np.matmul(A, b) in Python and when I use xtensor-blas's xt::linalg::dot(A, b) in C++.

I am investigating the reasons, as when saved and read from disk, A and b are identical when doing np.allclose(A, b) in Python.

The result of these 2 multiplications (in Py and in C++) is a 250 element 1D array. The first elements are identical, the elements from the middle are different, the elements from the end are identical. As below:

 [ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False  True  True False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True  True]

I identified what exactly is called by xtensor-blas: it is gemv from BLAS.

I want to see what is called by np.matmul(A, b). Does anyone know how to find what exactly is called from BLAS when doing np.matmul(A, b)?

Thank you!

This seems to be the function in question, so the CBLAS functions called seem to be cblas_?syrk (in special cases) and cblas_?gemm (in general).