The issue is that np.dot(a,b) for multidimensional arrays makes the dot product of the last dimension of a with the second-to-last dimension of b:

np.dot(a,b) == np.tensordot(a, b, axes=([-1],[2]))

As you see, it does not work as a matrix multiplication for multidimensional arrays. Using np.tensordot() allows you to control in which axes from each input you want to perform the dot product. For example, to get the same result in c_mat1 you can do: