多次元配列の内積を計算するプログラムの高速化を検討しております。

Pythonで作成したソースをもとにCythonでも実装しました。
私の環境では両者の速度は以下の通りです。

Python : 0.17599 [sec]
Cython : 0.17589 [sec]

Cythonにより速くなることを期待しておりましたが殆ど変わりませんでした。

Cythonについては初心者です。コードにおかしなところ等ありますでしょうか?
PythonおよびCythonのソースコードを以下に示します。
なお、恐縮ですが、少なくとも10倍以上の処理速度向上を望んでおります。

Python

$ python cal_area_from_dem.py

cal_area_from_dem.py

import numpy as np
import time as t

start = t.time()
np.random.seed(0)
xyz = np.random.rand(4, 100000, 3)
vec = np.random.rand(3, 100000)

def main():

    x = np.array(xyz[:, :, 0])
    y = np.array(xyz[:, :, 1])
    z = np.array(xyz[:, :, 2])

    cvec = vec[:, :]

    p0 = np.array([x.T, y.T, z.T])
    p0 = [p0[:, i, :] for i in range(100000)]
    p_ref = np.array([x[0], y[0], z[0]])
    d_vec = p0 - np.repeat(p_ref, 4).reshape(100000, 3, 4)

    dat0 = np.einsum("ij,ijk->ik", cvec.T, d_vec)

    print('Python : ' + str(round((t.time() - start),5)) + ' [sec]')

if __name__ == "__main__":
    main()

Cython

$ python setup.py build_ext --inplace; python main_cal_area_from_dem.py

main_cal_area_from_dem.py

import cal_area_from_dem as c
import numpy as np
import time as t

def main():
    start = t.time()
    np.random.seed(0)
    xyz = np.random.rand(4, 100000, 3)
    vec = np.random.rand(3, 100000)
    test = c.cal_area_from_dem(xyz, vec)
    print('Cython : ' + str(round((t.time() - start),5)) + ' [sec]')

if __name__ == "__main__":
    main()

cal_area_from_dem.pyx

import cython
import numpy as np
cimport numpy as np

ctypedef np.float64_t np_float_t
ctypedef np.int32_t np_int_t

@cython.boundscheck(False)
@cython.wraparound(False)

cpdef np.ndarray[np_int_t, ndim=1] cal_area_from_dem(
    np.ndarray[np_float_t, ndim=3] xyz,
    np.ndarray[np_float_t, ndim=2] vec):

    cdef np.ndarray[np_float_t, ndim=2] x, y, z, cvec, p_ref, new_dat0
    cdef np.ndarray[np_float_t, ndim=3] p0, d_vec

    x = np.array(xyz[:, :, 0])
    y = np.array(xyz[:, :, 1])
    z = np.array(xyz[:, :, 2])

    cvec = vec[:, :]

    p0 = np.array([x.T, y.T, z.T])
    p0 = np.array([p0[:, i, :] for i in range(100000)])
    p_ref = np.array([x[0], y[0], z[0]])
    d_vec = p0 - np.repeat(p_ref, 4).reshape(100000, 3, 4)

    new_dat0 = np.einsum("ij,ijk->ik", cvec.T, d_vec)

setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
    cmdclass = {'build_ext': build_ext},
    ext_modules = [Extension("cal_area_from_dem", ["cal_area_from_dem.pyx"])]
)