我试图从Cython使用Scipy的BLAS函数,我想在2D类型的内存视图上调用dnrm2的l2规范。根据BLAS文档,dnrm2返回向量的欧几里德范数。它等价于矩阵的Frobenius范数。我们将矩阵平坦化,然后传递给函数。
from scipy.linalg.cython_blas cimport dnrm2
ctypedef np.float32_t dtype_t
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef dtype_t l2_norm(dtype_t[:, ::1] array):
cdef int m = array.shape[0]
cdef int n = array.shape[1]
cdef np.intp_t lda = m
cdef np.intp_t incx = 1
cdef dtype_t result
cdef np.float64_t *data = <np.float64_t*>&array[0, 0]
result = <dtype_t>dnrm2(&n, data, &incx)
return result
dnrm2只接受fp64数字,所以我必须使用默认数据类型进行类型转换。
不幸的是,函数返回错误的结果。是我把数组平展错了,还是我用错了?cdef np.float64_t *data = <np.float64_t*>&array[0, 0]
对于以下矩阵
array([[0.65411323, 0.73329186, 0.74279535],
[0.617243 , 0.6950498 , 0.99457264],
[0.5428702 , 0.34277368, 0.73931944]], dtype=float32)
返回0.007486129179596901
而实际值为2.0807159881759274
cdef np.float64_t *data = <np.float64_t*>&array[0, 0]
这行不通。你有一个指向32位数组的指针。将其转换为指向64位数组的指针。这不会改变底层数据。它只是在它指向的地方撒谎。
我怀疑你想要snrm2
代替。如果不是,那么你需要as用不同的dtype(例如np.ndarray.astype
)复制底层数组